# Two stochastic versions of the Arps curve decline

Christian Paroissin

► To cite this version:

Christian Paroissin. Two stochastic versions of the Arps curve decline. Communications in Nonlinear Science and Numerical Simulation, 2026, 152 (Part D), pp.109400. 10.1016/j.cnsns.2025.109400. hal-04764329

HAL Id: hal-04764329

https://univ-pau.hal.science/hal-04764329v1

Submitted on 3 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Two stochastic versions of the Arps curve decline

Ch. Paroissin

Last version: November 3, 2024

This work is licensed under a Creative Commons “Attribution 4.0 International” license.

<figure><img src="image_2.png" /></figure>

**Abstract** Based on the Arps equation, we propose two stochastic models for curve decline useful in oil engineering context. Theoretical properties and simulations of these models are provided. The first passage time distribution of these stochastic models to a constant level is then studied. In conclusion, we discuss about statistical inference of the parameters from the observations of the oil production cumulative rate.

**Keywords** curve decline · Arps equation · Stochastic differential equation · Hitting time distribution · Stochastic comparison

**Mathematics Subject Classification (2000)** 60H10 · 60K30 · 62M05

## 1 Introduction

Arps [5] (see also Chapter 9 in [18])) introduced several decades ago an empirical decline equation during pseudo-steady state period. This curve decline is of the following form:

$$ \forall t \ge 0, \quad q_t = \frac{q_0}{(1 + b d_0 t)^{1/b}}, \qquad (1) $$

where $q_t$ is the oil production rate at production time $t$ and $q_0$ is the initial oil production rate and where $b$ and $d_i$ are two constants representing respectively a shape parameter and the initial flow. The case $b = 0$ represents an exponential decline:

$$ \forall t \ge 0, \quad q_t = q_0 e^{-d_0 t}, $$

and $b = 1$ represents a harmonic decline in oil production. In other cases ($0 < b < 1$), Equation (1) is defined as the hyperbolic model. Although this decline curve proposed by Arps is based only on an empirical study, this equation is still widely used in oil industry studies (see [10] or [14], for instance).

Here we propose stochastic models for curve decline analysis (section 2). These models are based on stochastic differential equations such that the expectation is the Arps equation. Uncertainty in decline curve analysis were already considered. For instance, Chang and Lin [8] have proposed a regression linear analysis of the decline curve after applying a transformation to fit data. We then study the important question of how long will the oil production rate be larger than a certain level (section

<footnote>Christian Paroissin</footnote>

<footnote>Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, LMAP, Pau, France. E-mail: christian.paroissin@univ-pau.fr</footnote>
3). Such problem is difficult and no exact computations can be carried out (mainly numerical ones). In conclusion, we discuss about statistical estimation of the parameters. Most of the time the oil production cumulative rate is observed.

## 2 Stochastic Arps model

We introduce here a set of stochastic differential equations (SDEs) modelling curve decline for oil production for instance. These SDEs are of the form:

$$dQ_t = - \frac{d_0 Q_t}{1 + b d_0 t} dt + \alpha(Q_t) dB_t, \quad (2)$$

with initial condition $Q_0 = q_0$ and where $(B_t)$ is a Brownian motion and $\alpha(\cdot)$ a real function. These SDEs involve the two parameters $d_0$ and $b$ that appear in the Arps equation.

When $b$ tends to zero, both the drift and the volatility do not depend no more on time (only on space); in such case the corresponding diffusion process is said to be homogeneous. In the general setting, if the function $\alpha(\cdot)$ satisfies the following assumptions:

(A1) (Lipschitz condition) there exists a constant $\kappa_1 > 0$ such that:

$$\forall t \in [0, T], \forall (x, y) \in \mathbb{R}^2, \quad |\alpha(x) - \alpha(y)| \le \kappa_1 |x - y|,$$

(A2) (linear growth bound) there exists a constant $\kappa_2 > 0$ such that:

$$\forall t \in [0, T], \forall x \in \mathbb{R}, \quad |\alpha(x)| \le \kappa_2^2 (1 + |x|^2),$$

then the corresponding SDE admits an unique strong solution on $[0, T]$, applying theorem 4.5.3 in [17] for instance. Two special cases, for which the above conditions are clearly satisfied, will be studied below. The first one is when the volatility is constant and the second is a case where the volatility is linearly increasing. It follows that these two special cases leads to scalar linear SDEs (see [17], page 110).

