## Abstract

This paper is concerned with multicomponent, two-phase, thermal fluid flow in porous media. The fluid model consists of component conservation equations, Darcy's law for volumetric flow rates and an enthalpy conservation equation. The model is closed with an equation of state and phase equilibrium conditions that determine the distribution of the chemical components into phases. The sequential formulation described in a previous article is used to build a second-order shock capturing scheme for the conservation equations using a primitive-variable-based linear reconstruction. The fluxes at the cell faces are calculated using an approximate Riemann solver. The method is validated and evaluated by means of one- and two-dimensional problems, including a gravity inversion test.

## 1. Introduction

In our previous work (van Odyck *et al.* 2009), a system of equations was derived to describe thermal flow in porous media. This system has applications to a variety of thermal oil recovery processes such as steam flooding and *in situ* combustion. Relevant work in this area can be found in Liu *et al.* (2007), Huang *et al.* (2007), Pasarai *et al.* (2005), Nilsson *et al.* (2005), Christensen *et al.* (2004) and Kristensen *et al.* (2009). Soil remediation is another important field in which thermal compositional solvers play a role (Class & Helmig 2002).

Many approaches to thermal simulations, including most of those cited earlier, are based on a fully implicit formulation with lower-order implicit spatial differencing. This approach offers significant stability benefits. However, here we are interested in capturing the sharp moving fronts that characterize many of the earlier-mentioned physical processes. For this reason, we explore a semi-explicit sequential approach, similar to the method proposed in Bell *et al.* (1986) for reducing numerical dispersion in a black oil model.

To capture fronts accurately, high-resolution shock-capturing schemes might be employed (Toro 1999; LeVeque 2002). Previous work using these types of methods for thermal simulation can be found in Trangenstein & Bell (1989) and Mallison *et al.* (2005). A further advantage of these types of methods is that they naturally lend themselves to implementation in adaptive mesh refinement (AMR) frameworks; see the recent work by Pau *et al.* (2012).

It was shown in van Odyck *et al.* (2009) that the system can be split into a hyperbolic part and a parabolic part, and a first-order upwind scheme was used to solve the hyperbolic system. In this work, we extend the order of accuracy of the hyperbolic solver to nominal second-order, using a MUSCLHancock primitive variable scheme, where MUSCL stands for monotonic upstream-centred scheme for conservation laws (Toro 1999). Because primitive variables are used for the second-order reconstruction, this scheme has the advantage of avoiding costly phase equilibrium calculations at cell interfaces.

The hyperbolic part of the system under consideration has highly nonlinear flux functions that exhibit eigenvector deficiencies and local linear degeneracies. The polymer flood model (Pope 1980; Isaacson 1989) exhibits the same kind of degeneracies and is used as a benchmark against which to test our new numerical method. Existing methods which can handle eigenvector deficiencies, such as the method described in Bell *et al.* (1989), are difficult to use in the case of the equations for thermal flow in porous media. These methods are based on the ability to monitor the change in value of eigenvalues of distinct wave families. In the case of the equations for thermal flow in porous media, there is no general analytical expression for the eigenvalues. It is therefore difficult to monitor the change of an eigenvalue of a certain wave family from one computational cell to the next.

A standard way of constructing a second-order scheme is to reconstruct conserved variables at cell faces. In the present case, using the conserved variables in this manner leads to a large amount of computationally expensive phase equilibrium calculations. A way around this is to reconstruct the primitive variables instead. To this end, the system of equations describing conservation of molar density and enthalpy is reformulated in primitive variables. This approach is demonstrated on two gas-injection test cases. In addition, to show that the method can handle waves travelling in opposing directions, a gravity inversion test is introduced wherein a cold multicomponent liquid is placed on top of a hot multicomponent vapour. The two phases tend to invert under the influence of gravity. The change of composition has an influence on the phase equilibrium and a mixture is created. In two dimensions, viscous fingering and complex phase transitions are observed. Similar instabilities have previously been studied in, for example, Tan & Homsy (1988), Rogerson & Meiburg (1993) and Islam *et al.* (2007).

The paper is structured as follows. In §2, the basic shock-capturing method (a primitive-variable MUSCL–Hancock scheme) is described. Section 3 applies this method to a model problem (the polymer flood model) which is chosen because it is simpler than the full thermal model but still shares similar numerical difficulties. In §§4 and 5, the thermal model, our main problem of interest, is introduced and a solution method developed. Section 6 derives the primitive formulation needed to apply the MUSCL scheme to the thermal model. Section 7 presents an overview of the numerical method, and results are shown in §8. Finally, §9 contains a summary of the paper.

