2.1 Flow equations

To analyse the flow of liquid with gas release we have to assume that the initial liquid contains two chemical components (at least): light (1) and heavy (2). When pressure decreases due to liquid extraction, the light component releases and creates the free gas. The gas will be assumed to consist only of the light component, i.e. the heavy component does not dissolve in gas.

For the sake of simplicity, thermal effects, capillary pressure, diffusion and deformations of the porous medium are neglected. We also accept the local phase equilibrium described by Henry’s law. Gas is ideal, while the liquid is incompressible (the liquid compressibility is neglected with respect to that of gas).

The mass and momentum balance equations of each chemical component are, for $i=l,g$ , in the one-dimensional domain $\mathfrak{D}=\{0<x<L\}$ , where $L$ is the reservoir thickness:

where $\unicode[STIX]{x1D70C}_{l},\unicode[STIX]{x1D70C}_{g}$ are the densities of liquid and gas ( $\text{kg}~\text{m}^{-3}$ ); ${\mathcal{S}}$ is the liquid saturation; ${\mathcal{C}}$ is the mass fraction of the light component in liquid; ${\mathcal{P}}$ is the pressure; $\boldsymbol{V}_{l}$ and $\boldsymbol{V}_{g}$ are Darcy’s velocities of the liquid and gas; $h$ is Henry’s constant ( $\text{Pa}^{-1}$ ); $\unicode[STIX]{x1D62E}$ is the molecular mass of the light component ( $\text{kg}~\text{mol}^{-1}$ ); $\boldsymbol{g}$ is the gravitational acceleration, $K$ is the absolute permeability, $k_{i}$ is the relative permeability of phase $i$ , $\unicode[STIX]{x1D707}_{i}$ is the dynamic viscosity of phase $i$ (assumed to be constant); $R$ is the universal gas constant; $q_{l}$ and $q_{g}$ are the volumetric rates of extraction of liquid and gas reported to $1~\text{m}^{3}$ of the medium [ $1/s$ ]; the vertical axis $x$ is directed upwards. Functions $k_{i}({\mathcal{S}})\in [0,1]$ are monotonic continuous and have the following properties:

where ${\mathcal{S}}_{\ast }$ and ${\mathcal{S}}^{\ast }$ are the percolation thresholds for liquid and gas respectively.

The system of $6$ equations (2.1) determines $6$ variables: ${\mathcal{P}},{\mathcal{S}}$ , $\unicode[STIX]{x1D70C}_{g}$ , ${\mathcal{C}}$ , $\boldsymbol{V}_{g}$ and $\boldsymbol{V}_{l}$ .

To formulate the boundary conditions we should take into account the following. The diffusive fluxes of gas through the bottom coming from the Earth depth and through the reservoir top (gas leakage) may be neglected during the period of reservoir exploitation. Therefore, the reservoir bottom is impermeable for both phases. This is not the case of the reservoir top, as both liquid and gas may be produced just from the reservoir top (and not through the overall reservoir thickness). In this case, the fluid extraction is better expressed through the boundary condition than through a distributed sources $q_{l}$ and $q_{g}$ .

At the initial state the fluid is assumed to be single-phase liquid and immobile, i.e. the initial pressure is hydrostatic.

Consequently, we obtain the following general boundary and initial conditions:

where ${\mathcal{P}}_{\ast }$ is the initial pressure at the reservoir bottom; $V_{l\ast }$ and $V_{g\ast }$ are the velocity of liquid and gas extraction from the reservoir top.

2.2 Introduction of the fractional flow

The phase velocities $\boldsymbol{V}_{l}$ and $\boldsymbol{V}_{g}$ are usually replaced by the total velocity $\boldsymbol{V}$ through the following relationships:

where

Indeed, if we take the sum of $\boldsymbol{V}_{l}$ and $\boldsymbol{V}_{g}$ and use Darcy’s law (2.1c ), then: $\boldsymbol{V}=-\unicode[STIX]{x1D706}\unicode[STIX]{x1D735}{\mathcal{P}}+(\unicode[STIX]{x1D706}_{l}\unicode[STIX]{x1D70C}_{l}+\unicode[STIX]{x1D706}_{g}\unicode[STIX]{x1D70C}_{g})\boldsymbol{g}$ . Then:

Then from Darcy’s law:

which leads to (2.4).

Function ${\mathcal{F}}$ is the fractional flow of liquid in the case without gravity. It has the following properties: ${\mathcal{F}}({\mathcal{S}})\in [0,1]$ , it is a monotonic nonlinear function of saturation and ${\mathcal{F}}(0)\equiv 0$ when $0\leqslant {\mathcal{S}}\leqslant {\mathcal{S}}_{\ast }$ ; ${\mathcal{F}}(1)\equiv 1$ when ${\mathcal{S}}^{\ast }\leqslant {\mathcal{S}}\leqslant 1$ .

The physical meaning of ${\mathcal{F}}({\mathcal{S}})$ is the ratio between the liquid flow rate ( $Q_{l}$ ) and the total flow rate ( $Q_{l}+Q_{g}$ ) at zero gravity force. Indeed, let ${\mathcal{F}}=(Q_{l}/Q_{l}+Q_{g})=(V_{l}/V_{l}+V_{g})$ , then using Darcy’s law (in the case of zero gravity), we obtain (2.5).