### 2.1 Case with constant volatility

In this subsection we will assume that the volatility is constant: $\alpha(x) = \sigma > 0$. Thus Equation (2) turns to be:

$$dQ_t = - \frac{d_0 Q_t}{1 + b d_0 t} dt + \sigma dB_t, \quad (3)$$

This choice of the function $\alpha(\cdot)$ can be motivated as follows: this SDE can be interpreted as a stochastic perturbation of the ordinary differential equation satisfied by $q_t$.

When the parameter $b$ tends to 0, the SDE is the so-called Ornstein-Uhlenbeck process [22]. The solution is then given by:

$$\forall t \ge 0, \quad Q_t = q_0 e^{-d_0 t} + \sigma e^{-d_0 t} \int_0^t e^{d_0 u} dB_u. \quad (4)$$

As for this special case, the solution can be computed:

**Proposition 1.** The solution of Equation (3) is:

$$\forall t \ge 0, \quad Q_t = q_0(1 + b d_0 t)^{-1/b} + \sigma(1 + b d_0 t)^{-1/b} \int_0^t (1 + b d_0 u)^{1/b} dB_u. \quad (5)$$
Proof. For all $t \ge 0$, set $X_t = (1 + bd_0t)^{1/b}Q_t$ with initial value $X_0 = Q_0 = q_0$. It follows that:

$$ 
\begin{aligned}
dX_t &= d_0(1 + bd_0t)^{1/b-1}Q_tdt + (1 + bd_0t)^{1/b}dQ_t \\
&= \sigma(1 + bd_0t)^{1/b}dB_t.
\end{aligned}
$$

Thus it follows that:

$$ X_t = X_0 + \sigma \int_0^t (1 + bd_0u)^{1/b} dB_u. $$

Back to the initial process, one obtains the following solution of the SDE:

$$ \forall t \ge 0, \quad Q_t = q_0(1 + bd_0t)^{-1/b} + \sigma(1 + bd_0t)^{-1/b} \int_0^t (1 + bd_0u)^{1/b} dB_u. $$

□

The integral in the solution is a Wiener one, it follows that ($Q_t$) is a Gaussian process with mean:

$$ \forall t \ge 0, \quad \mathbb{E}[Q_t] = q_0(1 + bd_0t)^{-1/b} = q_t $$

and variance:

$$ \forall t \ge 0, \quad \mathrm{Var}[Q_t] = \sigma^2(1 + bd_0t)^{-2/b} \int_0^t (1 + bd_0u)^{2/b} \, du = \frac{\sigma^2(1 + bd_0t)}{d_0(2+b)} \left(1 - (1 + bd_0t)^{-2/b-1}\right). $$

In addition, the covariance between $Q_s$ and $Q_t$, for any $(s, t) \in \mathbb{R}_+^2$, is given by:

$$ \mathrm{Cov}[Q_s, Q_t] = \sigma^2(1 + bd_0s)^{-1/b}(1 + bd_0t)^{-1/b} \int_0^{\min(s,t)} (1 + bd_0u)^{2/b} \, du. $$

It follows that the expectation of this SDE is precisely the Arps equation. When considering Equation (5) and with $b = 0$, one recovers the solution given in Equation (4) for the Ornstein-Uhlenbeck process.

One can express ($Q_t$) as a Brownian motion with a time scale change. Indeed applying the Dambis-Dubins-Schwartz theorem [23], there exists a Brownian motion ($W_t$) such that:

$$ Z_t = \int_0^t (1 + bd_0u)^{1/b} dB_u \stackrel{(d)}{=} W_{\tau(t)}, $$

where $\stackrel{(d)}{=} \text{ stands for equality in distribution and where:}$

$$ \tau(t) = \mathrm{Var}[Z_t] = \frac{1}{d_0(b+2)} \left[ (1 + bd_0t)^{\frac{2}{b}+1} - 1 \right]. $$

This expression is the same as the one in [4] when the parameter $b$ tends to zero. In any case, as a consequence, the following representation holds:

$$ Q_t \stackrel{(d)}{=} (1 + bd_0t)^{-1/b} (q_0 + \sigma W_{\tau(t)}). \quad (6) $$