## 2. Description of numerical method

In this section, we describe the basic numerical method (a primitive-variable-based MUSCL–Hancock scheme) that is used to solve the hyperbolic part of the flow equations. This method will first be applied to a model problem: the polymer flood model (Isaacson 1989). This model problem contains the same difficulties as encountered in the porous media flow equations, namely eigenvector deficiencies and local linear degeneracies. We will then apply the method to the full porous media flow equations, which will require the derivation of a similarity transform between the conserved and primitive variables of that model.

Primitive variables are chosen in the reconstruction because of the nature of the porous media flow model (see §4). If instead the conserved variables were used, extra phase equilibrium calculations would be needed at cell faces to convert the (reconstructed) conserved variables back to the primitive variables that are needed to calculate fluxes. The use of this hybrid scheme therefore represents a significant computational saving over the use of a scheme that reconstructs conserved variables.

In general, we would like to solve the system of conservation laws
2.1
for the conserved variable **U**. We are interested in using a finite-volume update formula
2.2
where approximates the integral average of the solution in cell *j* centred at *x*_{i} at time *t*^{n} and **F**_{i+1/2} is the numerical flux at the cell face *x*_{i+1/2}. The flux at the cell face is calculated by solving a Riemann problem using the Harten, Lax, van Leer (HLL) solver (Toro 1999):
2.3
with **F**_{L}=**F**(**U**_{L}) and **F**_{R}=**F**(**U**_{R}). Here, L (R) indicates the state in the cell to the left (right) of the cell face *x*_{i+1/2}. The highest left and right wave speeds are calculated using
2.4
where the *λ*_{k} are the eigenvalues of the flux Jacobian ∂**F**/∂**U**.

Second-order accuracy is achieved by means of a MUSCL approach (see Toro 1999), using the non-conservative form of the equation, namely
2.5
to reconstruct the state at the cell face. Here **W** is the vector of primitive variables and **A**=(∂**U**/∂**W**)^{−1}(∂**F**/∂**U**)(∂**U**/∂**W**). The primitive variables are linearly extrapolated to the cell faces using limited slopes, and these extrapolated variables are evolved for half a time step. This gives the states to the left and the right of the cell face *x*_{i+1/2} as
2.6
where the Δ_{i} are limited slopes. These can be calculated in various ways. For our calculations, using the polymer flood model, we use
2.7
with (Toro 1999, eqn. 14.44). The choice *β*=1 represents the minmod slope limiter. For the thermal flow calculations, we use a van Leer-type slope limiter (Toro 1999, eqn. (14.49) and (14.54)).

Finally, the flux **F**_{i+1/2} for the conserved variables is computed with the HLL Riemann solver using the evolved variables **W**_{L} and **W**_{R} (instead of **U**_{L}, **U**_{R}) as initial data. Note that although primitive variables are used for the extrapolation and evolution steps, the final update is performed on the conserved variables. Therefore the scheme is conservative.

## 3. Polymer flood model

This section demonstrates the application of the MUSCL–Hancock primitive variable scheme to a model problem: the polymer flood model. The aim of this section is to demonstrate that the scheme can cope with challenging numerical features such as local linear degeneracies and eigenvector deficiencies. These features are shared with the compositional equations for thermal flow in porous media that are developed in subsequent sections.

The polymer flood model is used to describe the addition of a polymer into water, which is injected into a reservoir. The polymer increases the viscosity of the injected water and therefore enhances sweep efficiency and oil recovery. This model is described by Pope (1980) and analysed by Isaacson (1989). It is used in Bell *et al.* (1989) as a model system to show that their second-order scheme can handle non-strictly hyperbolic PDEs.

The equations for thermal flow in porous media contain eigenvector deficiencies and local linear degeneracies, features shared by the polymer flood model. If the numerical algorithm converges to the exact solution of the polymer flood Riemann problem then this is a good indication that the method will perform well for the porous media flow model (for which there is no exact Riemann solution). The test cases in this section are chosen to match those in Bell *et al.* (1989, figs 1–4).

