Kick Risk Forecasting and Evaluating During Drilling Based on Autoregressive Integrated Moving Average Model

Abstract

Timely forecasting of the kick risk after a well kick can reduce the waiting time after well shut-in and provide more time for well killing operations. At present, the multiphase flow model is used to simulate and forecast the pit gain and casing pressure. Due to the complexity of downhole conditions, calculation of the multiphase flow model is difficult. In this paper, the time series analysis method is used to excavate the information contained in the time-varying data of pit gain and casing pressure. A forecasting model based on a time series analysis method of pit gain and casing pressure is established to forecast the pit gain and casing pressure after a kick. To divide the kick risk level and achieve the forecasting of the kick risk before and after well shut-in, kick risk analysis plates based on pit gain and casing pressure are established. Three pit gain cases and one casing pressure case are studied, and a comparison between measured data and predicted data shows that the proposed method has high prediction accuracy and repeatability.

1. Introduction

During drilling operation, a permanent concern is kick risk forecasting and well control. A kick is defined as the uncontrolled influx of formation fluid or gas into the wellbore, which occurs when the wellbore pressure at a certain location is less than the formation pressure, and requires a well control emergency plan [1]. Influx fluids tend to have a high pressure; a gas kick is much more detrimental than a liquid kick due to gas expansion and, as a result, has a higher variation in pressure. Late kick detection and risk forecast increase the amount of formation fluid that enters the well borehole, which increases the kick pressure and makes it hard to control the kick and can even lead to an uncontrolled blowout. For well control safety, it is necessary to forecast the risk of a kick to reduce the shut-in time of a well and provide more time for the killing operation.

As a kick can pose a significant threat to safety drilling, a lot of work has been devoted to detecting kicks. Nowadays, kick-detection methods have developed from simple thresholding methods [2] into complex intelligent methods based on machine learning [3,4,5,6,7,8,9,10]. Sun et al. [3] proposed one integrated pattern-recognition model consisting of a dynamic multiphase wellbore flow model and improved the piecewise approximation and similarity measure algorithms for kick diagnosis. P. Andia and R. R. Israel [4] utilized a cyber-physical approach that combined physics-based modeling with Bayesian mathematics for detecting a kick. Shi et al. [5] adopted the random forest and support vector machine algorithms to detect a kick during drilling in real-time. Alouhali et al. [6] trained and evaluated five models—Decision Tree, K-Nearest Neighbor (KNN), Sequential Minimal Optimization (SMO) Algorithm, Artificial Neural Network (ANN), and Bayesian Network—to detect kicks based on actual kick cases. The Decision Tree and K-Nearest Neighbor were shown to be slightly more accurate than the others. Hargreaves et al. [7] used a Bayesian probabilistic framework to detect a kick early based on noisy drilling data. Pournazari et al. [8] proposed a pattern-recognition system based on Symbolic Aggregate Approximation (SAX) to identify a kick. Unrau et al. [9] undertook modern machine learning techniques to reduce false kick alarm rates while increasing the probability of kick detection. Brakel J et al. [10] applied modern noise reduction and pattern-recognition signal processing to produce earlier and more informative kick-detection notification. However, their work was mainly done to detect a kick, whether or not by machine learning, not to predict the development trend of kick risk.

Once a kick occurs, the first thing to do is to decide whether and when to close the well. According to specifications for well control technology in oil and gas drilling by the Sinopec Group, an alarm should be triggered when the pit gain is 2 m³, and the well should be closed when the pit gain is 3 m³. As pit gain can indirectly reflect the volume and rate of formation fluid intrusion into the wellbore, pit gain can be treated as an important indicator of kick risk before well shut-in. Once the well is shut in, the casing annulus pressure will rise because of the wellbore afterflow and gas slippage, thus, casing pressure can be treated as the indicator of kick risk after well shut-in. Therefore, the kick risk can be forecasted through forecasting the pit gain before well shut-in and the casing pressure after well shut-in.

At present, as it was in the past, gas–liquid two-phase fluid dynamics is the main method used by scholars to simulate and forecast kick parameters before and after well shut-in. John Billingham [11] considered the influences of wellbore afterflow and slippage and simulated and analyzed the change in wellbore pressure during shut-in by using a continuous gas column model. Avelar et al. [12], Gao et al. [13], and Sun et al. [14] used gas–liquid two-phase flow to simulate the overflow and wellbore pressure distribution in deepwater drilling. Ren et al. [15] established a gas–liquid two-phase flow model to determine the gas invasion velocity and location during drilling by using the relationship between the standpipe pressure and the gas invasion velocity and height. Aarsnes et al. [16] and Meng et al. [17] used a transient two-phase flow model to simulate and predict the annulus pressure profile in the gas invasion into the wellbore. Based on percolation mechanics and gas–liquid two-phase flow, Ren et al. [18] considered the wellbore afterflow effect and gas slippage effect to simulate the bottom hole pressure during shut-in. Due to the laws of formation fluid invading the wellbore and those of the multiphase fluid flow in the wellbore being very complex when kick occurs, it is difficult to describe the physical process accurately using a mathematical model. Considering the difficulty and efficiency of calculation, the wellbore flow model can be simplified, resulting in a reduction in the accuracy of prediction. Besides, because of the complex downhole conditions, it is difficult to determine the intrusion rate of formation fluid and the distribution of gas–liquid two-phase in the wellbore accurately at the initial time. Empirical formulas or hypothetical methods are usually used, so the initial and boundary conditions do not reflect the actual annulus situation of the wellbore, thus reducing the prediction accuracy.