The initial and boundary conditions (2.3) become:

where $V_{\ast }=V_{l\ast }+V_{g\ast }$ .

2.3 Characteristic times and scales

The process is characterized by three times:

1. (i) the time of propagation of a pressure perturbation: $\unicode[STIX]{x1D635}_{p}\equiv (\unicode[STIX]{x1D707}_{l}L^{2}\unicode[STIX]{x1D719}/K{\mathcal{P}}_{\ast })$ ;

2. (ii) the time of gas rise due to the buoyancy force: $\unicode[STIX]{x1D635}_{h}\equiv (\unicode[STIX]{x1D707}_{g}L\unicode[STIX]{x1D719}/K\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{0}\,g)$ , $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{0}\equiv \unicode[STIX]{x1D70C}_{l}-\unicode[STIX]{x1D70C}_{g}^{0}$ ;

3. (iii) the time of complete fluid withdrawal: $\unicode[STIX]{x1D635}_{\ast }\equiv (1/(q_{l}+q_{g}))$ in the case of distributed sources, or $\unicode[STIX]{x1D635}_{\ast }\equiv (L/(V_{l\ast }+V_{g\ast }))$ in the case of point sources at the reservoir top.

where ${\mathcal{P}}_{\ast }$ enters in (2.3).

The main three hypotheses of this paper, related to the introduced scales, consist of the following:

$\boldsymbol{{\mathcal{H}}}\mathbf{1}$

the ratio between the characteristic times: $\unicode[STIX]{x1D635}_{p}\ll \unicode[STIX]{x1D635}_{\ast }$ . Then

(2.9) $$\begin{eqnarray}0\leqslant \unicode[STIX]{x1D700}\equiv {\displaystyle \frac{\unicode[STIX]{x1D635}_{p}}{\unicode[STIX]{x1D635}_{\ast }}}\ll 1.\end{eqnarray}$$

$\boldsymbol{{\mathcal{H}}}\mathbf{2}$

the initial pressure distribution along $x$ is almost constant:

(2.10) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x1D70C}_{l}^{0}gL}{{\mathcal{P}}_{\ast }}}\sim \unicode[STIX]{x1D700}.\end{eqnarray}$$

This assumption is valid for thin and deep reservoirs, in which the hydrostatic pressure difference over the reservoir height $\unicode[STIX]{x1D70C}_{l}^{0}gL$ is negligible with respect to the initial pressure ${\mathcal{P}}_{\ast }$ .

$\boldsymbol{{\mathcal{H}}}\mathbf{3}$

small pressure gradients: the square of the pressure gradient, $(\unicode[STIX]{x1D735}{\mathcal{P}})^{2}$ and ( $\boldsymbol{V}_{i}\boldsymbol{\cdot }\unicode[STIX]{x1D735}{\mathcal{P}}$ ) is neglected. This hypothesis is compatible with the preceding one.

The third time $\unicode[STIX]{x1D635}_{h}$ is the true characteristic time of the analysed process. It may take any value. The separation of pressure from the saturation is possible in two cases:

1. (i) slow gas ascension: $\unicode[STIX]{x1D635}_{h}\sim \unicode[STIX]{x1D635}_{\ast }$ ;

2. (ii) fast gas ascension: $\unicode[STIX]{x1D635}_{h}\sim \unicode[STIX]{x1D635}_{p}$ .

2.4 Reduction to the asymptotic model of kinematic waves

Using the hypotheses of the previous section, it is possible to transform system (2.1) into the following system of first-order differential equations for saturation and velocity:

where $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{l}-\unicode[STIX]{x1D70C}_{g}({\mathcal{P}}),r=(\unicode[STIX]{x1D70C}_{l}/\unicode[STIX]{x1D70C}_{g}({\mathcal{P}}))$ , $\mathfrak{T}({\mathcal{P}})=-((1/{\mathcal{P}})\unicode[STIX]{x2202}_{t}{\mathcal{P}})^{-1}$ , $q=q_{l}+q_{g}$ .

The initial and boundary conditions result from (2.8):

In this system the pressure ${\mathcal{P}}$ is a parameter, which is determined as the solution of the independent boundary value problem, which is different for the two different cases of $\unicode[STIX]{x1D635}_{h}$ mentioned in the previous section.

For the case of fast gas ascension: $\unicode[STIX]{x1D635}_{h}\sim \unicode[STIX]{x1D635}_{p}\ll \unicode[STIX]{x1D635}_{\ast }$ :

we obtain for the function $p\equiv {\mathcal{P}}-{\mathcal{P}}_{\ast }+\unicode[STIX]{x1D70C}_{l}gx$ :

where $\unicode[STIX]{x1D714}=(\unicode[STIX]{x1D719}h(r_{\ast }-1)/(1-{\mathcal{C}}_{\ast }))$ , $\unicode[STIX]{x1D706}_{l0}=\unicode[STIX]{x1D706}_{l}({\mathcal{S}})|_{{\mathcal{S}}=1}=K/\unicode[STIX]{x1D707}_{l}$ , $r_{\ast }=(\unicode[STIX]{x1D70C}_{l}/\unicode[STIX]{x1D70C}_{g}({\mathcal{P}}_{\ast }))$ , ${\mathcal{C}}_{\ast }=h{\mathcal{P}}_{\ast }$ .