On Figure 1 we present simulations of this SDE for various values of the parameter $b$ and with $d_0 = 3 \times 10^{-4}$, $q_0 = 380$ and $\sigma^2 = 1$. To do so, we have considered the Euler-Maruyama method (see [17], for instance). We have used the same simulations of the normal distribution in order to compare the different paths. For fixed value of $t$, $q_t$ is decreasing with parameter $b$. Hence, using the same simulations of normally distributed random variables, it follows that, for any $t$, sample paths are increasing with $b$.
<figure><img src="image_3.png" /><figcaption>Fig. 1: Simulations of the first stochastic models</figcaption></figure>

## 2.2 Case with linear volatility

Now we will consider a case where the volatility is increasing. More precisely we assume that $\alpha(x) = \sigma x$ with $\sigma > 0$. Thus Equation (2) turns to be:

$$dQ_t = -\frac{d_0 Q_t}{1 + b d_0 t} dt + \sigma Q_t dB_t, \qquad (7)$$

This choice of the function $\alpha(\cdot)$ can be motivated as follows: this SDE can be interpreted as a stochastic perturbation of the ordinary differential equation satisfied by $\ln(q_t)$.

When the parameter $b$ tends to 0, the SDE is the geometric Brownian motion [22] and the solution is given by:

$$\forall t \ge 0, \quad Q_t = q_0 e^{-d_0 t} \exp \left( -\frac{1}{2} \sigma^2 t + \sigma B_t \right). \qquad (8)$$

As for this special case, the solution can be computed:

**Proposition 2.** *The solution of Equation (7) is:*

$$\forall t \ge 0, \quad Q_t = q_0 (1 + b d_0 t)^{-1/b} \exp \left( -\frac{1}{2} \sigma^2 t + \sigma B_t \right). \qquad (9)$$
Proof. Let us divide Equation (7) by $Q_t$. It gives:

$$ \frac{dQ_t}{Q_t} = -\frac{d_0}{1+bd_0t}dt + \sigma dB_t. $$

It follows that:

$$ \int_{0}^{t} \frac{dQ_u}{Q_u} = -\frac{1}{b} \ln(1 + bd_0t) + \sigma B_t. \quad (10) $$

Thus one has to evaluate the left hand-side of this equation. Let us apply the Itô formula [22] to $\ln Q_t$:

$$ d \ln Q_t = \frac{dQ_t}{Q_t} - \frac{1}{2} \sigma^2 dt $$

Integrating this equation, one gets:

$$ \int_{0}^{t} \frac{dQ_{u}}{Q_{u}} = \ln\left(\frac{Q_{t}}{q_{0}}\right) + \frac{1}{2}\sigma^{2}t. \quad (11) $$

Joining Equations (10) and (11), it gives:

$$ \ln\left(\frac{Q_t}{q_0}\right) = -\frac{1}{b}\ln(1 + bd_0t) - \frac{1}{2}\sigma^2t + \sigma B_t. $$

Finally one gets:

$$ Q_t = q_0(1 + bd_0t)^{-1/b} \exp\left(-\frac{1}{2}\sigma^2t + \sigma B_t\right). $$

Clearly, for any $t$, $Q_t$ has the log-normal distribution with parameters $q_t - \frac{1}{2}\sigma^2t$ and $\sigma^2t$. Thus one gets that:

$$ \forall t \ge 0, \quad \mathbb{E}[Q_t] = q_0(1 + bd_0t)^{-1/b} $$

and variance:

$$ \forall t \ge 0, \quad \mathrm{Var}[Q_t] = q_0^2(1 + bd_0t)^{-2/b} (e^{\sigma^2t} - 1). $$

Thus, as for the previous model, the expectation of this SDE is the Arps equation.

On Figure 2 we present the simulation of this SDE for various values of the parameter $b$ and with $d_0 = 3 \times 10^{-4}$, $q_0 = 380$ and $\sigma^2 = 10^{-2}$. As previously, we have considered the Euler-Maruyama method [17] and we have used the same simulations of the normal distribution in order to compare the different paths.

### 3 Hitting time problem

An important question that could raise in petroleum context is the following one: when the production rate at a given site will be lower than a certain threshold? Thus one can be interested in the hitting time $T_{Q,x}$ of level $x < q_0$ for the stochastic process $Q = (Q_t)$:

$$ T_{Q,x} = \inf\{t; Q_t \le x\} = \inf\{t; Q_t = x\}, $$

since $Q$ has continuous sample paths. Such random variable has been well studied in the case of exponential decline ($b = 0$) or harmonic decline ($b = 1$) since it corresponds to some classical processes: the Ornstein-Uhlenbeck process (see [4] or [20] for instance) in the first case and the geometric Brownian motion in the second one. In fact, for any value of $b$, the problem can be reduced to the first-passage time of a Brownian motion through a one-sided moving barrier. Thus the first subsection is devoted to a (non-exhaustive) review of some results about it. In addition, we provide a lemma that will be used later. Then in the two last subsections we discuss about the hitting time distributions for our problem. In particular we will compare first passage time distribution in function of the parameter $b$. To this purpose we will also denote $T_{Q,x}(b)$ instead of $T_{Q,x}$.
<figure><img src="image_4.png" /><figcaption>Fig. 2: Simulations of the second stochastic models</figcaption></figure>

## 3.1 Brownian motion and hitting times

The first-passage time of a Brownian motion through a boundary is not an easy problem and the distribution of this random variable is not known explicitly except for a constant or a linear barrier. In these cases the hitting time distribution is called the inverse Gaussian distribution.

Let $(B_t)$ be a standard Brownian motion and $c(\cdot)$ a real function. Consider the first-passage time of $(B_t)$ through the boundary $c(\cdot)$:

$$T_{B,c(\cdot)} = \inf\{t; B_t = c(t)\}.$$

If the $c(\cdot)$ is constant or linear, then $T_{B,c(\cdot)}$ has the inverse Gaussian distribution [13]. More precisely, if $c(t) = at + b$, then the probability distribution function of $T_{B,c(\cdot)}$ is given by:

$$f(t) = \frac{|b|}{\sqrt{2\pi t^3}} \exp\left(-\frac{(at+b)^2}{2t}\right).$$

It corresponds to the inverse Gaussian distribution with parameters $(|b|/a, b^2)$. The expectation of $T_{B,c(\cdot)}$ is equal to $|b|/a$ and the variance to $|b|^3$. The Laplace transform of such probability distribution has been obtained by Abdel-Hameed and Nakhi [2] (see also [1], page 19):

$$\mathcal{L}_{T_{B,c(\cdot)}}(s) = \mathbb{E}[e^{-sT_{B,c(\cdot)}}] = \exp\left(-|b|\left(\sqrt{a^2 + 2s} - a\right)\right).$$
Based on these two special cases (constant and linear boundaries), several authors study the case of constant piecewise and linear piecewise boundaries (see [3] and references therein). A general expression (not necessary leading to an explicit formula, excepted for some cases) was given by Salminen [25] in the case where $c(\cdot)$ is twice continuously differentiable on $[0, +\infty)$.

We will compare hitting time distribution according to parameters values (especially the parameter $b$). We first recall the definition of the stochastic comparison in the usual sense [21]:

**Definition 1.** The random variable $X$ is said to be stochastically smaller than the random variable $Y$ if and only if:

$$ \forall t \in \mathbb{R}, \quad \mathbb{P}[X > t] \le \mathbb{P}[Y > t]. $$

It will be denoted as follows: $X \preceq_{st} Y$.

We recall that $X \preceq_{st} Y$ implies that $\mathbb{E}[X] \le \mathbb{E}[Y]$. Now we give the lemma that will be used later (proofs are omitted since there is no real difficulty):

**Lemma 1**

1. If $X$ and $Y$ are two positive random variables defined on the same probability space such that $X \preceq_{st} Y$ and if $f$ is a one-to-one map defined on $\mathbb{R}^+$, then $f(X) \preceq_{st} f(Y)$.

2. If $f$ and $g$ are two functions such that for any $t \in \mathbb{R}$ $f(t) \le g(t)$, then for any random variable $X$ $f(X) \preceq_{st} g(X)$.

3. Let $f$ and $g$ be two continuous functions defined on $\mathbb{R}^+$ such that:

$$ \forall t \ge 0, \quad f(t) \ge g(t), $$