The polymer model describes the flow in a porous medium of a two-phase incompressible fluid consisting of oil, water and polymer. The model takes the form of two conservation laws with
3.1
where *s* is the volume fraction of the aqueous phase, *c* is the concentration of polymer in the aqueous phase and *v*(*s*,*c*) is the phase velocity of water. The primitive variables are
3.2
The eigenvalue matrix and the right-eigenvectors of the Jacobian ∂**F**/∂**U** are, respectively,
3.3
The system is not strictly hyperbolic. If *λ*_{1}=*λ*_{2} then there is an eigenvector deficiency and *R* is singular. The 1-wave can have a local linear degeneracy:
3.4
The 2-wave is linearly degenerate:
3.5
The phase velocity of water *v*(*s*,*c*) is modelled as follows (Bell *et al.* 1989):
3.6
Here, *v*_{t} and *g* are constants representing, respectively, the total flow and acceleration due to gravity. All variables are dimensionless.

The first test we consider is a gravity inversion problem in which a heavy fluid is placed above a light fluid. We take *v*_{t}=0 and *g*=−1.4. The initial data are *c*=0.1,*s*=0.1 for *x*<0 and *c*=0.9,*s*=0.95 for *x*>0. The simulation is performed with 400 cells and 400 time steps using the primitive-variable MUSCL scheme with a fixed Δ*t*/Δ*x*=0.25 and minmod limiter. The results are shown in figure 1.

The analytic solution (solid line) consists of a shock moving to the left, separated by a constant state from a compound wave consisting of a contact discontinuity moving to the left, a transonic rarefaction and a shock moving to the right. There is a loss of strict hyperbolicity at the point where the eigenvalues collapse, at the contact discontinuity. There is a good agreement between analytic and numerical solution: shock positions are correct and well-resolved, and compound waves are captured correctly. In figure 1*b*, the spike in the numerical wave speeds at *x*/*t*=1.1 is a characteristic signature of a local linear degeneracy; it arises because the smeared shock passes through a region of phase-space in which the wave speed is higher than the shock speed (Bell *et al.* 1989).

The second problem we consider is designed in such a way that the solution contains a stationary eigenvector deficiency. The computations are performed in a moving frame of reference; both wave speeds are zero at *x*/*t*=0. For this problem, we take *v*_{t}=1, *g*=0. The solution produced with the use of the HLL flux does not satisfy the Oleinik–Liu entropy condition in the neighbourhood of the stationary eigenvector deficiency. Therefore, we add an artificial viscosity term to the HLL flux,
3.7
Here *ν*=0.1 was used. The initial data are *c*=0.9,*s*=1 for *x*<0 and *c*=0.1,*s*=0.1 for *x*>0. The simulation is performed with 400 cells and 400 time steps with Δ*t*/Δ*x*=0.5. The results are shown in figure 2. The analytic solution (solid line) consists of a shock moving to the right, separated by a constant state from a compound wave consisting of a sonic contact discontinuity and a rarefaction wave moving to the left.

Finally, the grid refinement study in figure 3 shows convergence for both polymer test cases considered. We have included a first-order result for comparison.

## 4. Mathematical formulation of thermal flow in porous media

In this section, we move on to the main problem of interest, and present equations describing non-isothermal, multiphase, multicomponent flow in heterogeneous porous media. A solution method for this model is developed over subsequent sections.