For the case of fast depletion ( $\unicode[STIX]{x1D635}_{h}\sim \unicode[STIX]{x1D635}_{\ast }$ ), the pressure is calculated as the solution to the following Cauchy problem for an ordinary differential equation:

where $\overline{{\mathcal{S}}}\equiv (1/L)\int _{x=0}^{x=L}{\mathcal{S}}(x,t)\,\text{d}x$ . The value $\overline{{\mathcal{S}}}$ should be considered as known and may be accepted in an approximate way. In particular, if we assume that $\overline{{\mathcal{S}}}=\text{const}.$ , $\unicode[STIX]{x1D70C}_{l}\gg \unicode[STIX]{x1D70C}_{g}$ and ${\mathcal{C}}\ll 1$ , then equation (2.14) has an analytical solution:

The derivation of all these relationships is given in the next section.

The physical meaning of various terms of equation (2.11) is as follows: term I is the saturation transport forced by the non-zero total velocity; term II is the transport by the buoyancy force, while III, IV and V are the saturation ‘production’ due to phase transition, gas expansion and non-proportionate liquid/gas extraction. The extraction is called proportionate if the production rates of the phases are proportional to their fractional flow. It is seen that for a proportionate extraction, term V disappears.

Note that the transport terms in the left-hand side are not affected by the phase transitions, which have an impact only on the source terms. Therefore, in the case without phase transitions, the model transforms into the classical Buckley–Leverett with the same convective part.

Equation (2.11) is a non-homogeneous extension of the Buckley–Leverett model known in the theory of oil recovery. The extension consists of two elements:

1. (i) first of all, the source terms appear in the right-hand side;

2. (ii) secondly, the total velocity is non-constant (it is constant in the classical Buckley–Leverett model), and is determined, in turn, from a differential equation.

2.5 Derivation of the model. Step I: reduction to a non-separated wave form

The derivation of the model (2.11)–(2.14) is done by two steps. At the first step we obtain the first-order equations (2.11) for ${\mathcal{S}}$ and $\boldsymbol{V}$ and the following problem for pressure

1. (1) It results from Henry’s law and $\boldsymbol{{\mathcal{H}}}\mathbf{3}$ :

   (2.17) $$\begin{eqnarray}\displaystyle d_{l}{\mathcal{C}}\equiv \unicode[STIX]{x2202}_{t}{\mathcal{C}}+\boldsymbol{U}_{l}\unicode[STIX]{x1D735}{\mathcal{C}}=h\unicode[STIX]{x2202}_{t}{\mathcal{P}}+h\boldsymbol{U}_{l}\unicode[STIX]{x1D735}{\mathcal{P}}=h\unicode[STIX]{x2202}_{t}{\mathcal{P}}+\cdots \,. & & \displaystyle\end{eqnarray}$$