$$ \text{with } g(0) \le f(0) \le 0. \quad \text{Then } T_{B,f(\cdot)} \preceq_{st} T_{B,g(\cdot)}. $$

## 3.2 Case with constant volatility

In the case of exponential decline, the hitting time of a fixed value has been studied recently. There is no closed form expression for the probability distribution function [4]. However expressions were obtained that yield numerical computations. For instance, the Laplace transform of $T_{Q,x}$ can be expressed in term of Hermite functions $H_\nu$ or in term of parabolic cylinder functions $D_\nu$ [7,27]:

$$ \mathbb{E} [e^{-uT_{Q,x}}] = \frac{H_{-u/d_0} \left( \frac{q_0}{\sigma} \sqrt{d_0} \right)}{H_{-u/d_0} \left( \frac{x}{\sigma} \sqrt{d_0} \right)} = \frac{e^{d_0 q_0^2 / 2\sigma^2} D_{-u/d_0} \left( \frac{q_0}{\sigma} \sqrt{2d_0} \right)}{e^{d_0 x^2 / 2\sigma^2} D_{-u/d_0} \left( \frac{x}{\sigma} \sqrt{2d_0} \right)}. $$

Then one can get numerically the p.d.f. of $T_{Q,x}$ applying for an appropriate algorithm (see [6] for instance). Otherwise a way to compute the hitting time distribution is to consider the representation of $(Q_t)$ by a time rescaled Brownian motion. The following lemma states the relationship between these two hitting time problems:

**Proposition 3.** For any $t \ge 0$, let $c_1(t) = \frac{q_0}{\sigma} \left( \frac{x}{q_{\tau^{-1}(t)}} - 1 \right)$. Let $W$ be a Brownian motion. Then $T_{Q,x}$ and $\tau^{-1}(T_{W,c_1(\cdot)})$ have the same distribution.

*Proof.* It is an application of the representation of $(Q_t)$ given by Equation (6). Indeed, for any $t \ge 0$, one gets:

$$ \begin{aligned} \mathbb{P}[T_{Q,x} > t] &= \mathbb{P}[\forall u \in [0, t], Q_u \ge x] \\ &= \mathbb{P}[\forall u \in [0, t], q_u \left( q_0 + \frac{\sigma}{q_0} W_{\tau(u)} \right) > x] \\ &= \mathbb{P}[\forall r \in [0, \tau(t)], q_{\tau^{-1}(r)} \left( q_0 + \frac{\sigma}{q_0} W_r \right) > x] \\ &= \mathbb{P}[\forall r \in [0, \tau(t)], W_r > \frac{q_0}{x} \left( \frac{x}{q_{\tau^{-1}(r)}} - 1 \right)] \\ &= \mathbb{P}[T_{W,c_1(\cdot)} > \tau(t)] \\ &= \mathbb{P}[\tau^{-1}(T_{W,c_1(\cdot)}) > t]. \end{aligned} $$
using that $\tau$ is one-to-one map with $\tau(0) = 0$. $\square$

<figure><img src="image_5.png" /></figure>

From the expression of $\tau$ given in the previous, one can easily compute the inverse function:

$$ \forall t \ge 0, \quad \tau^{-1}(t) = \frac{1}{bd_0} \left[ (1 + t(b + 2)d_0)^{\frac{2}{b+2}} - 1 \right]. $$

Hence the boundary involved in the previous lemma is:

$$ c_1(t) = \frac{q_0}{\sigma} \left( \frac{x}{q_{\tau^{-1}(t)}} - 1 \right) = \frac{1}{\sigma} \left[ x \left( 1 + d_0(b+2)t \right)^{\frac{1}{b+2}} - q_0 \right]. $$

Notice that $c_1(0) = (x - q_0)/\sigma < 0$. As $b$ tends to zero, $c_1$ turns to be a square-root boundary:

$$ c_1(t) = \frac{1}{\sigma} \left[ x \sqrt{1 + 2d_0t} - q_0 \right]. $$

We now compare stochastically the hitting time distributions according to the value of the parameter $b$:

**Proposition 4.** If $b_1 \le b_2$, then $T_{Q,x}(b_1) \le_{st} T_{Q,x}(b_2)$.

*Proof.* Let $b_1 \le b_2$. Applying the third property of Lemma 1, it follows that