An alternative method is to construct a prediction model based on the parameters themselves and to use the already measured parameters to predict the subsequent parameters. For the prediction of short-term and medium-term time-varying parameters, the time series analysis method (e.g.; ARIMA) can reveal the development laws of the time series from the autocorrelation point of the parameters and has the advantage of having a high fitting degree [19,20]. Characteristic parameters, such as the pit gain and casing pressure, change regularly with time in a short period. Therefore, this method can extract the feature information in a short time, mine the valuable information hidden in the historical data in a short time, and construct the time series prediction model of pit gain and casing pressure to predict the pit gain and casing pressure after a kick. Time series analysis methods are widely used in the field of petroleum prediction, such as crude oil price fluctuations [21], drilling site risk probability prediction [22], oil and gas well production prediction [23,24], but they are not applied in the prediction of kick risk. Therefore, this paper uses the summation autoregressive moving average model (ARIMA) to predict pit gain before shut-in and casing pressure after shut-in and uses the predicted values of pit gain and casing pressure to analyze the kick risk before and after shut-in.

2. Mathematical Models

In the field of statistics, a set of random variables arranged in chronological order (X₁, X₂, …, X_t, …) is often used to represent the time series of a random event, denoted as {X_t, t$\in$T}. The n ordered observations of the series are represented by {x_t, t = 1, 2,..., n}, referred to as the series of observations of series length n. The purpose of studying time series is to reveal the properties of a random time series {X_t}. To get the properties of a random time series {X_t}, we need to analyze the properties of its observation series {x_t} and infer the properties of the random time series {X_t} from the properties of the observation series [25].

As the series of kick characterization parameters, such as the pit gain or casing pressure, are mostly a non-stationary time series, the autoregressive moving average model (ARIMA) is the most suitable [26]. The expression of ARIMA (p,d,q) is