2. (2) The sum of equations (2.1a ) and (2.1b ) yields

   (∗) $$\begin{eqnarray}\unicode[STIX]{x1D719}\,\unicode[STIX]{x2202}_{t}(\unicode[STIX]{x1D70C}_{l}{\mathcal{S}}+\unicode[STIX]{x1D70C}_{g}(1-{\mathcal{S}}))+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{l}\boldsymbol{V}_{l}+\unicode[STIX]{x1D70C}_{g}\boldsymbol{V}_{g})=-\unicode[STIX]{x1D719}\unicode[STIX]{x1D70C}_{g}q_{g}-\unicode[STIX]{x1D719}\unicode[STIX]{x1D70C}_{l}q_{l}.\end{eqnarray}$$

   Differentiating by parts (2.1b ) and (∗) we obtain

   (2.18) $$\begin{eqnarray}\displaystyle (1-{\mathcal{C}})(\unicode[STIX]{x1D719}\unicode[STIX]{x1D70C}_{l}\,\unicode[STIX]{x2202}_{t}{\mathcal{S}}+\unicode[STIX]{x1D70C}_{l}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{V}_{l}+\unicode[STIX]{x1D719}\unicode[STIX]{x1D70C}_{l}q_{l})-\unicode[STIX]{x1D719}\unicode[STIX]{x1D70C}_{l}{\mathcal{S}}\,d_{l}{\mathcal{C}}=0, & & \displaystyle\end{eqnarray}$$
   (2.19) $$\begin{eqnarray}\displaystyle & & \displaystyle (-\unicode[STIX]{x1D719}\unicode[STIX]{x1D70C}_{g}\unicode[STIX]{x2202}_{t}{\mathcal{S}}+\unicode[STIX]{x1D719}(1-{\mathcal{S}})\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D70C}_{g}+\unicode[STIX]{x1D70C}_{g}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{V}_{g}+\unicode[STIX]{x1D719}\unicode[STIX]{x1D70C}_{g}q_{g})\nonumber\\ \displaystyle & & \displaystyle \qquad \qquad +\,(\unicode[STIX]{x1D719}\unicode[STIX]{x1D70C}_{l}\unicode[STIX]{x2202}_{t}{\mathcal{S}}+\unicode[STIX]{x1D70C}_{l}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{V}_{l}+\unicode[STIX]{x1D719}\unicode[STIX]{x1D70C}_{l}q_{l})+\cdots =0.\end{eqnarray}$$
   This leads to
   (2.20a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}\unicode[STIX]{x2202}_{t}{\mathcal{S}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{V}_{l}+\unicode[STIX]{x1D719}q_{l}=\frac{\unicode[STIX]{x1D719}{\mathcal{S}}}{(1-{\mathcal{C}})}\,d_{l}{\mathcal{C}}=\left\{(2.17)\right\}=\frac{\unicode[STIX]{x1D719}{\mathcal{S}}h}{(1-{\mathcal{C}})}\unicode[STIX]{x2202}_{t}{\mathcal{P}}\equiv -{\mathcal{M}}\quad & \displaystyle\end{eqnarray}$$

   (2.20b ) $$\begin{eqnarray}\displaystyle & \displaystyle -\unicode[STIX]{x1D719}\unicode[STIX]{x2202}_{t}{\mathcal{S}}+\frac{\unicode[STIX]{x1D719}(1-{\mathcal{S}})}{{\mathcal{P}}}\unicode[STIX]{x2202}_{t}{\mathcal{P}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{V}_{g}+\unicode[STIX]{x1D719}q_{g}=-\frac{\unicode[STIX]{x1D719}r{\mathcal{S}}}{(1-{\mathcal{C}})}\,d_{l}{\mathcal{C}}=r{\mathcal{M}}.\quad & \displaystyle\end{eqnarray}$$

3. (3) Summing (2.20a ) and (2.20b ) we obtain the second equation in (2.11).

4. (4) Replacing $\boldsymbol{V}_{l}$ in (2.20a ) by $\boldsymbol{V}$ through (2.4), we obtain

   (∗∗) $$\begin{eqnarray}\unicode[STIX]{x1D719}\unicode[STIX]{x2202}_{t}{\mathcal{S}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }({\mathcal{F}}\,\boldsymbol{V}-\unicode[STIX]{x1D706}_{g}{\mathcal{F}}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\,\boldsymbol{g})=\frac{\unicode[STIX]{x1D719}{\mathcal{S}}h}{(1-{\mathcal{C}})}\unicode[STIX]{x2202}_{t}{\mathcal{P}}-\unicode[STIX]{x1D719}q_{l}.\end{eqnarray}$$

5. (5) Multiplying the second equation (2.11) by ${\mathcal{F}}$ and subtracting it from (∗∗), we obtain the first equation (2.11).

6. (6) Using the second equation in (2.11) and the relationship (2.6) between the velocity and pressure, we can remove the velocity and obtain the differential equation (2.16a ) for the pressure. Conditions (2.16b ), (2.16c ) and (2.16d ) result from (2.8).

2.6 Derivation of the model. Step II: separation of pressure

At the second step we separate completely the pressure by applying hypotheses $\boldsymbol{{\mathcal{H}}}\mathbf{1}$ and $\boldsymbol{{\mathcal{H}}}\mathbf{2}$ and the asymptotic technique.

2.7 Hyperbolicity of the model (2.11)

Equations (2.11) may be presented as

where ${\mathcal{M}}\equiv \unicode[STIX]{x1D719}\,{\mathcal{C}}\,{\mathcal{S}}/(\mathfrak{T}(1-{\mathcal{C}})),{\mathcal{H}}\equiv -{\mathcal{M}}(1-{\mathcal{F}}+r{\mathcal{F}})-(\unicode[STIX]{x1D719}{\mathcal{F}}(1-{\mathcal{S}})/\mathfrak{T})+\unicode[STIX]{x1D719}q_{g}{\mathcal{F}}-\unicode[STIX]{x1D719}q_{l}(1-{\mathcal{F}}),{\mathcal{I}}\equiv (r-1){\mathcal{M}}+(\unicode[STIX]{x1D719}(1-{\mathcal{S}})/\mathfrak{T})-\unicode[STIX]{x1D719}q$ .

This system of two equations of first order with respect to ${\mathcal{S}}$ and $V$ can be represented in the form of Cauchy–Kowalevski, explicit for the derivatives $\unicode[STIX]{x2202}_{x}$ :

where

where the symbol ‘prime’ stands for the derivative with respect to ${\mathcal{S}}$ .

If the denominator in $\unicode[STIX]{x1D6FC}$ is non-zero for all ${\mathcal{S}}$ , then the eigenvalues of matrix $\unicode[STIX]{x1D63C}$ exist and are real: $\unicode[STIX]{x1D709}_{1}=0$ and $\unicode[STIX]{x1D709}_{2}=\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D719}$ , and there exists an orthonormal basis of the left eigenvectors of matrix $\unicode[STIX]{x1D63C}$ : $l^{1}=(1,0)$ and $l^{2}=(0,1)$ . Then system (2.28) is hyperbolic (Rhee, Aris & Amundson Reference Rhee, Aris and Amundson1986).

3.1 Problem set-up

Consider a closed reservoir, such that $V_{g\ast }=V_{l\ast }=0$ . The production of the fluid occurs throughout the overall thickness and is determined by the rates $q_{l}$ and $q_{g}$ . The main feature of such a system consists of the fact that the total velocity $\boldsymbol{V}$ is zero, which will be shown below.

Let the liquid and gas extraction rates be proportionate, i.e. the production rate of each phase is proportional to its fractional flow: $(q_{l}/q_{g})=((1-{\mathcal{F}})/{\mathcal{F}})$ . Then the term $V$ in (2.11) disappears. In practice, this means that gas and liquid are extracted by the same producing well, therefore their rates depend on each other. For pressure, we will use the asymptotic solution (2.13).

The problem of segregation resulting from (2.11) and (2.13) is reduced to a single equation for the saturation:

where $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ , $r$ and ${\mathcal{C}}$ depend on pressure, which is a parameter of the problem defined through (2.13), which has, in this case, an explicit analytical solution independent of $x$ : $p=p(t)=-((1-{\mathcal{C}}_{\ast })/h(r_{\ast }-1))\int _{0}^{t}q\,\text{d}t^{\prime }$ , so that

Parameter $\mathfrak{T}$ is:   $\mathfrak{T}=-((1/{\mathcal{P}})\unicode[STIX]{x2202}_{t}{\mathcal{P}})^{-1}={\mathcal{C}}(r_{\ast }-1)/(q(1-{\mathcal{C}}_{\ast }))$ .

Function ${\mathcal{G}}({\mathcal{S}})$ describes the gas transport due to the buoyancy force, while the right-hand side ( $-{\mathcal{H}}$ ) describes the decrease in liquid saturation in the overall domain due to liquid degassing and gas expansion. Function ${\mathcal{F}}({\mathcal{S}})$ is the fractional flow of liquid without gravity.

Functions ${\mathcal{F}}({\mathcal{S}})$ and ${\mathcal{G}}({\mathcal{S}})$ are shown in figure 1. Note that their derivatives are zero at points ${\mathcal{S}}=0$ and ${\mathcal{S}}=1$ , according to the properties of the relative permeabilities. The derivative $\text{d}{\mathcal{G}}/\text{d}{\mathcal{S}}$ determines the speed of the wave propagation and its direction. One distinguishes two families of waves:

1. (i) propagating upwards ( ${\mathcal{G}}^{\prime }<0$ ), which corresponds to high liquid saturation. Physically, this is caused by gas rise due to the buoyancy force.

2. (ii) propagating downwards ( ${\mathcal{G}}^{\prime }>0$ ), which corresponds to high gas saturation. Such a wave is the result of gas expansion.

Figure 1. Diagrams ${\mathcal{F}}({\mathcal{S}})$ and ${\mathcal{G}}({\mathcal{S}})$ for the problem of gas–liquid segregation.

3.2 Steady-state solution

Equation (3.1) has steady-state solutions which satisfy the ordinary differential equation:

and which are the result of the equilibrium between the wave propagation (term $\unicode[STIX]{x2202}_{x}{\mathcal{G}}$ ) and the gas generation (term ${\mathcal{H}}$ ).

Its solutions are shown in figure 2.

Figure 2. Steady-state solutions of the problem of gas–liquid segregation.

These solutions are non-single valued and have a singularity at point ${\mathcal{S}}={\mathcal{S}}_{m}$ of the maximum of the function ${\mathcal{G}}({\mathcal{S}})$ , in which the derivative $\text{d}{\mathcal{S}}/\text{d}x$ becomes infinite. Only the upper part of these curves above the dashed line has a physical meaning. It represents the equilibrium between gas rise, which leads to the increase of the local liquid saturation, and the gas generation, which leads to the decrease of ${\mathcal{S}}$ .

The red curve, which corresponds to the boundary condition ${\mathcal{S}}|_{x=0}=1$ , represents the main interest. We will show that this curve makes part of the non-stationary solutions.

3.3 General structure of the solution of (3.1)

The boundary condition in (3.1) means that the saturation at point $x=L$ is either ${\mathcal{S}}=1$ or ${\mathcal{S}}=0$ . The case ${\mathcal{S}}=1$ has no meaning (no gas at the reservoir top), consequently it is expected that ${\mathcal{S}}|_{x=L}=0$ , $\forall t>0$ . In other words, the free gas appears immediately at the reservoir top, by creating the initial shock of saturation between ${\mathcal{S}}=1$ on the left of it and ${\mathcal{S}}=0$ on the right. Therefore, one expects to have discontinuous solutions.

Consequently, the solution of the boundary initial problem for quasi-linear hyperbolic equation (3.1) is expected to consist of four elements:

1. (i) A rarefaction/compression wave (‘rc-wave’), which is a non-trivial continuous function presented by a fragment of the curve ${\mathcal{G}}({\mathcal{S}})$ and propagating at the velocity ${\mathcal{G}}^{\prime }({\mathcal{S}})$ along the characteristics $x=x(t)$ defined by the system of equations:

   (3.5a,b ) $$\begin{eqnarray}\frac{\text{d}{\mathcal{S}}}{\text{d}t}=-{\mathcal{H}}({\mathcal{S}}),\quad \frac{\text{d}x}{\text{d}t}=-{\mathcal{G}}^{\prime }({\mathcal{S}}),\end{eqnarray}$$
   which are different for various ${\mathcal{S}}$ ;

2. (ii) a stationary rc-wave determined by the problem (3.4), which is the equilibrium between gas rise and gas generation. It does not move but lengthens in time.

3. (iii) A plateau, which is constant in $x$ but varies in time. It satisfies the problem:

   (3.6) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}{\mathcal{S}}=-{\mathcal{H}}({\mathcal{S}}),\quad {\mathcal{S}}|_{t=0}=1.\end{eqnarray}$$
   For a fixed time moment it is represented by a point in the diagram ${\mathcal{G}}(S)$ .

