## Abstract

This paper is concerned with the formulation and numerical solution of equations for modelling multicomponent, two-phase, thermal fluid flow in porous media. The fluid model consists of individual chemical component (species) conservation equations, Darcy's law for volumetric flow rates and an energy equation in terms of enthalpy. The model is closed with an equation of state and phase equilibrium conditions that determine the distribution of the chemical components into phases. It is shown that, in the absence of diffusive forces, the flow equations can be split into a system of hyperbolic conservation laws for the species and enthalpy and a parabolic equation for pressure. This decomposition forms the basis of a sequential formulation where the pressure equation is solved implicitly and then the component and enthalpy conservation laws are solved explicitly. A numerical method based on this sequential formulation is presented and used to demonstrate some typical flow behaviour that occurs during fluid injection into a reservoir.

## 1. Introduction

This paper is concerned with the mathematical formulation and numerical solution of systems of partial differential equations for multicomponent, two-phase flow in porous media that include thermal effects. The latter are important in the numerical simulation of enhanced oil recovery processes such as steam flooding and *in situ* combustion. Relevant methodologies have been reported in recent papers by Christensen *et al*. (2004), Nilsson *et al*. (2005), Pasarai *et al*. (2005), Huang *et al*. (2007) and Liu *et al*. (2007). Thermal process simulation also plays an important role in the modelling of geothermal reservoirs (the reader is referred to the recent survey by O'Sullivan *et al*. (2001)). Modelling the behaviour and evaluating clean-up strategies for contaminants can also require the treatment of thermal processes (see Class *et al*. (2002) and references therein for a survey of this type of application).

The development of solution methodologies for the simulation of complex subsurface thermal processes often leads to a compromise between accuracy and stability. Many of these thermal processes are characterized by sharp fronts propagating through the medium, which places stringent requirements on spatial and temporal resolution. However, the simulations also rely on an accurate coupling of a number of different processes, which places stringent requirements on process coupling and stability. The trend in recent years has been to place an emphasis on stability rather than spatial accuracy. Consequently, most approaches to thermal simulations are based on a fully implicit formulation with lower order implicit spatial differencing. Most of the approaches in the literature cited above are based on this type of strategy. It is also worth mentioning the work of Christensen *et al*. (2004) and Nilsson *et al*. (2005), who use adaptive mesh techniques to improve spatial resolution.

The goal of the present work is to develop an algorithm for modelling multiphase, multicomponent non-isothermal reacting flow that achieves a better balance between stability and accuracy. As part of this development, we examine the mathematical structure of these types of flows, which provides the basis for applying contemporary high-resolution upwind methodology in this context. Examples of this type of approach for the black-oil model are discussed in Bell *et al*. (1989). Similar types of discretization procedures are considered for compositional models in Mallison *et al*. (2003).

The approach taken here was used previously to analyse the structure of black-oil (Trangenstein & Bell 1989*a*) and compositional (Trangenstein & Bell 1989*b*) models for reservoir simulation. It represents a generalization of the classical derivation of the Buckley–Leverett equation (e.g. Peaceman 1977). As in the development of compositional models by Trangenstein & Bell (1989*b*), the following analysis relies heavily on the mathematical structure of the multiphase multicomponent phase behaviour. The treatment of phase behaviour for non-isothermal systems using an optimization framework is discussed in Brantferger (1991), as part of the development of a fully implicit thermal–compositional solver, and in Michelsen (1999).

The rest of the paper is structured as follows. In §2, we review the basic equations for multiphase, non-isothermal flow in porous media and discuss the structure of the phase equilibrium problem. In §3, the sequential splitting of the flow equations is introduced. This sequential form separates the system into a parabolic pressure equation and a system of conservation laws for the chemical constituents of the fluid and the enthalpy of the rock–fluid system. In §4, we show that, in the absence of capillary pressure and diffusive transport, this system of conservation laws is hyperbolic. We develop a numerical algorithm based on this formulation in §5. Finally, we present numerical results based on the sequential formulation that illustrate the behaviour of the interacting waves encountered in the system. These examples serve to validate the sequential formulation and illustrate the types of wave behaviour associated with this type of system.

## 2. Mathematical formulation

In this section, we present the basic flow equations describing non-isothermal, multiphase, multicomponent flows in heterogeneous porous media. The porosity of the medium is denoted by *ϕ*, and the phase volume of each phase per pore volume is denoted by *u*_{α}. Greek subscripts refer to mobile phases. (The medium itself, which can be viewed as a solid phase, is given a distinguished treatment since it is immobile.) The multicomponent mixture is composed of *N* components (or lumped pseudo-components) and we define *n*_{α} as the vector of moles of each component in phase *α* divided by the pore volume. Thus, is the total number of moles per pore volume of these components in the combined fluid system. (Mineral content of the solid is, of course, important for geochemistry but is not included in these definitions.)

### (a) Flow equations

The flow is governed by the equations of mass and energy conservation and by Darcy's law, which gives the volumetric flow rate, *v*_{α}, of each phase in terms of the phase pressure, *p*_{α},(2.1)where *K* is the permeability of the medium; *k*_{r,α} is the relative permeability, which expresses the modification of the flow rate from multiphase effects; *η*_{α} and *ρ*_{α} are the phase viscosity and density, respectively; and ** g** is the gravitational force. Here,

*λ*

_{α}≡

*Kk*

_{r,α}/

*η*

_{α}is the phase mobility. The pressure in each phase is related to a reference pressure,

*p*(typically taken to be the most wetting phase), by a capillary pressure,

*p*

_{c,α}=

*p*

_{α}−

*p*, which is a function of saturation.

Conservation of mass for each component is given by(2.2)Here, ** D** are diffusive terms that include multiphase molecular diffusion and dispersion and

*R*_{n}are reaction terms. Both the diffusion and reaction terms can be quite complex, depending on the particular problem (see Chen

*et al*. (2006) for a detailed discussion of these terms).

For non-isothermal systems, it is necessary to include the energy conservation in the system of equations. The overall energy balance must include energy in the solid phase. If we assume that the porous medium and the fluids are in thermal equilibrium, the energy balance is of the form(2.3)where is the total enthalpy of the system; *h*_{α} are the partial molar phase enthalpies; *H*_{r} is the enthalpy of the medium; and *ρ*_{r} is the density of the medium. Here, ** q** represents diffusive energy transport processes such as thermal conduction and

*R*

_{H}represents energy release from reactions and external heating. If one does not assume the porous medium and the fluid to be in thermal equilibrium,

*R*

_{H}will contain relaxation terms that equilibrate the two temperatures. In the formulation (2.3), we have omitted a term from the right-hand side of the formOmitting this term implicitly assumes that the change in phase pressures is slow so that these terms can be ignored. We note that formulations based on internal energy make a similar assumption by dropping a term of the formfrom the energy equation. For a more detailed derivation and discussion, the reader should refer to Burger

*et al*. (1985).

### (b) Phase behaviour

The component conservation and energy equations express the change in the total mass of each component due to advection, diffusion and chemical reactions. Since components are being transported in phases, it is necessary to know the composition of the phases before we can solve the flow equations. This decomposition is referred to as the phase behaviour of the system. The phase behaviour is determined by saying that the equilibrium state of the mixture occurs at the point of maximum entropy. The entropy of the phases is derived from the chemical potential *μ*_{α}. These chemical potentials are typically specified in terms of an equation of state to model the dependence of pressure on temperature, composition and specific volume of the phase. For a more detailed discussion, see Brantferger *et al*. (1991) and Michelsen & Mollerup (2004).

The chemical potentials *μ*_{α} are functions of *p*_{α}, *T* and phase composition *n*_{α}. Note that we have preserved the role of capillary pressure in determining the thermodynamic behaviour of the system. One simplification is to define the thermodynamics in terms of the reference pressure and retain capillary pressure effects only in the definition of phase velocities (see Brantferger *et al*. 1991). Furthermore, the major thermodynamic variables describing each phase can all be expressed in terms of the phase's chemical potential. In particular, the partial molar entropies are given by(2.4)where *x*_{α}=*n*_{α}/*e*^{T}*n*_{α} are the mole fractions; ** e**=(1, 1, …, 1)

^{T}; and the partial molar enthalpies are given by(2.5)

Here, the phase equilibrium problem is to determine the composition of the phases *n*_{α}, given the total moles ** n**, pressure P and the total enthalpy

*H*

_{t}. The equilibrium distribution of the components is given by minimizing the negative entropy of the system. In a two-phase liquid and vapour case, this becomessubject toandalong with inequality constraints (

*n*_{α}≥0) guaranteeing non-negativity of the compositions and thermal stability of the fluid . As noted above, if we do not consider the solid and the fluids to be in thermal equilibrium, then we hold the total fluid enthalpy, , constant rather than the total enthalpy.

Treatment of this minimization problem, i.e. the so-called isenthalpic flash calculation, has been presented in the literature (Brantferger 1991; Michelsen 1999; Michelsen & Mollerup 2004) and will not be discussed in detail here. There are, however, a couple of key observations about the structure of the minimization that will play a role later in the development of the sequential discretization. First, the Hessian of the negative entropy is a rank-one perturbation of the Hessian of the Gibbs free energy, which is minimized in an isothermal flash calculation. Consequently, it is computationally simple to compute one of these matrices given the other. Second, at equilibrium, the chemical potentials are equal, i.e.(2.6)Here, the chemical potentials are typically specified in terms of mole fractions rather than moles; however, we can use the definition of mole fractions to express them in this form.

In addition to determining the composition of the phases and the temperature, the phase behaviour also determines the properties of the phases. In particular, given pressure, temperature and component mass densities, we can compute the volume occupied by the phases. To complete the mathematical formulation of the system, we require that the sum of the phase volumes matches the available pore volume, which we represent as(2.7)This equation plays the role of an equation of state that constrains the evolution of the component conservation and energy equations. Here, we have implicitly used the capillary pressure to relate the phase pressures to the reference pressure.

## 3. The sequential formulation

In this section, we present a sequential formulation of the above thermal–compositional model. This approach, based on a total velocity splitting with a pressure equation determined by differentiating the equation of state (2.7), was first introduced by Acs *et al*. (1985). The methodology introduced here follows the development in Trangenstein & Bell (1989*b*) and, as such, represents a generalization of that work to a non-isothermal setting. As in Brantferger (1991), we use pressure, total enthalpy and molar densities as the primary unknowns. The use of enthalpy as a primary unknown, instead of temperature, reflects its role as the conserved variable in the energy equation. This choice also avoids issues related to the Gibbs phase rule that can occur if there are more phases than components, in which case pressure and temperature cannot vary independently. However, it also introduces some complexity because many of the thermodynamic variables are given in terms of the temperature, *T*. In particular, the functional form of the enthalpy is given by(3.1)We view this relationship as implicitly defining the temperature as a function of the molar densities, the pressure and the enthalpy. We also note that, for some of the analysis, it will be convenient to derive an equation for temperature; however, in the numerical method, the enthalpy is the quantity that is advanced in time.

Formally, the dynamics evolves conservation equations for the chemical components and energy with phase velocities given by Darcy's law. The entire evolution is constrained by the equation of state.

### (a) Pressure equation

As in Acs *et al*. (1985) and Trangenstein & Bell (1989*b*), we first derive a pressure equation to satisfy the total volume balance by taking the first-order Taylor expansion in time of the pore volume constraint (equation of state)(3.2)Using the functional dependence of *U* on ** n**,

*p*and

*H*

_{t}, we can express the constraint as(3.3)The pressure equation expresses how the pressure needs to be changed to enforce the equation of state (2.7). Owing to splitting errors,

*U*is not necessarily unity at time

*t*and the right-hand side includes a forcing term that attempts to correct errors from previous times. It is more natural to rewrite equation (3.3) in terms of (

*T*,

*p*,

**) since many of the thermodynamic quantities are explicit functions of (**

*n**T*,

*p*,

**). Relationships between the partial derivatives can be derived by inverting the Jacobian matrix of the change of variables from (**

*n**H*

_{t},

*p*,

**) to (**

*n**T*,

*p*,

**). This transformation is described in detail in appendix A. Using this transformation, we can rewrite equation (3.3) as(3.4)where we have dropped the indices that indicated which variables are constant in the partial derivatives since it is always some combination of**

*n**p, T*and

**.**

*n*If we use the evolution equations for the enthalpy and molar densities to substitute for the time derivatives on the right-hand side of (3.4), we obtain the pressure equation(3.5)If we now express the phase velocities in terms of the pressure gradient, we obtain(3.6)This expresses the pressure equation in a form that is formally parabolic. A more detailed analysis of this equation is given below.

### (b) Component and energy conservation equations

In the incompressible case where *ϕ* and *U* are independent of pressure, it is possible to define the notion of a ‘total velocity’ that allows us to decompose the dynamics into an elliptic pressure equation and, ignoring diffusion, a system of hyperbolic conservation laws. To preserve this behaviour in the more general compressible case, we again split the equations by computing the total velocity *v*_{T}. We can then re-express the phase velocities in terms of the total velocity and eliminate the explicit dependence of the phase velocities on the pressure gradient in the conservation equations. First, we define *v*_{T} as(3.7)Then, solving for ∇*p* in terms of *v*_{T}, we express the phase velocity in terms of the total velocity,(3.8)where is the total mobility. Finally, writing the conservation equations in terms of the total velocity yields the fractional flow form of the chemical component conservation equations(3.9)and a fractional flow form of the energy conservation equation(3.10)In this form, the component conservation and energy equations form a system of nonlinear advection–diffusion equations, where the diffusion is typically small relative to the advection.

## 4. Analysis of the sequential formulation

In the limit of no capillary pressure, diffusion or reactions, the sequential formulation discussed above splits the dynamics into a formally parabolic pressure equation and a system of first-order conservation laws for the component conservation and energy equations. For this decomposition to be well-posed, we need to show that the pressure equation is, in fact, parabolic and show that the conservation laws for components and energy form a hyperbolic system. The wave structure associated with the hyperbolic system is also useful in the construction of high-resolution discretizations.

### (a) Thermodynamic derivatives

The structure of the equations is tightly coupled to the phase behaviour. Consequently, we need to understand how perturbations in the state variables are reflected in changes in the phase behaviour. In this section, we derive a number of useful relationships that will be needed later in the analysis. For a more detailed discussion of these derivations, see Trangenstein & Bell (1989*b*) and Brantferger (1991).

As in the case of the derivation of the pressure equation, the analysis is simplified by treating perturbations of the system in terms of *p*, ** n** and

*T*rather than

*p*,

**and**

*n**H*

_{t}. When

*T*is fixed instead of

*H*

_{t}, phase equilibrium perturbations in

**correspond to minimization of the Gibbs free energy rather than negative entropy. However, the similarity in the structure of the different types of equilibrium computations makes it easy to obtain the necessary structure about the isothermal flash and no additional minimization need to occur.**

*n*From the phase equilibrium, we defineandFrom the Gibbs–Duhem equation, *M* is symmetric, *M**n*_{l}=*n*_{l} and *M**n*_{v}=0.

Since the phase potentials are equal at equilibrium, we can differentiate (2.6) with respect to *T* to obtain(4.1)whereSimilarly,where *v*_{α} are the phase-specific volumes.

Other useful thermodynamic relationships are a consequence of first-degree homogeneity. A thermodynamic property is homogeneous of first degree if a multiplication of the composition by a scalar changes the property by the same scalar. Such quantities are typically expressed in the formand the first-degree homogeneity impliesThe first-degree homogeneity holds for enthalpy, phase volumes and phase-specific volumes.

### (b) Parabolicity of the pressure equation

Using these relationships, we can now analyse the pressure equation and the component conservation equations. Certain technical assumptions, which are typically true, are needed for the analysis. These additional assumptions will be noted in the discussion. To analyse the pressure equation, we need to first show that the coefficient multiplying in (3.6),is positive.

Mechanical and thermal stability of the fluid guarantee thatandWe expect both and to be non-negative. Thus, for the validity of the present model, we need to assume that(4.2)It is mentioned (e.g. O'Connell & Haile 2005) that the thermal expansion coefficient is, in general, positive for gases and liquids, except for water under 4°C and at atmospheric pressure. For low-density gases, , and it decreases with increasing pressure. For liquids, *α* has values of an order of magnitude smaller than 1/*T*, and they are nearly constant over modest changes of temperature and pressure. One can also derive, using the Maxwell relations, that , and thus for most liquids and gases and .

Using the first-degree homogeneity, the coefficient of can be rewritten assince , we also need this term to be positive or dominated by the first term. *U* is approximately 1, but there is not enough information to deduce the sign of the r.h.s. term. Thus, we need to assume that this term is positive as well, or that the medium is incompressible, or that the compressibility is sufficiently small to guarantee the parabolicity of the pressure equation.

We must also show that the coefficient of the second-order term on the right-hand side of (3.6),(4.3)defines a second-order elliptic equation. First, we note thatwhich is a lower order term that does not affect parabolicity. Thus, the leading order term on the right-hand side of the equation isThus, the leading term on the right-hand side of the pressure equation is ∇.*λ*_{T}∇*p*, so that right-hand side of the pressure equation (3.6) is elliptic.

### (c) Hyperbolic structure

If we ignore the capillary pressure terms, the fractional flow form of the equations forms a nonlinear system of conservation laws for the enthalpy and the molar densities. The capillary pressure terms introduce a nonlinear diffusion term that is typically fairly small so that the behaviour of the system is typically transport dominated. For incompressible two-phase, two-component systems, the fractional flow form of the equations reduces to the familiar Buckley–Leverett equation. For a solution and discussion, see Collins (1961). The work of Trangenstein & Bell (1989*a*,*b*) shows that for both the three-phase black-oil model and a two-phase compositional model, the conservation laws resulting from a total velocity splitting form a hyperbolic system in the limit of vanishing capillary pressure (subject to well-posedness conditions on the three-phase fractional flow modelling in the black-oil case). In this section, we demonstrate that the total velocity splitting leads to hyperbolicity in the compositional–thermal model as well. As in the case of isothermal compositional flow, the hyperbolicity is directly linked to thermodynamic consistency of the phase behaviour.

In this analysis, we look at the component conservation equations (2.2) and the energy equation (2.3) in the absence of capillary pressure and other diffusion terms, as well as reactions. In the sequential formulation, we view the pressure and total velocity as specifying a prescribed spatial dependence of the flux. In one space dimension, we can write the conservation equations in the form(4.4)and(4.5)whererepresent the phase velocities and l.o.t. denotes the lower order terms that arise from the dependence of the flux on the spatial coordinate either directly or through the spatial dependence on pressure and total velocity. These lower order terms do not affect hyperbolicity.

To show that the system is hyperbolic, we need to show that the eigenvalues of the linearized system, which correspond to wave speeds, are real. Although the dependent variables are the molar densities and enthalpy, it is more natural to perform the characteristic analysis in terms of the molar densities and the temperature. The characteristic structure can also be used to develop numerical discretization for the conservation laws.

When we recast the equations in terms of ** n** and

*T*, we need to be able to perturb the species at constant temperature, not at constant enthalpy. As noted above, having already determined temperature from an isenthalpic flash calculation, we note that an isothermal flash at that specified temperature will produce the same phase compositions and will also satisfy equality of the phase potentials.

We can now linearize equation (4.4), ignoring lower order terms, to obtain(4.6)

Similarly, linearization of the enthalpy equation without lower order terms leads to(4.7)where we have used homogeneity of degree 1 of the partial molar enthalpies. We now substitute for in equation (4.7) to obtain(4.8)whereHere, we have used *M**n*_{l}=*n*_{l} and *M**n*_{v}=0. (The definition of *c*_{p} above does not correspond to a standard definition of specific heat at constant pressure because it includes a mass scaling.)

The system can now be rewritten in the following form:(4.9)whereHyperbolicity depends on *A* having real eigenvalues. There are two cases to consider, one in which only one phase is formed and a second in which both phases are formed. Before discussing these cases, we note that in specifying the system, we have overdetermined the system. One could, for example, eliminate one of the molar densities from the system and solve for the missing variable from the equation of state using the pressure, temperature and other densities. However, owing to the splitting errors, this approach is not conservative. By carrying all of the molar densities and *H*_{t}, we are carrying some redundant information. This redundancy manifests itself as a fictitious wave that essentially only carries information about the consistency of the over-specified system. In the present case, following Trangenstein & Bell (1989*b*), we have defined the decomposition so that the fictitious wave moves with zero speed.

#### (i) Single-phase fluid

For the case of a single-phase flow, which we arbitrarily pick to be liquid here, the matrix *A* takes the form(4.10)where *n*_{l}=** n**. The upper left-hand corner is a rank-one perturbation of the identity matrix. This corresponds to a projection onto the space orthogonal to

**since , which follows from the homogeneity of degree 1 of the phase volumes. This corresponds to the fictitious wave that reflects the redundancy in the equation. Thus,**

*n**A*has a single eigenvalue of 0 corresponding to a right eigenvector (

**, 0) and**

*n**N*−1 eigenvalues of corresponding to right eigenvectors (

*n*^{⊥},0), with . Given the scaling in (4.9), this corresponds to wave speeds of .

The eigenvalue corresponding to the lower right-hand corner of *A* is with eigenvector (*α**n*_{l},1)^{T}, where . Given the scaling in (4.9), this corresponds to a wave speed of . This implies that the thermal wave will be lagged as a result of the heat capacity of the medium. We also note that in the case in which the fluid and the porous medium are not assumed to be in thermal equilibrium, the thermal wave speed is .

#### (ii) Two-phase fluid

The characteristic analysis of the two-phase flow case requires showing that *A* can be transformed into a symmetric matrix. This is not straightforward, but the derivation can be considerably simplified by noting that, since we are doing the characteristic analysis in terms of ** n** and

*T*, the upper left-hand corner of

*A*corresponds to the quasilinear form of the molar conservation equations in the compositional model analysed by Trangenstein & Bell (1989

*b*). This motivates the following similarity transformation of

*A*:where

*R*

_{M}is the matrix of right eigenvectors of

*M*so thatand

*Λ*

_{M}is the diagonal matrix of eigenvalues.

Following Trangenstein and Bell, *R*_{M} can be defined asIn addition, we know that , andWe note that the eigenvectors corresponding to the eigenvalues in represent perturbations in composition, *δ**n*_{l}, that are split between liquid and vapour phases, such thatwhereStability of the mixture with respect to composition implies 0≤*λ*_{k}≤1, so that all of the entries in *Λ*_{M} are between 0 and 1. Applying this similarity transformation to *A* yields

Finally, *Ã* can be rewritten as(4.11)whereand *P*_{R}=[**0**,**0**,*I*] is a ((*n*_{c}−2)×*n*_{c}) projection operator.

In Trangenstein & Bell (1989*b*), it was shown that there is a similarity transformation that diagonalizes *V*, and its eigenvalues are real. One of these eigenvalues is zero and corresponds to the fictitious wave. The other is , where *S*_{v}=*u*_{v}/(*u*_{l}+*u*_{v}) is the vapour saturation. This corresponds to the Buckley–Leverett-type wave. We now need to prove that there is a similarity transformation that symmetries the sub-matrix,from which we conclude that the eigenvalues of this sub-matrix should also be real. Thus, all eigenvalues of *Ã* would be real and the hyperbolicity is proved.

Substituting into the expression for *d*^{T}, we obtainNext, we make use of the following result, which can be derived by using several similarity transformations described in Trangenstein & Bell (1989*b*):where *D* isThis result is then used to further rewrite as*Ã*_{sub} has now becomewhere and . Since the matrix is non-singular, we can use it to define a final similarity transformation,which results in(4.12)*Â*_{sub} is a symmetric matrix, therefore it has real eigenvalues and the system of linearized conservation laws is hyperbolic. Although we are not including capillary pressure, thermal conduction or other diffusive processes in the analysis, these effects can be important in some applications. If included, these effects appear as second-order diffusive terms in the conservation laws that can be easily modelled using standard finite-difference techniques. The analysis here focuses on the limiting behaviour as these terms become small. Including these terms would serve to smooth the nonlinear hyperbolic waves associated with the underlying hyperbolic system. This completes the derivation of hyperbolicity for the two-phase flow.

## 5. Discretization issues

In this section, we present a numerical method for solving the thermal–compositional model described above. The intent here is to use the method to elucidate some of the wave phenomena associated with this type of system. For this reason, we focus on a simplified version of the model. We consider only one-dimensional flows without gravity. We also assume that there is no capillary pressure or diffusive transport and that there are no reactions. With these assumptions, we solve the conservation laws and the pressure equation on a uniform grid using a finite volume approach. A mesh is defined in the (*x*,*t*)-plane. The computational domain is restricted to *x*∈[0,*L*]. The points on the mesh are at locations (*x*_{i}=*i*Δ*x*,*t*^{n}=*n*Δ*t*), with *i*=0, …, *N*_{x} and *n*=0, …, *N*_{t}. The discrete values of at (*i*Δ*x*+Δ*x*/2,*n*Δ*t*) will be denoted by .

Pressure, enthalpy and component molar densities are defined on cell centres and the total velocity is defined on cell faces. In outline semi-discrete form, the algorithm for solving the conservation laws for components and enthalpy, equations (3.9) and (3.10), and the pressure equation (3.6) is as follows.

Given a solution at , calculate the phase equilibrium and hence .

Solve the pressure equation implicitly,

Calculate

*v*_{T}using*p*^{n+1}.Solve component and enthalpy conservation equations,and

Set (

*p*^{n},*H*^{n},*n*^{n})=(*p*^{n+1},*H*^{n+1},*n*^{n+1}) and go to step (i).

The different steps are described in more detail as follows.

*Step* (i). For the phase equilibrium calculation, an isenthalpic flash algorithm developed by Michelsen (1999) is used. The algorithm currently considers the fluid enthalpy, not the total enthalpy, in the minimization of negative entropy, therefore we have incorporated an outer iteration to solve the full minimization problem. In particular, we added an outer iteration to find *T*, such that(5.1)where *H*_{f} is at the minimum negative entropy. This nonlinear equation is solved using the secant method. Details of this additional step are described in appendix B.

*Step* (ii). The pressure equation is treated implicitly, thus, using a central difference approximation, equation (3.5) becomesQuantities at the cell faces, i.e. , are calculated as an average of their respective cell centre values , hence no extra phase equilibrium calculation is needed. The discretized system is solved with a tridiagonal matrix solver.

*Step* (iii). To calculate *v*_{T}, the variables at time *n* are used, except for calculating the pressure gradient in the definition of *v*_{T}, for which the pressure at time *n*+1 is used. Consequently, the total velocity is discretized as

*Step* (iv). This is a system of hyperbolic equations and we have used a simple first-order (in time and space) conservative upwind scheme to solve it since our numerical examples are designed in such a way that all the eigenvalues are positive. Hence, the system of conservation lawswhereis discretized using the first-order Godunov approachwhere strictly upwind fluxes are used,

The above forms a basis for higher order discretizations; we will present appropriate methodologies in a future communication.

## 6. Numerical examples

In order to give a more quantitative description of the types of waves that can form in a realistic situation where flow velocities depend on pressure gradients, four test cases are calculated. They correspond to particular situations occurring during thermal recovery processes: hot gas; hot liquid; and hot two-phase mixture injection.

In all four test cases, the following parameters are used: number of grid points *N*_{x}=400; reservoir length *L*=7620.0 (cm); permeability *k*=2×10^{−8} (cm^{2}); relative permeability ; rock reference temperature ; liquid viscosity *η*_{l}=0.001 (Pa s); and vapour viscosity *η*_{v}=6×10^{−5} (Pa s). The values for other parameters such as binary interaction parameters and critical properties, etc. used in the Peng–Robinson equation of state can be found in Danesh (1998).

It is sometimes necessary to dampen the influence of the pressure correction term in the pressure equation (3.6) by a factor *f*_{press}, e.g. use *f*_{press} (1−*u*_{l}−*u*_{v}) for the correction term. Normally, we take *f*_{press}=1. In table 1, the values for pressure correction factor *f*_{press}, time step Δ*t* and simulation time *t* are given.

The eigenvalues can be calculated in two ways. First, one can use the matrix *A* defined in equation (4.9), rescaled by *ϕ*, *c*_{p}. The other option is to use matrix *Ã* in equation (4.11). The eigenvalues of matrix *V* are zero and the Buckley–Leverett wave, which is . Subsequently, only the eigenvalues of the matrix *Â*_{sub} in equation (4.12) are needed. The second system is a two-by-two system and can easily be solved analytically.

In the single-phase case, there are *N*−1 essentially linearly degenerate waves, with wave speed , and one linearly degenerate wave, with wave speed zero. The ‘energy’ wave with wave speed is also nearly linearly degenerate, but propagates at a reduced speed. If *ϕ*=1, then the thermal wave speed is also .

In figures 1–4, a number of quantities (component densities, total enthalpy, pressure, eigenvalues, total velocity, saturation and deviation of the equation of state from unity) are presented as a function of position in the reservoir, for test cases 1–4, respectively. The error in the equation of state, which takes the form of a volume discrepancy, remains small in all cases as the calculations show.

In test 1, a hot vapour containing 95 per cent methane and 4.9 per cent butane, 0.1 per cent nonane is injected into a liquid containing 20 per cent methane, 20 per cent butane and 60 per cent nonane. The details of the initial and boundary conditions of the reservoir are summarized in table 2. The results are plotted in figure 1. The solution structure for this problem consists of two discontinuous waves that connect the single-phase vapour injection mixture to the single-phase liquid initially in the reservoir. The slower wave corresponds to the transition from a single-phase region to a two-phase region. At these transition regions, the flux function is not differentiable and there are discontinuous changes in the eigenstructure. This type of structure is fundamentally different from the behaviour of Buckley–Leverett owing to mass transfer effects between phases and, consequently, we do not see the long rarefaction associated with the simpler system. The faster discontinuity is a shock wave from the *λ*_{1} family.

In test 2, the same configuration as in test 1 is considered, but now with *ϕ*=0.4. The results are plotted in figure 2. The introduction of the rock heat capacity results in significant changes in the flow field. The extent of the two-phase region is much larger and shows a typical Buckley–Leverett rarefaction-shock pattern in *λ*_{1} for *x*>2000 (cm). Again, we see a discontinuous wave separating the single-phase vapour region from the two-phase region; however, in this case, it is considerably slower. Within the two-phase region, we also see an additional weak wave at approximately *x*=1600 (cm). The interaction between the waves inside the two-phase region is also more complex. Inside the two-phase region, the compound wave consists of a constant state, a contact discontinuity and a rarefaction wave. We also note that *λ*_{1} crosses the other waves at approximately *x*=1000 (cm), indicating a loss of strict hyperbolicity.

In test 3, a hot liquid containing 10 per cent carbon dioxide and 20 per cent butane, 70 per cent nonane is injected into a liquid containing 80 per cent carbon dioxide, 10 per cent butane and 10 per cent nonane, and the details of the initial reservoir and boundary conditions are summarized in table 3. The results are plotted in figure 3. This is a single-phase liquid simulation. Owing to the presence of the porous medium, there are two distinct wave speeds. This results in a slow wave in the *λ*_{3} at approximately *x*=1500 (cm) separated by a constant state and followed by a faster wave at approximately *x*=1500 (cm). The faster wave corresponding to *λ*_{1}, *λ*_{2} is essentially a contact discontinuity. Curiously, the slower wave is also essentially linearly degenerate. The structure in this case is carried by the change in the total velocity and reflects a change in the density resulting from changes in the compressibility as a function of temperature.

Test 4 is an example of a hot two-phase mixture containing 40 per cent carbon dioxide and 10 per cent butane, 50 per cent nonane being injected into a colder two-phase mixture containing 80 per cent carbon dioxide, 10 per cent butane and 10 per cent nonane. The details of the initial reservoir and boundary conditions are summarized in table 3. The results (figure 4) show a rather complex behaviour. The hot two-phase mixture first condenses in a small region ahead of the propagating hot front. This is a result of the lagged behaviour of the thermal wave arising due to the heat capacity of the rock. The transition again shows the compound wave behaviour typical of Buckley–Leverett. The mixture then transitions back into a two-phase mixture with a transition between phases accompanied by another Buckley–Leverett shock. This multiplicity of waves is made possible by the loss of strict hyperbolicity shown by the crossing of the eigenvalue curves.

## 7. Summary and conclusions

A system of equations describing the flow of non-isothermal multicomponent two-phase fluids in a porous medium has been developed and analysed. The equations governing the system are an extension of the compositional solver developed in Trangenstein & Bell (1989*b*). It is shown that the system of equations can be split into a hyperbolic system of conservation laws for the component density and the enthalpy and a parabolic pressure equation that constrains the volume of fluids to the available pore volume. This procedure of decoupling the fundamental equations forms the basis of a sequential numerical algorithm. The sequential method is then applied to study four test cases of multiphase flows that can be encountered in thermal recovery processes. They demonstrate that the method can resolve the interacting waves present in the flow field and can shed new light into the understanding of the solution structure. The analysis presented here can also form the basis for the development of higher-order approaches for the component and energy conservation equations and for including additional physical phenomena such as capillary pressure, diffusion and reactions.

## Acknowledgments

D.E.A.v.O. was funded by Schlumberger Cambridge Research. The work of J.B.B. was supported by the Director, Office of Science, of the US Department of Energy under contract DE-AC02-05CH11231.

## Footnotes

- Received July 1, 2008.
- Accepted October 6, 2008.

- © 2008 The Royal Society