$$\{\begin{array}{l}\Phi (B){\nabla }^d{x}_t=\Theta (B){\epsilon }_t\\ E({\epsilon }_t)=0,Var({\epsilon }_t)={\sigma }_{\epsilon }^2,E({\epsilon }_t{\epsilon }_s)=0,s\ne t\\ E{x}_s{\epsilon }_t=0,\forall s<t\end{array}$$

(1)

where $\left\{x_t\right\}$ is the observation series, ${\varphi }_p$ is the autoregressive coefficient, $\left\{{\epsilon }_t\right\}$ is the random interference series, $E$ is the mean value, $Var$ is the variance, $B$ is the delay operator, d is the differential order, $\Phi (B)=1-{\varphi }_{1}B-\cdots -{\varphi }_p{B}^p$ is the autoregressive coefficient polynomial for the stationary reversible ARMA (p,q) model, and $\Theta (B)=1-{\theta }_{1}B-\cdots -{\theta }_q{B}^q$ is the moving smoothing coefficient polynomial for the stationary reversible ARMA (p,q) model.

3. Methodology

The steps to establish the model are as follows (Figure 1):

A gas well named SC-X1 in Sichuan Basin encountered a gas kick during the drilling of the gas zone. Figure 2 is a time sequence diagram of the measured pit gain data before well shut-in. This is used to introduce how to use this methodology in combination with the theoretical model.

3.1. Stationarity Test of Parameter Series

There are three main methods that test the stability of the parameter series—the time series diagram test, the autocorrelation diagram test, and the unit root test [27]:

3.1.1. Time Series Diagram Test

For a stationary parametric time series, the series’ mean and variance are constant, so the sequence shown in the sequence diagram fluctuates around a constant value and the range of fluctuation is bounded. If the sequence presents an obvious trend and periodicity, the sequence is not stable.

3.1.2. Autocorrelation Diagram Test

The stationary series has a short-term correlation, that is, the autocorrelation coefficient (ACF), ${\widehat{\rho }}_k$, of the parameter sequence decays to zero quickly with the increase in the number of delay periods, k. However, for non-stationary series, the autocorrelation coefficient decays to zero faster than the stationary series:

$$\mathrm{ACF}:{\widehat{\rho }}_k=\frac{\widehat{\gamma }(k)}{\widehat{\gamma }(0)},\forall 0<k<n$$

(2)

where $\widehat{\gamma }(k)$ is the estimation of the delayed k-phase auto-covariance function and $\widehat{\gamma }(0)$ is the estimation of the population variance:

$$\widehat{\gamma }(k)=\frac{{\displaystyle \sum _{t=1}^{n-k}(x_t-\overline{x})(x_{t+k}-\overline{x})}}{n-k},\forall 0<k<n$$

(3)

$$\widehat{\gamma }(0)=\frac{{\displaystyle \sum _{t=1}^n{(x_t-\overline{x})}^2}}{n-1}$$

(4)

3.1.3. Unit Root Test

By calculating the augmented Dickey–Fuller (ADF) test statistic of the series, it can be compared with the critical value of the $\tau$ test statistic to calculate the probability of the critical value being less than the $\tau$ test statistic. A probability, P, of less than 0.05 indicates that the series is stable.

The hypotheses H₀: X_t is non-stationary and H₁: X_t is stationary and are tested with the regression equation $\nabla x_t=({\varphi }_1+\cdots +{\varphi }_p-1)x_{t-1}-({\varphi }_2+\cdots +{\varphi }_p)\nabla x_{t-1}-\cdots -{\varphi }_{p-1}\nabla x_{t-p+1}+{\epsilon }_t$. H₀ is accepted with a P-value > 0.05, otherwise H₁ is accepted.

From Figure 2 it can be seen that there is a marked downward trend in the sequence diagram. Figure 3 shows the autocorrelation coefficient and the partial autocorrelation coefficient of the pit gain time series, and it can be seen that the autocorrelation coefficient decreases slowly to zero and the autocorrelation coefficient map shows obvious triangular symmetry. The unit root test results of the pit gain time series are shown in

Table 1. Unit root test results of the pit gain time series.

		$\mathit{\tau }$ Statistics	Probability P
Test Statistics of ADF		3.49	1
Critical value of the test	1% Level	−3.81	- -
	5% Level	−3.02	- -
	10% Level	−2.65	- -

ADF is the augmented Dickey–Fuller.

, and the probability (P) of the ADF test statistic is greater than 0.05. Therefore, the pit gain series is not a stationary sequence and can be derived from Figure 2, Figure 3, and Table 1. To establish the time series prediction model, differential processing is needed.

3.2. Processing of Non-Stationary Parameter Series

For non-stationary parameter series, it is necessary to process and extract the effective information in the sequence, that is, the deterministic information contained in the series. Generally, the difference method is used to extract the deterministic information of non-stationary series [28].

Figure 4 shows the autocorrelation coefficient and partial autocorrelation coefficient of the pit gain time series after the first-order difference.

Table 2. Unit root test results of the pit gain time series after the first-order difference.

		$\mathit{\tau }$ Statistics	Probability P
Test Statistics of ADF		−2.32	0.17
Critical value of the test	1% Level	−3.78	- -
	5% Level	−3.01	- -
	10% Level	−2.64	- -

shows the unit root test results of the pit gain time series after the first-order difference. As the probability (P) of the ADF test statistic is greater than 0.05, the series is still not a stationary sequence after the first-order difference, and second-order differential processing is needed.

Figure 5 shows the autocorrelation coefficient and partial autocorrelation coefficient of the pit gain time series after the second-order difference; the second-order difference series has a strong short-term correlation.

Table 3. Unit root test results of the pit gain time series after the second-order difference.

		$\mathit{\tau }$ Statistics	Probability P
Test statistics of ADF		−4.54	0.0025
Critical value of the test	1% Level	−3.86	- -
	5% Level	−3.04	- -
	10% Level	−2.66	- -

shows the unit root test results of the pit gain time series after the second-order difference. The probability (P) of the ADF test statistic is less than 0.05. Therefore, the series is stable after the second-order difference, and the model can be established.

3.3. Model Identification and Order Determination

The values of the autocorrelation coefficient (ACF) and partial autocorrelation coefficient (PACF) are calculated to determine p and q, that is, to identify and order the model [29]:

$$\mathrm{PACF}:{\rho }_{x_t,x_{t-k}|x_{t-1},\cdots ,x_{t-k+1}}=\frac{E[(x_t-\widehat{E}{x}_t)(x_{t-k}-\widehat{E}{x}_{t-k})]}{E[{(x_{t-k}-\widehat{E}{x}_{t-k})}^2]}$$

(5)

where $\widehat{E}{x}_t=E[x_t|x_{t-1},\cdots ,x_{t-k+1}]$ and $\widehat{E}{x}_{t-k}=E[x_{t-k}|x_{t-1},\cdots ,x_{t-k+1}]$.

According to the autocorrelation and partial autocorrelation diagram of the second-order difference of the pit gain series, shown in Figure 5, the preliminary order of the model can be determined. It can be seen from the figure that the autocorrelation coefficients fall within two standard deviations after delaying the first-order difference and exhibit tailing. The partial autocorrelation coefficients are within two standard deviations after delaying the first-order difference and exhibit tailing. Therefore, it can be preliminarily determined that p = 4 or 5 and q = 1, and the model can be preliminarily determined to be either ARIMA (4,2,1) or ARIMA (5,2,1).

3.4. Model Validity Test

After the model is fitted, it is necessary to test or evaluate the effectiveness of the fit of the model. If the fitting effect does not meet the evaluation requirements, it needs to be remodeled. This mainly includes a residual analysis and significant analysis of the parameters [30]:

3.4.1. Residual Analysis

For a model with a good degree of fit, the useful information in the measurement sequence should be extracted. The fitting residual term should have no information related to the sequence, that is, its residual term should tend toward a white noise sequence.

The white noise test process is as follows: it is assumed that the sequence values of the delay period, less than or equal to the m period, are independent of each other, and the Q-statistic method is used for analysis:

$$Q=n{\displaystyle \sum _{k=1}^m{\widehat{\rho }}_k^2}$$

(6)

where n is the measurement sequence period and m is the delay period.

The Q-statistic approximately obeys the chi-square distribution with a degree of freedom of m:

$$Q=n{\displaystyle \sum _{k=1}^m{\widehat{\rho }}_k^2}\sim {\chi }^2(m)$$

(7)

It can be seen from Equation (7) that when the Q-statistic is greater than the quantile ${\chi }_{1-\alpha }{}^2(m)$ or if the P-statistic is greater than $\alpha$, then the hypothesis is accepted and the residual sequence is a white noise sequence.

3.4.2. Significance Test of Parameters

Through the significance test of the parameters, it can be tested whether the unknown parameter is significant to zero, which makes the fitted model more compact and more convenient for predictive calculation.

The test begins with the following assumptions:

$$H_0:{\beta }_j=0\leftrightarrow H_1:{\beta }_j\ne 0,\forall 1\le j\le m$$

(8)

$$E(\widehat{\beta })=E[{(X^{\prime }X)}^{-1}{X}^{\prime }\tilde{y}]={(X^{\prime }X)}^{-1}{X}^{\prime }X\tilde{\beta }=\tilde{\beta }$$

(9)

$$Var(\widehat{\beta })=Var[{(X^{\prime }X)}^{-1}{X}^{\prime }\tilde{y}]={(X^{\prime }X)}^{-1}{X}^{\prime }X{(X^{\prime }X)}^{-1}{\sigma }_{\epsilon }^2={(X^{\prime }X)}^{-1}{\sigma }_{\epsilon }^2$$

(10)

where $\tilde{\beta }$ is the least-squares estimation.

${\widehat{\beta }}_j$ obeys a normal distribution, i.e., ${\widehat{\beta }}_j\sim N(0,a_{jj}{\sigma }_{\epsilon }{}^2),1\le j\le m$, where ${\sigma }_{\epsilon }{}^2=\frac{Q(\tilde{\beta })}{n-m}$.

Thus, a test statistic (T-test statistic) can be constructed for testing the significance of unknown parameters:

$$T=\sqrt{n-m}\frac{{\widehat{\beta }}_j}{\sqrt{a_{jj}Q(\tilde{\beta })}}\sim t(n-m)$$

(11)

When $\left|T\right|\ge t_{1-\alpha }(n-m)$, the parameters are significant. The probability (P-value) of the test statistics can also be used for analysis. When P is less than $\alpha$, the parameters are significant.

Table 4. Results of the model parameter significance $\tau$ statistic test (P-value).

Model		C	AR(1)	AR(2)	AR(3)	AR(4)	AR(5)	MA(1)
ARIMA (4,2,1)	a	0.0001	0.4269	0.0954	0.1898	0.0092		0
	b	0.2333		0.1009	0.7544	0.0031		0.001
	c	0.0372			0.921	0.0073		0.001
	d	0.0151				0.0044		0.0003
	e	0.0001	0.1524		0.2888	0.0282		0.0001
	f	0.0004	0.2145			0.0062		0
	g	0.0004	0.6243	0.1266		0.0027		0
	h	0.0001		0.049		0.0009		0
ARIMA (5,2,1)	a	0.0001	0.8719	0.0311	0.5338	0.0106	0.2114	0
	b	0		0.0372	0.479	0.0067	0.1029	0
	c	0.1186	0		0.6073	0.0724	0.2626	0.0068
	d	0	0.6756	0.0325		0.004	0.0833	0
	e	0.6893	0.0486	0.4462	0.3219		0.8929	0.0002
	f	0.6508	0.0482	0.7756			0.2628	0.0002
	g	0			0.8892	0.0401	0.0279	0.0003
	h	0.343	0			0.0958	0.2713	0
	i	0.4383	0			0.1029	0.945	0.0072
	j	0.0007		0.2377	0.2332		0.055	0.0009
	k	0		0.0854		0.0035	0.0447	0
	l	0				0.0331	0.0224	0
	m	0.0001			0.6014		0.0166	0.0001
	n	0.0013		0.4441			0.0072	0.0003
	o	0.5813	0				0.1888	0
	p	0					0.0033	0

shows the results of the model parameter significance $\tau$ statistic test (P-value), and models with P-values less than 0.05 are ARIMA (4,2,1)_d, ARIMA (4,2,1)_h, ARIMA (5,2,1)_l, and ARIMA (5,2,1)_p, that is, all the parameters of the these four models are significant.

3.5. Model Optimization

Since more than one model passes the test, to select the optimal model, the model needs to be optimized and selected. The criteria for optimal selection are the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (SBC)—the optimal model being when the AIC and SBC of the fitted model are the smallest [31]:

$$AIC=-2\ln MLE+2NUP$$

(12)

$$SBC=-2\ln MLE+\ln NUP$$

(13)

where MLE is the maximum likelihood estimate and NUP is the number of unknown parameters in the model.

Table 5. White noise test results and Akaike Information Criterion (AIC) and Bayesian Information Criterion (SBC) values of the residual sequences.

Model	AIC	SBC	Delay Order	Q-Statistics	P-Value	Whether to Pass the Test
ARIMA (4,2,1)_d	−4.69	−4.55	6	16.11	0.003	No
			12	21.42	0.018
ARIMA (4,2,1)_h	−3.9	−3.7	6	7.91	0.05	Yes
			12	13.04	0.16
ARIMA (5,2,1)_l	−4	−3.81	6	10.01	0.018	No
			12	16.85	0.051
ARIMA (5,2,1)_p	−3.7	−3.6	6	6.68	0.15	Yes
			12	10.91	0.36

shows the white noise test results and the AIC and SBC values of the four models’ residual sequences. It can be seen that the P-values of the Q-statistics of the ARIMA (4,2,1)_h and ARIMA (5,2,1)_p models are greater than 0.05, and their residual sequences are white noise sequences, which satisfy the white noise test. Compared with the AIC and SBC of the ARIMA (4,2,1)_h model and the ARIMA (5,2,1)_p model, the AIC (−3.9) and SBC (−3.7) of ARIMA (4,2,1)_h model are smaller, so the ARIMA (4,2,1)_h model is the best model.

3.6. Parameter Estimation

The parameter estimation determines the model structure and order to establish the ARIMA (p,d,q) and then estimate the value of the unknown parameter in the model. In this paper, the least-squares estimation is used for the parameter estimation [32].

For the ARIMA (p,d,q) model,

$$\tilde{\beta }=({\varphi }_1,\cdots ,{\varphi }_p,{\theta }_1,\cdots ,{\theta }_q{)}^{\prime }$$

(14)

$$F_t(\tilde{\beta })={\varphi }_1{x}_{t-1}+\cdots +{\varphi }_p{x}_{t-p}-{\theta }_1{\epsilon }_{t-1}-\cdots -{\theta }_q{\epsilon }_{t-q}$$

(15)

Its residual term is

$${\epsilon }_t=x_t-F_t(\tilde{\beta })$$

(16)

The sum of the squared residuals is

$$Q(\tilde{\beta })={\displaystyle \sum _{t=1}^n{\epsilon }_t^2}={\displaystyle \sum _{t=1}^n{(x_t-{\varphi }_1{x}_{t-1}-\cdots -{\varphi }_p{x}_{t-p}+{\theta }_1{\epsilon }_{t-1}+\cdots +{\theta }_q{\epsilon }_{t-q})}^2}$$

(17)

The least-squares estimate of $\tilde{\beta }$ is a set of parameters that minimizes the sum of the squared residuals.

The model parameters are shown in

Table 6. Model parameters of ARIMA (4,2,1)_h.

Variable	C	AR(2)	AR(4)	MA(1)
Coefficient	0.0005	0.348	−0.684	−0.954

.

3.7. Forecasting

Future parameters are predicted based on the principle of the minimum mean square error [33]. The ARIMA (p,d,q) model can be represented by a linear function of random perturbation terms, namely,

$$x_t={\epsilon }_t+{\psi }_1{\epsilon }_{t-1}+{\psi }_2{\epsilon }_{t-2}+\cdots =\psi (B){\epsilon }_t$$

(18)

where ${\psi }_1,{\psi }_2,\cdots$ are determined by $\Phi (B){(1-B)}^d\psi (B)=\Theta (B)$.

${\psi }_1,{\psi }_2,\cdots$ satisfy the following formula:

$$\{\begin{array}{l}{\psi }_1={\varphi }_1-{\theta }_1\\ {\psi }_2={\varphi }_1{\psi }_1+{\varphi }_2-{\theta }_2\\ ⋮\\ {\psi }_j={\varphi }_1{\psi }_{j-1}+\cdots +{\varphi }_{p+d}{\psi }_{j-p-d}-{\theta }_j\end{array},\mathrm{where}{\psi }_j=\{\begin{array}{l}0,j<0\\ 1,j=0\end{array}$$

(19)

Thus, the true value of $x_{t+l}$ can be obtained with

$$x_{t+l}=({\epsilon }_{t+l}+{\psi }_1{\epsilon }_{t+l-1}+\cdots +{\psi }_{l-1}{\epsilon }_{t+1})+({\psi }_l{\epsilon }_t+{\psi }_{l+1}{\epsilon }_{t-1}+\cdots )$$

(20)

However, the parameters ${\epsilon }_{t+l}$ and ${\epsilon }_{t+l-1}$ cannot be directly obtained, so the estimated value of $x_{t+l}$ is

$${\widehat{x}}_t(l)={\psi }_0^{*}{\epsilon }_t+{\psi }_1^{*}{\epsilon }_{t-1}+{\psi }_2^{*}{\epsilon }_{t-2}+\cdots$$

(21)

The mean square error between the true and estimated values is

$$E{[x_{t+l}-{\widehat{x}}_t(l)]}^2=(1+{\psi }_1{}^2++{\psi }_{l-1}{}^2){\sigma }_{\epsilon }^2+{\displaystyle \sum _{j=0}^{\infty }{({\psi }_{l+j}-{\psi }_j{}^{*})}^2{\sigma }_{\epsilon }^2}$$

(22)

In the case of minimizing the mean square error (when ${\psi }_j^{*}={\psi }_{l+j}$), the estimated value of the l-period is

$${\widehat{x}}_t(l)={\psi }_l{\epsilon }_t+{\psi }_{l+1}{\epsilon }_{t-1}+{\psi }_{l+2}{\epsilon }_{t-2}+\cdots$$

(23)

The l-period estimation error is

$$e_t(l)={\epsilon }_{t+l}+{\psi }_1{\epsilon }_{t+l-1}+\cdots +{\psi }_{t-1}{\epsilon }_{t+1}$$

(24)

Therefore, the estimated value plus the estimated error is the true value of the prediction:

$$x_{t+l}=({\psi }_l{\epsilon }_t+{\psi }_{l+1}{\epsilon }_{t-1}+{\psi }_{l+2}{\epsilon }_{t-2}+\cdots )+({\epsilon }_{t+l}+{\psi }_1{\epsilon }_{t+l-1}+\cdots +{\psi }_{t-1}{\epsilon }_{t+1})$$

(25)

Figure 6 is the result of forecasting pit gain after 5 min using the ARIMA (4,2,1) model. It can be seen from the figure that the pit gain has not reached 2 m³ at 12 min, and the well has not been closed at this time. By predicting the pit gain after 5 min, the predicted pit gain has reached 3.18 m³ and it is also showing a gradually increasing trend.

3.8. Forecast Performance Measure

The mean relative error (MRE) is used for assessing the quality of the forecast. MRE can be considered as the mean of the relative values of the forecast errors and can be expressed with the following formula:

$$MRE=\frac{{\displaystyle {\sum }_{t=1}^N\left|y_t-{y^{\prime }}_t\right|/y_t}}{N}$$

(26)

where $y_t$ is the realized value at time period t, $y^{\prime }$ is the forecast about that period, and N is the number of data points.

3.9. Kick Severity Analysis

3.9.1. Analysis of the Kick Severity before Closing the Well

It can be seen from the permeation fluid mechanics and natural gas reservoir engineering that when drilling a gas-bearing reservoir the negative pressure difference at the bottom of the well has not yet reached the boundary of the formation, so the fluid flow in the formation can be regarded as the plane radial flow in an infinitely homogeneous isotropic gas reservoir. At this point, the rate at which natural gas invades the wellbore is [34]

$$Q_g=\frac{\left[p_i^{}-p_{wf}^{{}^{}}\left(t\right)\right]Kh{B}_g}{7.33\times {10}^{-3}\mu ZT\lg \left(8.085\eta t/r_w^2\right)}$$

(27)

where $Q_g$ is the gas flow rate from the formation into the wellbore in m³/d, $P_i$ is the supply pressure in MPa, $P_{wf}(t)$ is the bottom hole pressure in MPa, $K$ is the effective permeability in μm², $h$ is the effective thickness of the gas reservoir layer in m, $B_g$ is the volume coefficient of gas, $\mu$ is the gas viscosity in mPa·s, $Z$ is the natural gas compression factor, T is the formation temperature in K, $\eta$ is the pressure transmitting coefficient in μm²·MPa/(mPa·s), t is the invasion time in h, and $r_w$ is the borehole radius in m.

The amount of natural gas that invades the wellbore in t days is

$$V={\displaystyle \sum (Q\times t})$$

(28)

where $V$ is the amount of natural gas that has invaded the wellbore in t time in m³ and t is the invasion time in h.

From Equation (27) and Equation (28), it can be seen that the amount of gas intrusion is related to the reservoir parameters, bottom hole pressure, and time. The pit gain and the growth rate of the pit gain (the change rate of pit gain with time) can reflect the severity of a kick. In the curve of pit gain over time, several horizontal lines are used to distinguish the pit gain level and several straight lines with different slopes are used to distinguish the pit gain rate. Therefore, the kick severity analysis chart before closing the well includes several horizontal lines of pit gain and several oblique lines of the pit gain change rate. According to the well control standard and related operation specifications of the Chuandong area in China [35], the maximum pit gain should not exceed 3 m³ when drilling the reservoir and the well should be shut in when the overflow is 2 m³. Moreover, according to relevant statistics in more than 100 wells, more than two-thirds of the wells had a time interval from the discovery of overflow to the blowout of less than 30 min with more than half of the two-thirds being within 10 min [36]. The kick severity analysis chart based on well control standards and statistical data before closing the well is shown in Figure 7.

Figure 7 is divided into nine regions according to the pit gain increment and the kick development time. Additionally, in Figure 7, there are pit gain vs. time curves for two kicks that are used to evaluate the kick severity. Kick 1 and Kick 2’s pit gains are located in Area ① and Area ②, respectively. Comparing the Kick 1 curve and the Kick 2 curve, the pit gain of Kick 1 is 3.5 m³, while Kick 2 has a pit gain of 2.5 m³ in the same time. The pit gain before the well closes is the direct parameter of the kick tolerance and is used to judge whether the well can be shut down safely and to select killing methods. The kick tolerance is a quantity that represents the amount of residual energy in the process of well killing after a kick occurs. The larger the pit gain’s volume, the smaller the kick tolerance, the greater the well shut-down risk, and the greater the risk of conventional killing. Thus, the kick severity of Area ① is greater than that of Area ②. The results are as follows: in Area ①, the pit gain is very large, with a rapid development trend and the greatest risk of kick; in Area ②, the pit gain is large, with a rapid development trend and a very high kick risk; in Area ③, the pit gain is less, but its development trend is rapid, and the kick risk is greater; in Area ④, the pit gain is very large, there is a certain development trend, and the kick risk is very high; in Area ⑤, the pit gain is large, there is a certain development trend, and the kick risk is relatively large; in Area ⑥, the pit gain is less, there is a slower development trend, and the kick risk is relatively small; in Area ⑦, the pit gain is very large, but the pit gain develops slowly, and there is a large kick risk; in Area ⑧, the pit gain is large, but the development is slow, and there is a certain kick risk; and in Area ⑨, the pit gain is small, the pit gain develops slowly, and the kick risk is small. The order of the kick severity before closing the well is ① > ② > ④ > ⑤ > ⑦ > ③ > ⑧ > ⑥ > ⑨.

3.9.2. Analysis of the Kick Severity after Shut-in

After well shut-in, the casing pressure will gradually increase under the influence of wellbore afterflow. When the bottom hole pressure and formation pressure are balanced, the casing pressure tends to be stable, and this point is the "analysis point" that is used to analyze the severity of the kick risk.

The increase in value of the casing pressure during shut-in depends on the difference between the formation pressure and the bottom hole pressure—the greater the difference, the greater the casing pressure, and the more serious the kick risk. Through the analysis of the casing pressure, the size of the unbalance between the bottom hole pressure and the formation pressure can be evaluated, and then the size of the kick risk can be analyzed.

When the underbalanced pressure difference between the formation pressure and the bottom hole pressure is constant, the formation coefficient is proportional to the increase in the casing pressure, that is, the slope of the casing pressure and time can be used to characterize the formation coefficient. The larger the slope is, the larger the formation coefficient is, and the faster the casing pressure increases, the more serious the kick risk is [37]. Through the analysis of the casing pressure recovery curve, the size of the kick risk can be analyzed. The kick severity analysis chart based on the casing pressure recovery curve during the well shut-in period is shown in Figure 3.

In the shut-in casing pressure recovery curve, the pressure value of the analysis point represents the degree of unbalance between the formation pressure and the bottom hole pressure and several pressure value horizontal lines are used to distinguish the casing pressure level. The slope from the analysis point to the starting point indicates the speed at which the casing pressure recovers, which is affected by the formation coefficient. Therefore, the kick severity analysis chart during the well shut-in period includes several horizontal lines dividing the pressure level and several oblique lines dividing the rate of change in the casing pressure. According to information on some wells in the Sichuan Basin in China, the low-medium pressure limit is 4.5 MPa, the medium-high pressure limit is 9 MPa, the lower slope is tan20°, and the medium-high slope is tan45°, as shown in Figure 8. It should be noted that the specific values of the pressure limit and the curve slope in other regions were adjusted according to the different stratigraphic conditions of the corresponding regions [37].

Figure 8 is divided into nine regions according to the casing pressure and the casing pressure increase rate. Additionally, in Figure 8, there are casing pressure vs. time curves for two kicks that were used to evaluate the kick and well control severity. Kick 3 and Kick 4’s casing pressures are located in Area ① and Area ⑦, respectively. Comparing the Kick 3 curve and Kick 4 curve, it takes 10 min and 45 min to restore a stable casing pressure, respectively. Although the underbalanced pressure difference between the formation pressure and the bottom hole pressure is the same, the velocity of pressure transfer varied greatly. As the velocity of pressure transfer depends largely on the formation coefficient, such as permeability, porosity, and reservoir thickness, the larger the formation coefficient, the faster the pressure transfer velocity. During a well killing operation, by adjusting the opening of the ground valve, the bottom hole pressure should be controlled within a safe range above the formation pressure. However, in a high temperature and high pressure deep well, the complex law of bottom hole pressure changes making it difficult to maintain a safe range, which could cause a secondary kick and increases the difficulty of well killing. Owing to being under the same time and underbalanced pressure conditions, the larger the formation coefficient, the more gas is intruded, and the greater the well control risk. Thus, the kick and well control severity of Area ① is greater than that of Area ⑦. The results are as follows: in Area ①, the kick is very serious and the risk to well control is very large; in Area ②, the kick is serious and the well control risk is large; in Area ③, the kick severity is medium and the well control has a certain difficulty risk; in Area ④, the kick is more serious and the risk to well control is larger; in Area ⑤, the overflow is small and the well control has some risks; in Area ⑥, there is a slight risk of a kick and no well control; in Area ⑦, the kick is small and the well control has some risks; in Area ⑧, there is a slight risk of a kick and no well control; and in Area ⑨, there is no overflow risk and no well control risk. The order of the kick severity during the well closed period is ① > ② > ④ > ⑤ > ⑦ > ③ > ⑧ > ⑥ > ⑨.

4. Results and Discussion

In this section, the proposed method is applied to the kick severity forecast in several wells, which is based on the prediction of pit gain and casing pressure.

4.1. Analysis and Calculation of the Kick Severity before Closing the Well

4.1.1. Case Study One

In Section 3, the pit gain forecasting model ARIMA (4,2,1) of well SC-X1 was established; Figure 9 shows the measured data and predicted values of pit gain vs. time in different periods.

Figure 9a–c show that with an increase in prediction time, the MRE increases gradually (from 8.3% to 8.7% and then to 14%). Specifically, the prediction accuracy decreases with an increase in prediction time, which can also be obtained from Figure 9d,e. By comparing Figure 9a,d, we can see that the more measured data there are, the smaller the MRE is (8.3% and 8%), which can also be obtained from Figure 9b,e. Overall, the MRE based on the measurement data period of 0 to ~8.5 min is small (less than 8.2%).

Figure 10 is the result of forecasting the pit gain after 5 min using the ARIMA (4,2,1) model. The forecasted curve is plotted in the kick severity analysis chart based on the pit gain volume increment. It can be seen from the figure that the data are located in Area ⑥ before the prediction, the pit gain has not reached 2 m³ at 12 min, and the well has not been closed at this time. By predicting the pit gain after 5 min, it can be seen that the pit gain curve is already in Area ④. At this time, the predicted pit gain has reached 3.18 m³ and is also showing a gradually increasing trend. The well control risk is already large—the well needs to be immediately shut down and the appropriate well control measures need to be selected.

4.1.2. Case Study Two

In this case, one well in [3] is used for analysis; the time sequence diagram of the pit gain before well shut-in is shown in Figure 11. Using the same modeling process as described in Section 3, the model ARIMA (2,2,1) is established. Figure 12 shows the measured data and predicted values of the pit gain vs. time in different periods.

From Figure 12, the same rule as Figure 9 can be obtained. Figure 12a–c show that with an increase in the prediction time, the MRE increases gradually (from 8.6% to 13.4% and then to 20.2%). Specifically, the prediction accuracy decreases with the increase in prediction time, which can also be obtained from Figure 12d,e. By comparing Figure 12a,d, it can be seen that the more measured data there are, the smaller the MRE is (8.6% and 3.6%), which can also be obtained from Figure 12b,e. Overall, the MRE based on the measurement data period 0 to ~129 s is small (less than 5.4%).

Figure 13 is the result of the forecasted pit gain in the next 54 s using the ARIMA (2,2,1) model. The forecasted curve is plotted in the kick severity analysis chart based on the pit gain volume increment. It can be seen from the figure that the data are located in Area ⑥ before the prediction, the pit gain is less, there is a slower development trend, and the kick risk is relatively small. By predicting the pit gain after 54 s, it can be seen that the pit gain curve is already in Area ③, and although the pit gain is less, its development trend is rapid and the kick risk is greater. Thus, it needs to be immediately shut down and the appropriate well control measures need to be selected.

4.2. Analysis and Calculation of the Kick Severity after Well Shut-In

Figure 14 is the time sequence diagram of the casing pressure during well shut-in, and it can be seen that there is a marked upward trend in the sequence diagram. Figure 14 shows the autocorrelation coefficient and partial autocorrelation coefficient of the casing pressure time series. It shows that the autocorrelation coefficient decreases to zero slowly, and the autocorrelation coefficient graph shows obvious triangular symmetry. Therefore, the casing pressure sequence is not a stationary sequence and needs differential processing.

Figure 15 shows the autocorrelation coefficient and partial autocorrelation coefficient of the casing pressure time series after the first-order difference; the first-order difference series has a strong short-term correlation.

Table 7. Unit root test results of the casing pressure time series after the first-order difference.

		$\mathit{\tau }$ Statistics	Probability P
Test statistics of ADF		−5.5	0.0019
Critical value of the test	1% Level	−4.3	- -
	5% Level	−3.2	- -
	10% Level	−2.75	- -

shows the unit root test results of the casing pressure time series after the first-order difference; the probability (P) of the ADF test statistic is less than 0.05. Therefore, the series after the first-order difference is stable, and the model can be established.

Based on the autocorrelation and partial autocorrelation diagram of the first-order difference of the casing pressure shown in Figure 16, the preliminary order of the model can be determined. It can be seen from the figure that both the autocorrelation coefficient and partial autocorrelation coefficient fall within two standard deviations after the delay of the first-order, and the model is established by ARIMA (1,1,1).

Table 8. White noise test results and AIC and SBC values of the residual sequences.

Model	AIC	SBC	Delay Order	Q-Statistics	P-Value	Whether to Pass the Test
ARIMA (1,1,1)	−3.26	−3.15	6	2.82	0.59	Yes

shows the residual white noise test of the ARIMA (1,1,1) model. The probability P-value of the Q-statistic is greater than 0.05, so the residual sequence is a white noise sequence, which satisfies the model’s requirements.

Table 9. Model parameters of ARIMA (1,1,1).

Variable	C	AR(1)	MA(1)
Coefficient	0.19	0.76	−0.99
Probability	0.006	0	0

shows the ARIMA (1,1,1) model’s unknown parameter prediction and the $\tau$ statistic’s probability P-value. The P-value is less than 0.05, so the model meets the requirements.

From Figure 17, the same rules as Figure 9 and Figure 12 can be obtained. Figure 17a,b show that with an increase in prediction time, the MRE increases gradually (from 2.1% to 4%). Specifically, the prediction accuracy decreases with the increase of prediction time. Generally speaking, the MRE is small (less than 4%).

Figure 18 shows the change in casing pressure after well shut-in. The ARIMA (1,1,1) model is constructed by the time series method to predict the future casing pressure. It is predicted that the reading point of the casing pressure will be in Area ④, the kick will be more serious, and the risk to well control will be larger, so it needs to be dealt with in time.

4.3. Discussion

The application of this paper can be divided into two stages—the prediction of pit gain before well shut-in and the prediction of casing pressure after well shut-in. After the occurrence of a kick, the prediction model is established based on the real-time data of the previous pit gain, the development trend of the pit gain is predicted, and the next step is taken with reference to the predicted pit gain. If the predicted pit gain is small (for example less than 2 m³), for the commonly used managed pressure drilling technology at present, the well may not need to be shut in. By adjusting the surface equipment of managed pressure drilling, the bottom hole pressure can be increased, the kick can be recycled out, and drilling time can be saved. If the predicted overflow is large (for example more than 2 m³), closing the well before 2 m³ can reduce formation fluid intrusion, which will be helpful to reduce the risk of a killing operation. After well shut-in, it is generally necessary to wait for casing pressure to return to a stable state on the drilling rig, then calculate the formation pressure, design a well killing construction report, and configure heavy mud. The prediction model in this paper can save killing time by predicting the casing pressure so that drilling operators can make preparations for killing the well in advance.

This method is suitable for short-term prediction but is not very good for long-term prediction. Our future work will consider the multimodel coupling method to improve long-term prediction accuracy.

5. Conclusions

Studies and field applications of measured data during drilling are kick warning and have been used to determine whether a kick occurs or not. However, there is no research on forecasting the severity of a kick by using measured data after a kick occurs.

In view of the time-varying characteristics of the characterization parameters after a kick, combined with the high degree of fitting a time series analysis method, the established time series model can be used for the prediction of pit gain and casing pressure.

Since pit gain and casing pressure can indirectly reflect the bottom hole condition, the changes in pit gain and casing pressure with time can be used to judge the severity of the kick. Combined with the prediction model, the severity of the kick can be predicted.

In this paper, the proposed method has been applied to pit gain forecasting in two wells and casing pressure forecasting in one well. The same rules can be seen in the case studies—an increase in prediction time, a gradual increase in MRE, and the more measured data, the smaller the MRE. Future work will consider the multimodel coupling method to improve long-term prediction accuracy.

The comparison between measured data and predicted data has shown that the proposed method has high prediction accuracy and repeatability. It has a certain guiding significance for the prediction of kick development trends, the analysis of kick severity, the selection of well control measures by engineers, and the reduction of the risk of a killing operation.