4. (iv) A shock, which is a strong discontinuity between two rc-waves, between an rc-wave and a plateau or between two plateaux. It is determined by the following conditions:

   1. (a) the mass balance equation (the Hugoniot–Rankine condition) relating the shock velocity $U_{s}$ and the saturation behind it ${\mathcal{S}}^{-}$ and ahead of it ${\mathcal{S}}^{+}$ :

      (3.7) $$\begin{eqnarray}U_{s}=\frac{{\mathcal{G}}^{+}-{\mathcal{G}}^{-}}{{\mathcal{S}}^{+}-{\mathcal{S}}^{-}}.\end{eqnarray}$$
      In the diagram ${\mathcal{G}}({\mathcal{S}})$ the shock velocity is determined by the tangent of the straight line relating the points behind the shock ( ${\mathcal{G}}^{-},{\mathcal{S}}^{-}$ ) and ahead of it ( ${\mathcal{G}}^{+},{\mathcal{S}}^{+}$ ).

   2. (b) The entropy conditions of Lax and Oleinik, which require the shock velocity to be between the velocities of the rc-waves that are connected to it and to not cross the curve ${\mathcal{G}}({\mathcal{S}})$ . If the shock is connected to only one mobile rarefaction wave, the entropy condition implies that the shock velocity must be equivalent to the rc-wave velocity at the contact point. Graphically this means that the straight line (3.7) is tangent to the curve ${\mathcal{G}}({\mathcal{S}})$ at point of contact with a mobile rc-wave.