The porosity of the medium is denoted by *ϕ*, and the phase volume of each phase per pore volume is denoted *u*_{α}. (Here, Greek subscripts refer to mobile phases. The medium itself, which can be viewed as a solid phase, is treated separately because 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 molar density in the combined fluid system. (Mineral content of the solid is, of course, important for geochemistry but is not included in these definitions.) Table 1 summarizes the main physical variables of importance for this model, together with their units.

### (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*_{α} as
4.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 respectively the phase viscosity and density, and ** g** is the acceleration due to gravity. Here

*λ*

_{α}≡

*Kk*

_{r,α}/

*η*

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

*p*, typically taken to be that of the most wetting phase, by a capillary pressure

*p*

_{c,α}=

*p*

_{α}−

*p*, which itself is a function of saturation.

Conservation of mass for each component is given by
4.2
where *u*_{α} is the volume occupied by phase *α* divided by the pore volume at thermodynamic equilibrium. 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.

To model non-isothermal effects, it is necessary to include 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 total energy balance is of the form
4.3
where is the total enthalpy of the system, **h**_{α} the partial molar phase enthalpies, *H*_{r} the enthalpy of the medium and *ρ*_{r} 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 are in thermal equilibrium, then

*R*

_{H}will contain relaxation terms that equilibrate the two temperatures.

### (b) Phase behaviour

The component conservation and energy equations express the change in total mass of each of the component due to advection, diffusion and chemical reactions. Since components are transported in phases, it is necessary to know the distribution into phases before we can solve the flow equations. The local phase equilibrium occurs at a phase composition of maximum entropy, where the entropy of each phase is derived from the chemical potential *μ*_{α}. For a more detailed discussion, see Michelsen & Mollerup (2004) and Brantferger *et al.* (1991).

The chemical potentials *μ*_{α} are functions of *p*_{α}, *T* and phase composition **n**_{α}. The major thermodynamic variables describing each phase can all be expressed in terms of the chemical potential of that phase. In particular, the partial molar entropies are given by
4.4
where are the mole fractions. The partial molar enthalpies are given by
4.5

Here the phase equilibrium problem is to determine the composition of the phases **n**_{α} given the total moles **n**, pressure *p* and total enthalpy *H*_{t}. The equilibrium distribution of the components is given by minimizing the negative entropy of the system. In the two-phase liquid and vapour case, this means minimizing the quantity
subject to
and
along with inequality constraints guaranteeing non-negativity of the compositions (**n**_{α}≥0) and thermal stability of the fluid . The treatment of this minimization problem, the so-called isenthalpic flash calculation, has been discussed extensively in the literature (Brantferger 1991; Michelsen 1999; Michelsen & Mollerup 2004) and will not be discussed in detail here.

To complete the mathematical formulation of the system, we require that the sum of the phase volumes match the available pore volume, which we will represent as
4.6
(Here we have implicitly used the capillary pressure to relate the phase pressures to the reference pressure.) This equation plays the role of an equation of state that constrains the evolution of the component conservation and energy equations. In a numerical method, this constraint will not necessarily be satisfied at every point in the domain; we refer to the value *U*−1 as the volume discrepancy.

## 5. Sequential formulation

In this section, we describe a sequential formulation of the thermal compositional model discussed above. This approach, based on a total velocity splitting with a pressure equation determined by differentiating the equation of state (equation 4.6), was first introduced by Acs *et al.* (1985). Formally, conservation equations for the chemical components and energy are evolved with phase velocities given by Darcy's Law. The entire evolution is constrained by the equation of state. We follow the work in van Odyck *et al.* (2009) and use the same expressions for the pressure equation and conservation laws for molar density and enthalpy.

### (a) Pressure equation

The pressure equation follows by taking the first-order Taylor expansion of the pore volume constraint (equation (4.6)) in time over a single time step of Δ*t*; see §3*a* of van Odyck *et al.* (2009). The partial derivatives are computed in (*T*,*p*,**n**) space:
5.1
(eqn. (3.6) in van Odyck *et al.* 2009). In that paper, it was shown that the pressure equation is formally parabolic subject to physical constraints.

### (b) Component and energy conservation equations

The total velocity
5.2
is introduced in order to eliminate the explicit dependence of the phase velocities on the pressure gradient in the conservation equations. Writing the conservation equations in terms of the total velocity yields the fractional flow form of the chemical component conservation equations:
5.3
where . A fractional flow form of the energy conservation equation is given by
5.4
(cf. van Odyck *et al.* 2009, eqns (3.7), (3.9) and (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.

### (c) 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 repeat a number of useful relationships that will be needed later in the definition of the characteristic structure of the component and energy conservation laws. For a more detailed discussion of these derivations, see van Odyck *et al.* (2009, §4a), Trangenstein & Bell (1989) and Brantferger (1991).

The computations are simplified by treating perturbations of the system in terms of *p*, **n** and *T* rather *p*, **n** and *H*_{t}. When *T* is fixed instead of *H*_{t}, phase equilibrium perturbations in **n** correspond to minimization of the Gibbs free energy rather than negative entropy. From the phase equilibrium, we define
From the Gibbs–Duhem equation, *M* is symmetric, *M***n**_{l}=**n**_{l} and *M***n**_{v}=0.

### (d) Characteristic structure

In van Odyck *et al.* (2009, §4*c*), it was shown that the conservation laws are a hyperbolic system. In this section, we recall the results therein derived concerning the characteristic structure of the system. The molar density and enthalpy conservation laws can be put into the form
5.5
where
5.6
with
5.7
and
5.8
The quantity **r**_{A} is a lower-order term arising from the transformation of variables, rather than a physical source term. Finally, the *σ*_{α} are defined by
5.9
The work in van Odyck *et al.* (2009, §4*c*(ii)) details how to obtain the eigenvalues of *A*. These eigenvalues contain a spurious zero-wave as a result of the splitting, and the Buckley–Leverett wave
Here the phase saturations *s*_{l,v}=*u*_{l,v}/**e**^{T}**u** are the fractions of the fluid volume occupied by the liquid and vapour phases, respectively. The other eigenvalues can be found numerically. Furthermore, it can be shown that there is a similarity transformation between the flow equations expressed in (**n**,*T*) variables and those expressed in (*ϕ***n**,*H*_{t}) variables. Therefore, we can use these eigenvalues in our numerical scheme in which we advance the conserved variables (*ϕ***n**,*H*_{t}).

## 6. Flux computation

To update the variables over a time step, the MUSCL approach is used as described in §2. This requires an expression for the primitive flux Jacobian, which is developed in this section.

The conserved variables are
and the fluxes in the *x*-direction follow from the conservation laws ((5.3) and (5.4):
In two dimensions, similar expressions can be derived for the fluxes in the *y*-direction. This work is focused on modelling systems with no capillary pressure, no reaction terms and no diffusion terms. With some additional work, these terms can be added to the numerical method explained in this article. For example, in Monmont *et al.* (2011), we extend this method with reactions and thermal diffusion to create a model for *in situ* combustion.

### (a) Primitive-variable formulation

In this subsection, we derive a primitive-variable formulation of the conservation laws for use in the second-order MUSCL scheme. The molar densities in the different phases depend on the total molar density, temperature and pressure. Therefore,
6.1
6.2
The above expression gives
6.3
and a similar equation can be derived for **n**_{v}:
6.4
The ∂*p*/∂*t* term can be replaced by terms containing only spatial derivatives with the help of the pressure equation (5.1). This is a second-order derivative in space and does not alter the wave structure. Therefore, the ∂*p*/∂*t* term is cast as a source term
The ∂**n**/∂*t* and the ∂*T*/∂*t* term in equations (6.3) and (6.4) are replaced by terms only containing spatial derivatives with the help of equation (5.5). The system of equations for the primitive variables now reads
6.5
The matrices (*A*_{1},*B*_{1}) are defined as
6.6
and
6.7
the vectors (**a**_{2},**b**_{2},**c**_{1}) as
6.8
6.9and
6.10
with the *A*_{ij} as defined in (5.7) and the scalar *c*_{2} as
6.11
The system (6.5) is used to compute conservative inter-cell fluxes using the second-order method described in §2.

In the MUSCL approach, the primitive values at the cell faces are linear extrapolations from the value at the cell centre. They cannot capture a phase change: if the cell-centred value represents a liquid, then the reconstructed value at the cell face will also be a liquid. Nevertheless, in §8, we show that the method still displays good numerical properties. The advantage of the method is that we only need one phase equilibrium calculation per computational cell per time step.

## 7. Discretization issues

In this section, we present an overview of our numerical method for solving the thermal compositional model described earlier. We 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. The exposition below is for one space dimension; for multiple dimensions, we use a dimensional splitting 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}) with *i*=0,…,*N*_{x} and *n*=0,…,*N*_{t}. The discrete values of quantities ** Q**(

*x*,

*t*) at (

*i*Δ

*x*+Δ

*x*/2,

*t*

^{n}) 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 form, the algorithm for solving the conservation laws for components and enthalpy (equations (5.3) and (5.4)), and the pressure equation (5.1) is as follows:

Given a solution at

*t*^{n}, (*p*^{n},*H*^{n}_{t},**n**^{n}), compute phase equilibrium calculation .Compute time-step Δ

*t*.Solve the pressure equation implicitly.

Calculate

*v*_{T}using*p*^{n+1}.Solve the 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 1.

We now describe the different steps in more detail.

Step 1. 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
7.1
where *H*_{f} is at the minimum negative entropy. This nonlinear equation is solved using a bracketed secant method. Details of this additional step are described in van Odyck *et al.* (2009).

Step 2. A time-step Δ*t* is calculated using the CourantFriedrichsLewy (CFL) condition based on the eigenvalues of the flux Jacobian, and is used for both the parabolic and hyperbolic steps. We have found that this is generally effective. We also limit the maximum growth of the time step. Suppose Δ*t*_{0} is the previous time-step and Δ*t*_{1} is the new time step calculated using the CFL condition. Then we simply replace Δ*t*_{1} with
7.2
where *f*_{Δt} is a constant growth-limiting parameter, which we take to be 1.5. The time step is also multiplied by a constant factor (usually 0.2) if the maximum volume discrepancy in the domain is larger than some value (usually 0.09). The purpose of these measures is to control the volume discrepancy and limit its growth, which improves the simulation stability and accuracy. See for example the work by Dicks (1993), where these techniques were also used.

Additional physical effects such as reaction or diffusion may impose other considerations on time-step selection. For example, stiff reaction mechanisms can lead to large volume changes in reaction zones, placing additional constraints on the time step. However, explicit time steps can still be used if an operator-splitting approach is used for the source terms. This allows for time-step sub-cycling on the source term ODEs, using any standard explicit or implicit ODE solver. We have found that this avoids many potential stability problems.

Step 3. The pressure equation is treated implicitly. Using a central-difference approximation, equation (5.1) becomes Quantities at the cell faces, i.e. , are calculated as an average of their respective cell-centred values, ; hence no extra phase equilibrium calculation is needed. The discretized system is solved with a tridiagonal matrix solver in one dimension, and a biconjugate gradient stabilized method (Van der Vorst 1992) with diagonal preconditioner in two dimensions.

The pressure equation is the result of a first-order Taylor expansion in time of the equation of state (see §5). It gives an approximation for the pressure which up to first-order fulfils the equation of state (equation 4.6). Because of the sequential formulation, the volume discrepancy is non-zero, and crucially this discrepancy is (approximately) compensated for by the (1−*U*)/Δ*t* term on the right-hand side of the pressure equation (which is a result of the Taylor expansion). The discrepancy could be eliminated by iterating between pressure and advection equations. However, this would require many expensive flash calculations and negate the advantages of the sequential method. We find that the time-step controls described earlier are sufficient to control the volume discrepancy.

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

Step 5. This is a system of hyperbolic equations and we use the numerical method described in §2 to solve it. The system of conservation laws
where *Q*=(*ϕ***n**^{T},*H*_{t})^{T} and
is discretized using a finite-volume update formula
where is the HLL flux at time-step *n* using the evolved primitive variables **W**_{L},**W**_{R} as defined in equation (2.6) and the Jacobian matrix **A** as defined in equation (6.5).

Dimensional splitting is used to obtain a numerical solution in two dimensions. In each dimension, the one-dimensional Godunov update (equation 2.2) is used with the appropriate flux function. For convergence in multiple dimensions, a Strang-type dimensional splitting will preserve the order of the one-dimensional advection scheme (Toro 1999). Note that this type of simple dimensional splitting can suffer from grid orientation effects in multiple dimensions (Shubin & Bell 1984). A possible solution is to use unsplit schemes as in Bell *et al.* (1988). Unsplit schemes using MUSCL-type approaches are possible; see for example, the work by Colella (1990) and LeVeque (1997). However, these methods involve solving additional Riemann problems in transverse directions and are thus generally more expensive than the splitting approach.

## 8. Numerical examples

This section presents numerical solutions of the porous media flow model using the primitive-variable MUSCL scheme developed in previous sections. For all computations, the term **r**_{A} in the primitive system (equation (6.5)) is omitted: it appears only in the extrapolation step, where it is constant, and we found in practice that it had a negligible effect. All computations use the physical parameters given in table 3, quadratic relative permeability and the Peng–Robinson equations of state, properties for which can be found in Danesh (1998).

### (a) Convergence properties

To test the convergence properties of the method, we consider a single-phase test case with a smooth solution. A one-dimensional domain of length *L*=1000 cm is initialized with methane (C1), butane (C4) and nonane (C9) with mole fractions 0.7, 0.15, 0.15, respectively. Let *x*′ be the spatial coordinate divided by the domain length *L*. A smooth non-uniform initial temperature distribution (in Kelvin) is specified as follows:
8.1
An inflow (injection) boundary condition is used for the left boundary and a transmissive boundary condition for the right boundary. A pressure gradient is enforced by specifying the pressure as 13.8 MPa on the left boundary and 13.7 MPa on the right boundary. Gravity is neglected and porosity is uniformly 1. The simulation is run until *t*=1.5×10^{4} s using the van Leer slope-limiter with a CFL number of 0.9.

Convergence data for this test are given in table 2. Errors were measured against a reference solution computed on a fine grid of 25 600 cells using the MUSCL scheme. The method displays a convergence rate of approximately 2 for the grid sizes considered. Error norms are also decreased compared with a first-order HLL method. For example, the first-order method gave an L2 error norm of 3.25×10^{−7} for *N*=1600, two orders of magnitude greater than the MUSCL method. However, because the pressure equation is derived using a first-order Taylor expansion in time of the pore-filling constraint (4.6), and is also discretized to first-order in time, we expect temporal errors for the MUSCL method to be first-order for very fine grids.

### (b) One-dimensional test cases

Three one-dimensional tests are now presented. Parameters for these tests can be found in tables 3–6.

We first consider an isothermal test case similar to the first test given in Trangenstein & Bell (1989). This test case simulates the injection of a gas consisting of 95 per cent methane, 4.9 per cent butane and 0.1 per cent nonane into a liquid reservoir consisting of 20 per cent methane, 20 per cent butane and 60 per cent nonane. The temperature is everywhere constant at the start of the simulation. Values for pressures, domain length and temperature are chosen to match with parameters in Trangenstein & Bell (1989).

Figure 4 shows the solution at *t*=5×10^{7} s (579 days) using the MUSCL scheme with 1000 cells and a CFL number of 0.9. Wave structures very similar and with comparable resolution to those in Trangenstein & Bell (1989) can be observed. The injected vapour mixes with the reservoir liquid to form a two-phase mixture, with about half of the reservoir fluid being swept out by the leading front. The temperature remains within 5 K of the initial uniform temperature. There is a small overshoot in the *n*_{1} component at the left-going shock, where the flow transitions from gaseous to two-phase, and similar small overshoots at the phase boundaries on the temperature plot. These are also present if the case is solved with a first-order HLL method, which suggests that they are due to the phase equilibrium computation. The left-going and right-going shocks are captured within two to three cells, but the two-phase region (which includes the contact wave) requires more cells to resolve. In particular, in the region near *x*=1000 cm the Buckley–Leverett eigenvalue (labelled as *λ*_{1} in figure 4) crosses the other two eigenvalues. This results in loss of strict hyperbolicity, making this region particularly difficult to resolve. We next consider a thermal injection test case. This test models the injection of a hot gas into a cold liquid reservoir. This test is identical to test 1, except now the injected mixture is at a temperature of 800 K, and the rock has porosity 0.4 and heat capacity 1.5 J cm^{−3} K^{−1}. This test was solved in our previous paper (van Odyck *et al.* 2009) using a first-order upwind method.

The solution at *t*=3×10^{7} s (347 days) is shown in figure 5. The MUSCL scheme was used with 400 cells and a CFL number of 0.9. Results are comparable with those in van Odyck *et al.* (2009). The contact wave in the two-phase region (at about *x*=2000 cm) is captured with noticeably better resolution by the method used here. Finally, we consider a one-dimensional gravity inversion test. A cold liquid containing 20 per cent methane, 20 per cent butane, 60 per cent nonane is placed on top of a hot vapour containing 95 per cent methane, 4.9 per cent butane and 0.1 per cent nonane. The left and right boundaries are no-flow boundaries; that is, mass and enthalpy flux is set to zero at these boundaries and the pressure gradient across the boundary is specified in such a way that the total velocity *v*_{T} is zero. Since we neglect capillary pressure, it is not possible to specify that both phase velocities are zero; however, since the phase velocities do not appear in the flux function this poses no numerical difficulty.

The system is subjected to a gravity of *g*_{x}=−980 cm s^{−2}. Porosity is uniformly 1. The computed solution at time *t*=4×10^{4} s using the MUSCL scheme with 12 000 cells is shown in figure 6. Numerical convergence to this solution is also shown for a first-order HLL scheme and the MUSCL scheme. The solution displays a complex wave structure in the mixing zone between the cold liquid and the hot vapour. At the liquid phase boundaries, we see a shock, which is directly connected to a rarefaction wave: this is a so-called compound wave. The Buckley–Leverett wave displays a spike just after the shock. This spike is a common feature for a locally linearly degenerate wave, and occurs because the smeared shock passes through a region of phase-space in which the wave speed is higher than the shock speed (Bell *et al.* 1989). In the mixture region, the Buckley–Leverett wave crosses the other waves, giving rise to eigenvector deficiencies. Nevertheless, the scheme handles these challenging regions well.

### (c) Work comparison

To test the efficiency of the method, we compared the computation expense required for the first-order and the MUSCL advection schemes. Figure 7 shows computation time against solution error for two of the one-dimensional test cases previously considered: the smooth test and the thermal gas injection test. Grid sizes ranged from 50 to 1600 cells. Curvature of the plots for low computational time can be attributed to static initialization costs. The figure demonstrates the increase in efficiency given by the MUSCL scheme, especially as applied to the smooth test. For the thermal injection test, the MUSCL scheme is only first-order owing to the presence of multiple discontinuities. However, the scheme still demonstrates good efficiency, and for small errors, it is about three to five times faster than the first-order scheme. This is mainly due to the improved discontinuity resolution.

### (d) Two-dimensional test case

In this section, the gravity inversion test (test 3) is repeated in two dimensions, using a medium with porosity 0.4 and heat capacity 1.5 J cm^{−3} K^{−1}. The domain is a square of side length 100 cm with no-flow boundary conditions at each edge. The interface between the two fluids is located at *y*=50 cm and is perturbed at *t*=0 with a sine wave of amplitude 3 cm. The system is initially at a constant pressure of 8 MPa and is allowed to evolve under a gravity of 980 cm s^{−2} in the negative *y*-direction. The MUSCL scheme was used with a grid of 400^{2} cells and a CFL number of 0.9.

Figure 8 shows the computed solution at *t*=10000 s≈2.8 h. The initial interface perturbation has resulted in the formation of long multiphase fingers which have themselves split into sub-fingers. Because the flow is unstable, the exact nature of these fingers depends on the numerics used, including how the initial conditions are implemented (here we used volume averaging). However, because the fingers are two-phase with only a small vapour component, their mobility ratio with the liquid is only about 1.4. For this reason, the flow is much less unstable than might be expected from a simple comparison of the mobilities of the liquid and vapour phases. For larger mobility ratios (above about 10), a simple dimensionally split scheme will generally suffer from strong grid orientation effects, and further techniques would be needed as discussed at the end of §7.

Finger detail varies if different methods, resolutions or CFL numbers are used. This is demonstrated in figure S1, electronic supplementary material. Note that only those simulations with the MUSCL scheme at high resolution demonstrate splitting of sub-fingers. In this case, lowering the resolution, using a low CFL number or using a first-order scheme tends, to increase numerical diffusion and eliminate splitting of sub-fingers.

Figure 8*f* shows the volume discrepancy at each point in the domain. This discrepancy is generally largest at phase boundaries. For this simulation, the discrepancy is *O*(10^{−2}) (maximum 0.028), which is acceptable for a two-dimensional simulation involving complex phase behaviour. See for example Dicks (1993, p. 95)], where this order of discrepancy is found typical of two-dimensional black-oil computations.

Finally, we note that mass and energy conservation are essential for thermal problems of this type. As the MUSCL scheme uses conservative variables for the advection step, we expect stable numerical conservation with floating-point imprecision the only source of conservation error. Figure S2, electronic supplementary material demonstrates this for the two-dimensional test using the MUSCL scheme.

## 9. Summary

In this paper, we derived a primitive-variable-based MUSCL-type numerical scheme for equations describing thermal multicomponent multiphase flow in porous media. This extends the work in van Odyck *et al.* (2009).

The basic approach was first tested on a model problem (the polymer flood model) in order to demonstrate its applicability to a system exhibiting local linear degeneracies and eigenvector deficiencies. The scheme was then derived for the porous media flow equations and applied to four test cases. A convergence study indicated that the scheme displays second-order convergence for moderate-sized grids (up to 1600 cells in one dimension), even though the overall method necessarily includes a first-order temporal error due to the first-order discretization of the pore volume constraint. Compared with a first-order scheme, we see an improvement in the value of the error norms of one to two orders of magnitude for smooth problems. For discontinuous problems, we see better resolution of flow discontinuities and thus greater computational efficiency (up to five times in our tests).

- Received March 9, 2012.
- Accepted June 6, 2012.

- This journal is © 2012 The Royal Society