$$ T_{W,c_1(\cdot)}(b_1) \le_{st} T_{W,c_1(\cdot)}(b_2). $$

Now, using the first property in this Lemma,

$$ T_{Q,x}(b_1) \stackrel{(d)}{=} \tau_{b_1}^{-1}(T_{W,c_1(\cdot)}(b_1)) \le_{st} \tau_{b_1}^{-1}(T_{W,c_1(\cdot)}(b_2)), $$

and the second property still in Lemma 1 gives

$$ \tau_{b_1}^{-1}(T_{W,c_1(\cdot)}(b_2)) \le_{st} \tau_{b_2}^{-1}(T_{W,c_1(\cdot)}(b_2)) \stackrel{(d)}{=} T_{Q,x}(b_2). $$

Joining these stochastic inequalities lead to the result stated in the proposition. $\square$

<figure><img src="image_6.png" /></figure>

This is consistent with the disadvantages pointed out in the literature when using the Arps equation: the exponential decline ($b = 0$) leads to an underestimation of reserves and production rates while the harmonic decline ($b = 1$) leads to an overestimation of reservoir performance [16,19].

Using the above proposition with $b_1 = 0$ towards with the expression for the mean first-passage time obtained by Sato ([26], see also [24]), we have the following corollary (using the same notations as in [26]):

**Corollary 1.** The mean hitting time for ($Q_t$) to reach the level $x < q_0$ is bounded as follows:

$$ \mathbb{E}[T_{Q,x}] \ge \phi^{(1)}(0, x) - \phi^{(1)}(0, q_0), $$

where:

$$ \phi^{(1)}(0, z) = \frac{1}{d_0} \left( \frac{z\sqrt{d_0\pi}}{\sigma} F\left(\frac{1}{2}, \frac{3}{2}, \frac{d_0z^2}{\sigma^2}\right) + \sum_{m=0}^{\infty} \frac{2^m}{(m+1)(m+2)!!} \left(\frac{z\sqrt{d_0}}{\sigma}\right)^{2m+2} \right). $$

in which $F$ denotes the Kummer function (or confluent hypergeometric function).

Some approximations of the mean first-passage time of a fixed level by the Ornstein-Uhlenbeck process have been studied by Thomas [29].
### 3.3 Case with linear volatility

The exponential case, i.e. with $b = 0$, is easy to derive. In fact, in such case, $T_{Q,x}$ is the hitting time of a constant by a Brownian motion with a linear drift. More precisely, after some easy computations, one gets:

$$T_{Q,x} = \inf \left\{ t; -(d_0 + \frac{1}{2}\sigma^2)t + \sigma B_t \ge \ln \left( \frac{x}{q_0} \right) \right\}.$$

It follows that $T_{Q,x}$ has the inverse Gaussian distribution with parameter:

$$m = \frac{2}{2d_0 + \sigma^2} \ln \left( \frac{q_0}{x} \right) \quad \text{and} \quad \lambda = \left( \frac{1}{\sigma} \ln \left( \frac{q_0}{x} \right) \right)^2. \qquad (12)$$

Let us now consider the general case: $T_{Q,x}$ is the hitting time by a Brownian motion of barrier of the form:

$$\begin{aligned} c_2(t) &= \frac{1}{2}\sigma t - \frac{1}{\sigma}\ln\left(\frac{q_t}{x}\right) \\ &= \frac{1}{2}\sigma t + \frac{1}{\sigma}\ln\left(\frac{x}{q_0}\right) + \frac{1}{b\sigma}\ln(1 + bd_0t). \end{aligned}$$

Using our notations we have $T_{Q,x} \stackrel{(d)}{=} T_{B,c_2(\cdot)}$. The boundary $c_2(\cdot)$ can be bounded as follows by a linear boundary:

$$\forall t \ge 0, \quad c_2(t) \le \tilde{c}_2(t) = \frac{1}{\sigma} \ln \left( \frac{x}{q_0} \right) + \left( \frac{1}{2}\sigma + \frac{d_0}{\sigma} \right) t,$$

the equality holding if $b = 0$. Notice that $c_2(0) = \tilde{c}_2(0) = \log(x/q_0)/\sigma < 0$. More generally, using the third property of Lemma 1, the following proposition holds:

**Proposition 5.** If $b_1 \le b_2$, then $T_{Q,x}(b_1) \le_{st} T_{Q,x}(b_2)$.

As we recall previously, this proposition implies that $\mathbb{E}[T_{Q,x}(b_1)] \le \mathbb{E}[T_{Q,x}(b_2)]$. In particular, $T_{Q,x}(0)$ is inverse Gaussian distributed with parameters given in Equation (12). So we have:

$$\mathbb{E}[T_{Q,x}] \ge \frac{2 \ln(q_0/x)}{2d_0 + \sigma^2}.$$

This is also coherent with the disadvantages pointed out in the literature as we discussed above for the first stochastic model. This lower bound tends to the hitting time of level $x$ in the deterministic exponential curve decline as the parameter $\sigma^2$ tends to zero.

Since the function $c_2(\cdot)$ is continuously differentiable, $T_{Q,x}$ has a continuous probability density function [12]. The exact probability density function $f_{T_{Q,x}}(t)$ of $T_{Q,x}$ cannot be computed explicitly but can be approximated using the result by Durbin [11]. Indeed $f_{T_{Q,x}}(t)$ is approximatively equal to $p(t)f(t)$ with:

$$p(t) = \frac{1}{\sigma t} \ln \left( \frac{q_0}{x} \right) - \frac{1}{b\sigma t} \ln(1 + bd_0t) - \frac{d_0}{\sigma(1 + bd_0t)}$$

and

$$f(t) = \frac{1}{\sqrt{2\pi t}} \exp \left( -\frac{c_2(t)^2}{2t} \right).$$

For various values of $b$ and with $d_0 = 3 \times 10^{-4}$, $q_0 = 380$, $\sigma^2 = 10^{-2}$ and $x = 100$, we have computed numerically the approximation of $f_{T_{Q,x}}(t)$ which are plotted on the Figure 3.
<figure><img src="image_7.png" /><figcaption>Fig. 3: Approximation of the pdf of the hitting time for the second model</figcaption></figure>

## 4 Conclusion

We introduce here two stochastic versions of the Arps equation based on scalar linear SDEs. Solutions have been calculated explicitly. First passage time (FPT) distribution of a fixed level was also considered. Stochastic comparison properties show that the disadvantages pointed out in the literature when using the Arps equation hold also for these models. Further study on these stochastic models could be the statistical inference. In fact, in oil production context, available data are obtained rather from the oil production cumulative rate. It means that one does not observe ($Q_t$) at several times, but the stochastic process ($Q_t^c$) defined as follows:

$$ \forall t \ge 0, \quad Q_t^c = \int_0^t Q_u \, du. $$

On Figure 4 we present simulations of ($Q_t^c$) from the two stochastic models we studied here. Simulations on the left concerns the stochastic process ($Q_t^c$) when ($Q_t$) has a constant volatility while those on the right concerns the stochastic process ($Q_t^c$) when ($Q_t$) has a (positive) linear volatility. From the previous simulations of the two SDEs, it is not so surprisingly that sample paths on the left plot are more smooth than the ones of the right. If ($Q_t$) is a Gaussian process, then ($Q_t^c$) is still a Gaussian process and is a called integrated Gaussian process (this is clearly the case for the first stochastic Arps model we introduced). Statistical inference for such kind of models have been studied previously in the literature, see [15] for instance. Beside, the approach developed here could be adapted to other
<figure><img src="image_8.png" /><figcaption>Fig. 4: Simulations of the cumulative stochastic models</figcaption></figure>

models for flow behaviour of shale gas, like the Matthews-Lefkovits model, the Generalized Weng's prediction model and the Stretched Exponential Decline Production model, see [9] or [28].

References

1. M.S. Abdel-Hameed. Degradation processes: an overview. *Advances in degradation modelling. Applications to reliability, survival analysis and finance*, M.S. Nikulin N. Limnios N. Balakrishan W. Kahle C. Uber-Carol (Eds). Chapter 2: 17–25. Birkhäuser, 2010.

2. M.S. Abdel-Hameed Y. Nakhi. Optimal control of a finite dam using $P_{\lambda,\tau}^M$ policies and penalty cost: total discounted and long run average cost cases. *J. Appl. Probab.*, 28: 888–898, 1990.