3.4 Solution at the first stage

First of all, we consider the case when the coefficients of (3.1) do not depend on pressure, i.e. ${\mathcal{P}}={\mathcal{P}}_{\ast }$ , then $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ , $r$ and ${\mathcal{C}}$ are constant.

The solution of the problem is obtained analytically by applying the method of characteristics. Knowing all the structural elements of the solution mentioned in § 3.3, we can calculate their velocities and saturations at the points of junction by using analytical or semi-analytical relationships obtained previously: (3.4), (3.6), (3.7). The ordinary differential equations have been solved by using Mathlab.

The solution of the problem is presented in figure 3.

Figure 3. Variation of the liquid saturation at the first stage (a) and the corresponding diagram ${\mathcal{G}}({\mathcal{S}})$ (b). The dashed curve is the stationary wave.

The average value of the derivative ${\mathcal{G}}^{\prime }({\mathcal{S}})$ of the convective term in (3.1) tends to zero (according to figure 1), then the ‘average’ equation (3.1) has the form (3.6), whose solution is a plateau variable in time. This corresponds to the saturation decrease in the reservoir due to progressive liquid degassing. Consequently, the main element of the solution becomes a set of plateaux $aa$ , $bb$ , $cc$ , $dd$ , $ee$ , $gg$ , At the same time, the buoyancy force, governed by the convective term, maintains a high liquid saturation at the reservoir bottom ( ${\mathcal{S}}\rightarrow 1$ , at $x\rightarrow 0$ ) and low liquid saturation at the reservoir top ( ${\mathcal{S}}\rightarrow 0$ , at $x\rightarrow L$ ).

At the reservoir top, gas is accumulated, i.e. ${\mathcal{S}}=0$ , as mentioned in § 3.3. The connection between the low saturation and the plateau occurs through a shock, whose progressive positions in space are $bo$ , $cp$ , $dq$ , $er$ , $gu$ . This shock represents the wave reflected from the reservoir top and propagating downwards. It appears, since the rc-waves connecting a high and a low saturation produce non-single-valued solutions. Indeed, the rc-wave $c{-}d{-}e{-}g{-}u{-}r{-}q{-}p{-}o{-}0$ , which passes from point $c$ of the plateau $cc$ to point $({\mathcal{G}},{\mathcal{S}})=(0,0)$ , has negative transport velocity for high saturations (points $c,d,e$ ), but positive velocity for low saturations (points $u,r,q,p,o$ ), which implies the following formal solution in terms of ${\mathcal{S}}(x)$ (figure 4). Such a solution is not single valued, cannot exist and must be replaced by a shock.

Figure 4. Example of a non-single-valued solution described by the rc-wave $c{-}d{-}e{-}g{-}u{-}\cdots {-}o$ .

Shock $bo{-}cp{-}dq{-}er{-}gu$ marks the boundary of the gas cap. It has positive velocity, i.e. is transported downwards, which means the growth of the gas cap. The shock velocity increases in time. At points of connection with the rc-wave ( $o,p,q,r,u$ ) the straight line of this shock is tangent to the curve ${\mathcal{G}}({\mathcal{S}})$ , according to the entropy conditions.

The total system of equations that defines these shocks is: equation (3.7), in which ${\mathcal{S}}^{+}$ is the solution of (3.6), while ${\mathcal{S}}^{-}$ is the solution of ${\mathcal{G}}^{\prime }({\mathcal{S}}^{-})=U_{s}$ , which is the entropy condition. This system of equations is closed and determines the shock in a unique way.

At the reservoir bottom, gas saturation is low, the process is controlled by the right-hand part of the diagram ${\mathcal{G}}({\mathcal{S}})$ , which corresponds to the gas rise. The connection between the high saturation and the plateau occurs via a stationary rc-wave described by (3.4) at the boundary condition ${\mathcal{S}}|_{x=0}=1$ . It passes through points $a,b,c,d,e,g$ along the right-hand branch of the curve ${\mathcal{G}}({\mathcal{S}})$ .

The first stage continues up to the moment when the solution reaches the singular point $g$ , at which ${\mathcal{S}}={\mathcal{S}}_{m}$ , and which corresponds to the final point of the steady-state rc-wave.

3.5 Solution at the second and the third stages

At the reservoir bottom, after the stationary rc-wave $a{-}b{-}c{-}d{-}e{-}g$ reaches the singular point $g$ (at which $S={\mathcal{S}}_{m}$ ), the gas expansion becomes more important than the gas rise, which leads to the appearance of a descending wave. This wave is however non-single valued, similar to preceding paragraph, and transforms into a new shock propagating downwards (figure 5). Its consecutive positions are shown as $hk$ , $im$ and $jn$ . This shock connects the plateau and the stationary wave, which enables us to determine its parameters in a unique way by (3.7), in which ${\mathcal{S}}^{+}$ is the solution of (3.4) with the boundary condition ${\mathcal{S}}|_{x=0}=1$ , while ${\mathcal{S}}^{-}$ is determined by (3.6).

Figure 5. Variation of the liquid saturation at the second stage (a) and the corresponding diagram ${\mathcal{G}}({\mathcal{S}})$ (b).

At the reservoir top, the gas cap boundary becomes less sharp but more continuous, which means the transformation of the shock $gu$ into an rc-wave propagating downwards. The transport velocity of this wave is determined by the tangents lines to the left-hand part of the function ${\mathcal{G}}({\mathcal{S}})$ , which are, obviously, larger than the velocities of the shock $hk$ – $im$ – $jn$ . Consequently, the upper rarefaction waves should collide with the down shock.

After the collision, the plateau disappears. The solution has a single shock, which continue to propagate downward. This is the third stage. Curve $1{-}j{-}n{-}a$ marks the beginning of this stage.

3.6 Comparison with numerical solution

The comparison of the analytical solutions presented above has been made with numerical solutions obtained by two different methods. One of them is based on finite element discretization of the domain, which is realized in the simulator Comsol Multiphysics, while the second one is based on the finite volumes implemented in the open source code DuMux. Both methods use the fully implicit scheme of discretization in time (the fully implicit means that all coefficients and the right-hand part are implicit), which is unconditionally stable.

Comsol Multiphysics allows the use on non-structured grids, but the discrete schemes are non-necessary conservatives, as for any finite element method.

DuMux (Flemisch et al. Reference Flemisch, Darcis, Erbertseder, Faigle, Lauser, Mosthaf, Muthing, Nuske, Tatomir and Wolff2011) is based on the algebraic code DUNE (Distributed and Unified Numerics Environment) (Bastian et al. Reference Bastian, Blatt, Dedner, Engwer, Klofkorn, Kornhuber, Ohlberger and Sander2008) and the module MpNc (M-phases, N-components) that simulates multicomponent and multiphase flow and transport in porous media. The spatial discretization is done by the Box-method, which is a finite volume method (vertex-centred) based on non-structured non-orthogonal grids. The principle ideas of this method were already presented in Patankar (Reference Patankar1980) in terms of the ‘control volume techniques’. The more modern version may be found in Bank, Rose, Fichtner (Reference Bank, Rose and Fichtner1983) and Kramer, Nicholas & Hitchon (Reference Kramer, Nicholas and Hitchon1997). The finite volume method traditionally uses orthogonal grids, since in this case the expressions for fluxes through the faces of finite elements are greatly simplified. The use of non-orthogonal grids (as in the method of finite elements) is a significant advantage of the Box-method.

For the one-dimensional problem studied herein, both methods lead to close numerical schemes.

For hyperbolic systems having discontinuous solutions, both methods apply the technique of low diffusion.

Two cases are presented in figures 6 and 7.

Figure 6. Comparison of the analytical solution (the solid lines) with numerical simulations by DuMux (the black points) and Comsol Multiphysics (the grey points), for $\unicode[STIX]{x1D700}=0.2$ .

Figure 7. Comparison of the analytical solution (the solid lines) with numerical simulations by DuMux (the black points) and Comsol Multiphysics (the grey points), for $\unicode[STIX]{x1D700}=0.4$ .