3. M. Abundo. Some conditional crossing results of Brownian motion over a piecewise-linear boundary. *Statist. Probab. Letters*, 58(2): 131–145, 2002.

4. L. Alili P. Patie J.L. Pederson. Representations of the first hitting time density of an Ornstein-Ulhenbeck process. *Stochastic Models*, 21: 967–980, 2005.

5. J.J. Arps. Analysis of decline curves. *Trans. AIME*, 160: 228–247, 1945.

6. M.P. Atkinson and S.R. Lang. A comparison of some inverse Laplace transform techniques for use in circuit design. *The Computer Journal*, 15(2): 138–139, 1972.

7. L. Breiman. First exit times from a square root boundary. *In Fifth Berkeley Symposium, 2: 9–16*, 1967.

8. C.P. Chang Z.S. Lin. Stochastic analysis of production decline data for production prediction and reserves estimation. *Journal of Petroleum Science and Engineering*, 23: 149–160, 1999.

9. S. Coutry M. Tantawy S. Fadel. Assessing the accuracy of empirical decline curve techniques for forecasting production in unconventional reservoirs: a case study of Haynesville, Marcellus, and Marcellus Upper Shale. *Journal of Engineering and Applied Science*, 70:69, 2023.

10. L. Dean R. Mireault. Reservoir engineering for geologist. Part 4 - Production decline analysis. *Reservoir: 20-23*, January 2008.

11. J. Durbin. The first-passage density of a continuous Gaussian process to a general boundary. *J. Appl. Probab.*, 22: 99–122, 1985.

12. B. Ferebee. The tangent approximation to one-sided Brownian exit densities. *Z. Wahrscheinlichkeitsh*, 61: 293–302, 1982.

13. J.L. Folks R.S. Chhikara. The inverse Gaussian distribution and its statistical application - A review. *J. Roy. Stat. Soc. (B)*, 40: 263–275, 1978.

14. R. Gaskari S.D. Mohaghegh J. Jalali. An integrated technique for production data analysis with application to mature fields. *SPE Production & Operations*, 22(4): 403–416, 2007.

15. A. Gloter. Discrete sampling of an integrated diffusion process and parameter estimation of the diffusion coefficient. *ESAIM P&S*, 4: 205–227, 2000.

16. H.H. Khanamiri. A non-iterative method for decline curve analysis. *Journal of Petroleum Science and Engineering*, 73(1-2), 59–66, 2010.

17. P.E. Kloeden E. Platen. *Numerical solution of stochastic differential equations*. Springer, 1992.

18. J. Lee R.A. Wattenbargar. *Gas Reservoir Engineering* Society of Petroleum Engineers, 1996.

19. K. Li R.N. Horne. A decline curve analysis model based on fluid flow mechanisms. *SPE* (83470). Presented at the SPE/AAPG Western Regional Meeting, Long Beach, California, 2003.
20. V. Linetsky. Computing hitting time densities for CIR and OU diffusions: applications to mean-reverting models. *J. Comput. Finance*, 7: 1–22, 2004.

21. A. Müller D. Stoyan. *Comparison methods for stochastic models and risks*. Wiley, 2003.

22. B. Oksendal. *Stochastic differential equations*. Springer, 2000.

23. D. Revuz M. Yor. *Continuous martingales and Brownian motion*. Springer, 1999.

24. L.M. Ricciardi S. Sato. First-passage-time density and moments of the Ornstein-Uhlenbeck process. *J. Appl. Probab.*, 25: 43–57, 1988.

25. P. Salminen. On the first hitting time and the last exit time for a Brownian motion to/from a moving boundary. *Adv. Appl. Probab.*, 20: 411–426, 1988.

26. S. Sato. On the moments of the firing interval of the diffusion approximated model neuron. *Mathematical Biosciences*, 39(1-2): 53–70, 1978.

27. A.J.F. Siegert. On the first passage time probability problem. *Physical Rev.*, 81: 617–623, 1951.

28. L. Tan L. Zuo B. Wang; Methods of decline curve analysis for shale gas reservoirs. *Energies*, 11: 552, 20: 2018

29. M.L. Thomas. Some mean first-passage time approximations for the Ornstein-Uhlenbeck process. *J. Appl. Probab.*, 12: 600–604, 1975.