The characteristic times have been selected as follows: $\unicode[STIX]{x1D635}_{h}=0.5$ , $\unicode[STIX]{x1D635}_{\ast }=2$ . The dimensionless parameters that determine the process are: $\unicode[STIX]{x1D700}=0.2$ , $\unicode[STIX]{x1D70C}^{0}gL/{\mathcal{P}}_{\ast }=0.2$ , $\unicode[STIX]{x1D70C}_{l}/\unicode[STIX]{x1D70C}_{g}^{0}=5$ , ${\mathcal{C}}_{l}^{0}=0.2$ , $\unicode[STIX]{x1D6FE}/{\mathcal{P}}_{\ast }=0.2$ , $h/{\mathcal{P}}_{\ast }=0.2$ , $\unicode[STIX]{x1D719}=0.2$ . The perturbation time is  $\unicode[STIX]{x1D635}_{p}=0.4$ when $\unicode[STIX]{x1D700}=0.2$ , or $\unicode[STIX]{x1D635}_{p}=0.8$ when $\unicode[STIX]{x1D700}=0.4$ , and so on.

Despite the satisfactory results, we should note that the total velocity $V$ is not exactly constant, as predicted by the asymptotic model. Figure 8 shows the behaviour of the dimensionless total velocity for four moments of time.

Figure 8. Total velocity, exact numerical results for four moments of time.

Figure 9. Evolution of the gas cap structure in time: stages I, II and III.

Note that the maximal/minimal values of these fluctuations are very small, between $-0.001$ and $0.003$ .

3.7 Structure of the gas cap

According to the results obtained, the formation of a gas cap passes through three stages shown in figure 9.

As seen, the zone at the bottom, highly saturated with liquid, is separated from the gas cap by a gas cap underlayer, in which the gas saturation is almost invariable in space. The growth of the gas cap leads to the contraction of this underlayer, while keeping almost constant the thickness of the liquid zone.

4.1 Problem formulation

Let the fluid be extracted just from the reservoir top. Then in the problem (2.13): $V|_{x=L}=V\ast$ , and $q_{l}=q_{g}=0$ . Let, as previously, the liquid and gas extraction rates be proportionate, i.e. $(V_{l\ast }/V_{g\ast })=((1-{\mathcal{F}})/{\mathcal{F}})$ .

Then the problem of segregation resulting from (2.11) and (2.13) becomes:

where ${\mathcal{G}}({\mathcal{S}})=(1/\unicode[STIX]{x1D719})(\unicode[STIX]{x1D706}_{g}{\mathcal{F}}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\,g)$ , and parameter $\mathfrak{T}(x,t)$ is $\mathfrak{T}^{-1}=-(1/{\mathcal{P}})\unicode[STIX]{x2202}_{t}{\mathcal{P}}\sim -(1/{\mathcal{P}}_{\ast })\unicode[STIX]{x2202}_{t}{\mathcal{P}}$ .

In this system the pressure is determined from the independent problem (2.13) formulated for function $p={\mathcal{P}}-{\mathcal{P}}_{\ast }+\unicode[STIX]{x1D70C}_{l}gx$ :

which has an approximate analytical solution obtained by the integral method of Kármán–Pohlhausen:

where $\unicode[STIX]{x1D635}_{p}$ is determined as the solution to the equation:

Calculating the derivative $\unicode[STIX]{x2202}_{t}p$ from (4.3) we obtain the explicit relationship for $\mathfrak{T}$ :

In the case $V_{\ast }=\text{const}.$ , we obtain $\unicode[STIX]{x1D635}_{p}=(L^{2}\unicode[STIX]{x1D714}/6\unicode[STIX]{x1D706}_{l0})$ , while the functions $l(t)$ and $\mathfrak{T}(x,t)$ become

where ${\mathcal{H}}(t)$ is the Heaviside function.

4.2 Derivation of the solution to the problem for pressure

The solution of (4.2) may be obtained by the method of integral relationships of Kármán–Pohlhausen.

Let us assume that the pressure is perturbed in the finite zone $l(t)<x\leqslant L$ (the perturbation comes from point $x=L$ ), as shown in figure 10, where $l(t)$ is the unknown mobile coordinate.

Figure 10. Approximation of the function $p(x,t)$ , defined through (2.14), in the method of Kármán–Pohlhausen.

Within this zone, the pressure is a polynomial function which satisfies the boundary value conditions and the differential equation (2.14) in the average:

Then we obtain for functions $a(t),b(t),c(t)$ and $l(t)$ :

or

The last equation can be integrated, which gives: $(l-L)^{2}=-(6\unicode[STIX]{x1D706}_{l0}/\unicode[STIX]{x1D714}V_{\ast })\int _{0}^{t}V_{\ast }(t^{\prime })\,\text{d}t^{\prime }$ . This leads to (4.3).

4.3 Results of solution

The solution of the problem is presented in figure 11.

Figure 11. Evolution of the liquid saturation in the case of gas fluid production from the reservoir top, for seven moments of time.

As seen, the upper shock appears, as in the previous example, but remains immobile during all the period. The gas cap is determined by the lower shock, which progressively moves down.