Simulation of Three-Component Two-Phase Flow in Porous Media Using Method of Lines

Numerical simulation of compositional flow problems commonly involves the use of 1st- or 2nd-order Euler time stepping. Method of lines (MOL), using highly accurate and efficient ODE solvers, is an alternative technique which, although frequently applied to the solution of two-phase, two-component flow problems, has generally been overlooked for problems concerning more than two components. This article presents the development of a numerical simulator for 1D, compressible, two-phase, three-component, radially symmetric flow using the method of lines (MOL) and a 3rd-order accurate spatial discretization using a weighted essentially non-oscillatory (WENO) scheme. The MOL implementation enables application of the MATLAB ODE solver, ODE15s, for time integration. Simulation examples are presented in the context of CO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {CO}_2$$\end{document} injection into a reservoir containing a mixture of CH4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {CH}_4$$\end{document} and H2O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {H}_2\hbox {O}$$\end{document}. Following an assumption of constant equilibrium ratios for CO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {CO}_2$$\end{document} and CH4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {CH}_4$$\end{document}, a ternary flash calculator is developed providing closed-form relationships for exact interpolation between equations of state for CO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {CO}_2$$\end{document}–H2O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {H}_2\hbox {O}$$\end{document} and CH4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {CH}_4$$\end{document}–H2O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {H}_2\hbox {O}$$\end{document} binary mixtures. The numerical code is successfully tested and verified for a range of scenarios by comparison with an existing analytical solution.

an infinitesimal volume. These phases typically include aqueous liquid, non-aqueous liquid, gas and solid. Each phase can represent a mixture of multiple components (e.g., H 2 O, CO 2 and CH 4 ). The volume fraction of the jth phase is mathematically represented as a volume fraction of the pore space, often referred to as the saturation, S j [-]. Multiple components are then considered as mass fractions of the various phases. For example, X i j would denote the mass fraction of the ith component in the jth phase.
Mathematical simulation of such problems involves solving a coupled set of conservation equations for each component (Orr 2007). Equations of state (EOS) are required to understand how the components partition into different phases. This is particularly important because the mobilities of the components are strongly controlled by their distribution between the present fluid phases. Firstly, component mass fractions, X i j , affect the viscosity of the phases (Centeno et al. 2011). Secondly, the permeability available to each phase is a nonlinear function of the phase saturations, S j (Mathias et al. 2013).
In the absence of diffusion and capillary pressure, this problem is governed by a set of coupled hyperbolic and parabolic transport equations. Hyperbolic equations frequently give rise to the formation of shocks, leading to difficulties with regard to obtaining accurate solutions. Problems associated with one-dimensional transport of incompressible fluids, in the absence of capillary pressure and under isothermal conditions, can be solved exactly using the method of characteristics (Orr 2007). However, even under these restricting conditions, great care must be taken when considering nonzero initial conditions and non-unity boundary conditions (e.g., see Section 4.3 of Orr (2007)).
Alternative techniques involve the application of approximate methods. The spatial dimension is typically treated using conservative methods such as finite volume (Chen et al. 2006, p. 128). Alternatively, one can consider the use of finite elements (Chen et al. 2006, p. 94) or pseudospectral methods (van Reeuwijk et al. 2009). Such spatial schemes give rise to either stability problems or numerical diffusion due to truncation terms associated with the Taylor's expansion, the latter of which can be reduced using flux limiters or their variants (e.g., Mallison et al. 2005).
Handling of the temporal term, which is critical to resolving the nonlinear nature of the problem, generally revolves around the choice of explicit or implicit treatment. Fully explicit treatment, although easier to implement, can run into severe time-step limitations due to the well-known Courant-Friedrichs-Lewy (CFL) condition. Fully implicit treatment leads to an unconditionally stable solution (as far as time stepping is concerned), but leads to additional numerical diffusion. Furthermore, implementation of the solution is significantly more challenging.
Popular approaches for solving MCMP problems in this context are the so-called semiimplicit methods, the most common variant of which is referred to as implicit pressure explicit saturation (ImPES) (Chen et al. 2006) or implicit pressure explicit mass (ImPEM). In ImPES, the governing equations are rearranged to identify a transport equation of hyperbolic (or nearly hyperbolic) nature and a pressure equation (of parabolic or elliptic character). The pressure equation is solved implicitly which allows for larger time steps. The transport equation is solved explicitly, allowing easier implementation and reduced computational memory requirements, hence the semi-implicitness. Both the implicit time stepping and explicit time stepping typically employ simple first-order schemes.
Multi-step multi-order time integration algorithms (Shampine and Reichelt 1997) represent an alternative method, which treats the temporal term in a more accurate fashion. These techniques maintain a specific time integration error while maximizing the time-step size. Moreover, due to the wide availability (e.g., MATLAB or FORTRAN with NAG) of high-quality solvers and simplicity of implementation, there is no need to redevelop the sophisticated solution algorithms. Rather, the so-called method of lines (MOL) approach can be taken. In this case, the partial differential equations (PDE) are discretized in space to form a set of coupled ordinary differential equations (Wouwer et al. 2005). These can then be solved simultaneously using any ordinary differential equation (ODE) solver of choice.
For two-phase immiscible flow, where one of the phases is treated as inviscid, the MCMP problem reduces to a single PDE often referred to as Richards' equation (RE). This equation is commonly solved to better understand hydrological problems associated with unsaturated soils. Numerical studies by Farthing et al. (2003) and Kees and Miller (2002) have shown that applying MOL with higher-order time integration to the solution of RE leads to both improved accuracy and computational efficiency. Indeed, there are many recent articles (Mathias et al. 2006(Mathias et al. , 2008bIreson et al. 2009) reporting MOL solutions of RE using the MATLAB ODE solver, ODE15s, which is particularly suitable for stiff systems of ODEs (Shampine and Reichelt 1997). ODE15s has also been found to provide useful solutions to non-Darcian flow problems (Mathias et al. 2008a;Mathias and Wen 2015;Wen et al. 2009) and two-phase immiscible flow problems ).
Often when dealing with partial differential equations, it is useful to distinguish between dependent and independent variables (Stroud and Booth 2007, p. 122). In this case, time and space are independent variables. All other variables are dependent variables.
Application of MATLAB ODE solvers to multi-component partially miscible problems has proven more challenging. Consider N c components residing in N p phases. The problem will be defined by N c mass conservation equations. However, considering the various values of S j and X i j , it can be understood that there will be at least (N c + 1)N p dependent variables. It is therefore necessary to choose N c dependent variables to solve for. Special care should be taken to ensure that the dependent variables selected are persistent (Amaziane et al. 2012;Bourgeat et al. 2013). This selected set of dependent variables is hereafter referred to as the primary dependent variables (PDV).
When using an ODE solver, the user must construct an ODE function. Within this function, a scalar value of time is provided as an input along with an associated vector of the PDVs. The user must define the ODE function such that it calculates the derivatives of the PDVs with respect to time, which generally involves using a combination of chain rule and product rule differentiation. This results in a need to evaluate the partial derivatives of the bulk fluid mass per unit pore space, F [ML −3 ], with respect to each component mass fraction, z i [-] (these terms will be mathematically defined later in the article). For conventional first-order time stepping, it is arguably acceptable to evaluate these derivatives using first or second-order finite differencing. However, given the high accuracy associated with the use of MATLAB's ODE solvers, it is pertinent to obtain these derivatives in exact form wherever possible.
There are many detailed works concerning applications of MOL for immiscible two-phase flow and two-component two-phase flow problems (e.g., Amaziane et al. 2012;Vohralik and Wheeler 2013;Bourgeat et al. 2013;Mathias et al. 2014). Mallison et al. (2005) present a numerical simulation of an MCMP problem using MOL in conjunction with a 3rd-and 4th-order Runge-Kutta time integration method. However, Mallison et al. (2005) provide no discussion concerning the casting of equations in terms of PDVs. Indeed, little information is available as to how to obtain exact equations to describe the necessary partial derivatives, ∂ F/∂z i , needed to solve MCMP problems for situations concerning more than two components. In this article, we focus on obtaining such expressions for three-component two-phase problems. These are implemented within a radial flow simulator using MATLAB. Comparisons are then made, in the context of enhanced gas recovery by CO 2 injection, with an associated analytical solution, previously presented by Hosseini et al. (2012).

Governing Equations
Consider three components: CO 2 , CH 4 and H 2 O, denoted hereafter as i = 1, 2 and 3, respectively. The three components can partition into a gas phase and an aqueous liquid phases, denoted hereafter as j = 1 and 2, respectively. A horizontally orientated, homogeneous and isotropic cylindrical reservoir of radius, r e [L], and formation thickness, H [L], is invoked. The reservoir is initially filled with a mixture of CH 4 and H 2 O· CO 2 is injected into the center of the reservoir via a fully completed vertically orientated injection well of radius, r w [L]. Fluid flow is assumed to be one-dimensional such that the problem reduces to the following set of one-dimensional radially symmetric conservation equations: and q j [LT −1 ] are the density, saturation (a volume fraction of the pore space) and volumetric fluid flux of phase j, respectively, and X i j [-] is the mass fraction of component i in phase j. Note that i=1 Nc X i j = 1 and j=1 N p S j = 1. The volumetric fluxes are calculated using Darcy's law: where k [L 2 ] is the reservoir permeability and k r j [-], μ j [ML −1 T −1 ] and P j [ML −1 T −2 ] are the relative permeability, dynamic viscosity and pressure of the jth phase, respectively. For two-phase flow, without loss of generality, the relative permeability functions are assumed to take the form of power laws (Mathias et al. 2013): where S jc [-], k r j 0 [-] and n j [-] are the critical saturation, end-point relative permeability and power-law exponent for phase j, respectively. Furthermore, it is assumed that the following density mixing rule applies (Orr 2007) where ρ i j [ML −3 ] is the density of the ith component in the jth phase.

Recasting in Terms of Primary Dependent Variables
An appropriate choice of primary dependent variables (PDVs) to solve for is the global fluid pressure, P [ML −1 T −2 ], defined in this case by Chen et al. (2006), p. 342 and the bulk mass fraction of each component, z i [-], defined by where F [ML −3 ] is the bulk fluid mass per unit volume of rock Note that N c i=1 z i = 1. In some previous studies, the mass of each component per unit volume of rock, G i , has also proven effective as PDVs in this context (Amaziane et al. 2012;Bourgeat et al. 2013). However, an advantage of using P and z i (for i = 1, 2, . . . N c − 1) as PDVs (as opposed to say G i ) is that z i are independent of P. For a given volume of fluid mixture, the mass fractions of each component, z i , will not change with pressure. However, the associated mass of each component per volume of rock, G i , may change with pressure, depending on how the individual component mass densities, ρ i j , vary with pressure. Furthermore, z i are the variables used in the ternary diagram (discussed later in the article), which determine the equilibrium properties of the three-component fluid mixture.
Differentiating Eq. (9) with respect to time leads to where, from Eqs. (1) and (10) Application of the chain rule to Eq. (10) and rearranging leads to The main focus of this article is the derivation and application of exact formulae for the relationships defining ∂ F/∂ P and ∂ F/∂z i for i = 1, 2, . . . , N c − 1. Note that Eqs. (11)-(13) are important because they directly relate the time derivatives of the PDVs to the original mass conservation statements.

Differentiating the F Function
Considering the identity in Eq. (10), the total derivative of F can be written as When there are only two phases, S 2 = 1 − S 1 . From Eqs. (2), (9) and (10), it can then be understood that which on differentiation leads to Invoking Eq. (7), it can also be shown that It is now assumed that component densities are unaffected by composition such that ∂ρ i j /∂z i = 0. Additionally noting that the remaining challenge is to define the d X i j terms.
The terms X i j can be further defined by: where x i j [-] are the equilibrium mass fractions of the ith component in the jth phase within the two-phase region. Then d X i j is given by: It is important to note that for two-component two-phase problems, x i j only varies with pressure. However, for three-component two-phase systems, the problem is much more complicated because the equilibrium constants, x i j , are no longer constant with composition (i.e., z 1 and z 2 ).
EOSs for the binary mixtures of CO 2 -H 2 O and CH 4 -H 2 O are well characterized by relatively simple sets of equations. In this work, the binary mixture properties of CO 2 -H 2 O and CH 4 -H 2 O were obtained using expressions provided by Spycher et al. (2003) and Duan  (2006), respectively. These can be joined together to form a ternary system by assuming that the equilibrium ratios, K i [-], defined as for components 1 and 2 (but not 3) do not vary with composition. When K 1 and K 2 are constant, the two-phase region, for a given temperature and pressure, is defined by two straight lines on a ternary diagram (see Fig. 1). Consequently, following the work of Juanes (2008), it can be shown that values of x i j can be related back to their values obtained from the binary mixtures, a i j , by the set of linear equations: where a 1 j = x 1 j when z 2 = 0, a 2 j = x 2 j when z 1 = 0 and A [-] is a weighting parameter that linearly interpolates between the bounding tie-lines that coincide with the z 2 and z 3 axes of the ternary diagram. Writing out Eq. (23) for components 1 and 2 and eliminating A lead to: Similarly it can be said that, on a given tie-line in the two-phase region, where B [-] is a weighting parameter that linearly interpolates between the bounding lines that define the two-phase region of the ternary diagram. Writing Eq. (25) out for components 1 and 2 and then eliminating B lead to Note that the equilibrium ratio for H 2 O (component 3) is not necessarily constant and is defined by the relationship Using Eq. (24) to eliminate the x i2 and x 2 j terms yields the quadratic equation where T 0 = (a 21 − a 22 )a 11 z 1 (29) T 2 = (a 12 a 21 − a 11 a 22 )/a 11 (31) which has the solutions Equation (32) gives an explicit expression for x 11 with respect to z 1 and z 2 . Furthermore, from Eq. (28) we have the total derivative of x 11 : where  (21) and (20).

Additional Considerations Required for Accommodating Capillary Pressure
When capillary pressures can be assumed negligible, the information in this section can be ignored. However, for non-negligible capillary pressure, the fluid properties x i j and ρ i j should be calculated from the phase pressure, P j , as opposed to the global pressure, P. An equation of state can provide derivatives of these variables with respect to P j . But to obtain derivatives with respect to P, the following transformations must be applied: For two-phase flow systems, Eq. (8) reduces to where P c = P 1 − P 2 is the capillary pressure.
Noting that P c is generally expressed uniquely as a function of S 1 (e.g., van Genuchten 1980) and that d S 1 = −d S 2 , differentiating Eq. (39) with respect to P leads to Recalling Eq. (16) and that z i are independent of P, it can be said that where from which we obtain

Numerical Solution
A weighted essentially non-oscillatory (WENO) scheme (Shu 2009) was used to approximate the flux, H i , in Eq. (1). For more details on the WENO method and its applications, the reader is referred to Coralic and Colonius (2014), Noelle et al. (2007) and Zhang and Shu (2012) and references therein. The WENO scheme implemented here uses a 2nd-order stencil to approximate the flux at the mid-points (k + 1/2 and k − 1/2). A central difference approximation is then used to evaluate the derivative of the flux at grid point k, which is 3rd-order accurate in space (Shu 2009): The pressure derivative contained within the (H i ) k+ 1 2 term was approximated using the following central difference expression: The Jacobian matrix is a matrix describing the dependencies of each equation in the system of equations on every other equation. In this problem, we are solving for three PDVs, which, by virtue of the discretised conservation equations, are mutually dependent on one another. Conventional second-order finite difference approximations of second-order derivatives at a grid point, k, will depend on three points, namely k − 1, k and k + 1. However, when using the WENO scheme described above, additional information from k − 2 and k + 2 is required, leading to a penta-diagonal structure for the Jacobian matrix. For systems of equations with more than one PDV, a so-called block-diagonal structure is formed.
In MATLAB, it is possible to provide the ODE solver with the Jacobian pattern as a sparse matrix of zeros and ones to indicate where the Jacobian matrix is nonzero. This tells the solver to only evaluate the sparse system and not the full matrix. In stiff problems with non-constant  Fig. 2 Block penta-diagonal Jacobian pattern for the sparse system for n = 10 grid points and 3 PDVs (nz number of nonzero elements) Jacobian matrices, specifying the Jacobian pattern a priori can lead to significant savings in computation time. A pictorial representation of the Jacobian pattern for the problem discussed in this article is presented in Fig. 2.
Spatial grids were locally refined near the well bore and gradually coarsened away from the injection well to ensure simulation stability. Simulation time steps are refined automatically using MATLAB's ODE15s adaptive time stepping scheme. ODE15s is particularly suitable for stiff problems where the governing equations include a combination of terms that lead to rapid variation in the solution and terms that vary slowly. For the problem described in this article, compositions propagate very slowly thorough the reservoir compared to the pressure waves, which move much more rapidly.
ODE15s is a multi-order multi-step solver. The scalar relative error tolerance and the absolute error tolerance are set to MATLAs default values, 10 −3 and 10 −6 , respectively. For more information about how the solver is implemented, readers are referred to Shampine and Reichelt (1997) and Shampine and Thompson (2001).

Model Verification
Results from the numerical solution (implemented in MATLAB) were compared against results from the analytical solution by Hosseini et al. (2012). Both models assume a fully completed well at the center of a 1D radially symmetric flow field of radial extent, r e . Capillary pressure, molecular diffusion and gravity effects were neglected. Initial and boundary conditions were applied as follows: where S 10 [-] is the initial gas saturation and the other parameters are as specified in Tables 1 and 2. Figure 3 shows gas saturation profiles for CO 2 injection into deep (Fig. 3a, b) and shallow (Fig. 3c, d) reservoirs with 10 % residual CH 4 , during a simulation period of 1000 days. The solid lines and dashed lines represent results from the analytical solution and numerical solution, respectively.
The high gas saturation around the injection well, often referred to as a dry-out zone (Mathias et al. 2011), is due to the vaporization of the residual water saturation by the injected CO 2 . A CH 4 bank with about 22 % gas saturation in front of the CO 2 plume develops. The length of the CH 4 bank increases with time. Correspondence between the analytical solution and the numerical solution is very good for both gas phase saturation and pressure buildup. The numerical solution is seen to accurately locate the associated shock fronts while considering the partial miscibility of both CO 2 and CH 4 in H 2 O. Fig. 3 a Comparison of gas saturation profiles between the analytical and numerical simulation of CO 2 injection into a deep reservoir with initial gas saturation, S 10 = 0.1, after 10, 100 and 1000 days. b Corresponding pressure profiles for deep reservoir. c Gas saturation profiles for a shallow reservoir. d Corresponding pressure profiles for a shallow reservoir. A total of 300 nodes were used for the numerical simulation For the shallow reservoir, just ahead of the CH 4 bank, it can be seen that the numerical solution predicts a slightly lower gas saturation as compared to the analytical solution. This is due to the fact that the analytical solution assumes constant fluid properties and hence is not capturing volume change effects due to pressure change. This discrepancy is not noticeable for the deep reservoir because the gas compressibility is a lot lower at 34 MPa. The possibility that the above-mentioned discrepancy was due to insufficient mesh refinement was investigated by comparing simulation results using a 600 point grid and a 900 point grid. The comparison study confirmed that the discrepancy was not due to numerical error.
Similar simulations but with different initial gas saturations are compared in Fig. 4. It is found that the extent of the dry-out region is insensitive to the initial gas saturation. The extent of the dry-out region is smaller for the shallow reservoir, and the volume of the gas plume is larger. The reduced dry-out region, in this case, is due to the reduced evaporation that occurs at cooler temperatures. The increased gas volume is due to the reduced gas density that occurs at lower pressures. Again it can be seen that the numerical solution is able to accurately predict the analytical results of Hosseini et al. (2012).

Numerical Solution Performance
The performance of the MOL numerical solution was explored further through a grid convergence study. The shallow reservoir scenario, depicted in Fig. 3c, d, was repeated using different numbers of grid cells. Numerical solution performance was quantified by calculating the mean absolute normalized error (MANE) between results (from each grid cell) from the numerical solution and those from the analytical solution for a given time. Figure 5 shows a plot of MANE for gas saturation and pressure for the different times previously presented in Fig. 3c, d. Grid convergence can be seen to have been achieved at around 300 grid cells.
The converged MANE for pressure is around 0.02 %. However, the converged MANE for gas saturation is quite high at around 0.25 %. This is due to conceptual differences between the numerical solution and the analytical solution. Recall that the analytical solution assumes constant fluid properties, whereas the numerical solution is allowing for variations of fluid properties with pressure and composition, as discussed above. Hosseini et al. (2012) were able to achieve a very similar level of accuracy for a very similar set of simulations using the commercial reservoir simulator, CMG GEM (Computer Modeling Group Ltd. 2015) (see Fig.  5 of Hosseini et al. (2012)). GEM uses an adaptive-implicit solver as described by Collins et al. (1992).
In terms of computation time, the 300 grid cell numerical solution using our new MOL approach was able to simulate 1000 days of gas injection in around 5 min using an Intel Xeon CPU E5-2630 2.30 GHz (2 processors). For comparison, the GEM simulations undertaken by Hosseini et al. (2012) took around 2 h on an Intel Xeon CPU E5-2687 3.10 GHz (2 processors) (Hosseini, 2015, Personal Communications). Clearly, our new MOL solution has the potential to offer significant computation time savings in this context.

Effect of Initial Gas Saturation
As was shown in the previous section, numerical simulations of CO 2 injection into a reservoir containing CH 4 predict the accumulation of a CH 4 bank at the head of the CO 2 plume (Oldenburg and Doughty 2011;Battistelli and Marcolini 2009;Taggart 2010).
The system in discussion can be differentiated into three regions (Hosseini et al. 2012). These regions, starting from the injection point and moving outward, are: 1. a single-phase, dry-out region around the well bore filled with pure CO 2 . 2. a two-phase, two-component system containing CO 2 and H 2 O. 3. a two-phase, two-component system containing CH 4 and H 2 O.
Within the two-phase mixture, each phase propagates at a rate according to its mobility. The mobility of each phase varies from one region to another due to associated compositional changes. As a consequence, a trailing shock forms at the contact between regions (1) and (2) and a leading shock forms at the contact between regions (2) and (3).
The development of the CH 4 bank ahead of the CO 2 has been explained as follows (Taggart 2010;Oldenburg and Doughty 2011;Hosseini et al. 2012). As CO 2 is injected, it partitions into the gas phase and the aqueous phase. The initially dissolved CH 4 exsolves immediately and is then pushed ahead of the growing CO 2 plume leading to the development of a CH 4 bank (Oldenburg and Doughty 2011). Mathematically, the system is constrained to constantly enter and leave the two-phase region along the tie-lines representing the injection and initial compositions; therefore, the leading CH 4 bank is free from injected gas, CO 2 (Taggart 2010).
Intuitively, it is expected that the amount of CH 4 initially present should affect the methane bank saturation; the more the initial CH 4 saturation, the higher the bank saturation. However, numerical simulation of CO 2 injection for different initial gas saturations (everything else being the same), shows that the bank saturation is independent of the initial CH 4 saturation (see Figs. 4,7). In fact, Hosseini et al. (2012) showed mathematically that the CH 4 bank saturation is independent of the initial gas saturation.
This can be further explained using the principles of fractional flow theory (Pope 1980;Orr 2007). Because of the differences in phase viscosities in the two-phase region (i.e., between mixtures of CO 2 -H 2 O and CH 4 -H 2 O), flow occurs on different fractional flow curves in the two-phase region (Fig. 7). Figure 7a, c and e shows the fractional flow curves (plots of H i /ρ i1 against G i /ρ i1 ) for CO 2 and CH 4 along with the locations of the shock fronts for different initial gas saturations.
The partial derivative ∂ H i /∂G i represents the wave velocity of the system. The wave velocities of the shock fronts are found from the gradients of straight lines that link the two conditions on either side of the shock. Fractional flow theory dictates that valid solutions should satisfy both the velocity constraint and the so-called entropy constraint (Orr 2007). The velocity constraint implies that wave velocity should always decrease with increasing distance from the injection boundary. The entropy constraint implies that the shock wave velocity should be equal to the gradient of the fractional flow curve immediately upstream of the shock.
Due to the zero initial condition for G 1 , the only valid path on the CO 2 fractional flow curve is a tangent (i.e., (0,0) to L). On the other hand, velocities must be equal at the contact between a pair of different fluids (Pope 1980) (i.e., points G and L in Fig. 6). This means that the gas bank saturation (point G) is dictated by the intersection of the tangent to the CO 2 curve with the CH 4 curve. Figure 7b, d and f shows the corresponding saturation profiles for different initial gas saturations. The level of saturation at point G is always determined by the tangent from (0,0) to point L. Physically, this implies that the bank saturation is only dependent on how fast the injected gas propagates. The solid and dashed lines in Fig. 7b, d and f are from the analytical solution and numerical solution, respectively. There is an excellent correspondence between the two. The analytical solution was developed on the basis of the fractional theory described above. The numerical solution therefore further confirms the finding of Hosseini et al. (2012) that the CH 4 bank saturation is independent of the initial gas saturation.
However, if the initial gas saturation is sufficiently high it may not be possible to build a solution that travels from the initial condition (point I in Fig 6) to the intersection point of the tangent line with the CH 4 fractional flow curve (point G in Fig 6). As a consequence, a CH 4 bank will no longer be present. Therefore, it should be noted that the existence of the CH 4 bank is dependent on the initial gas saturation (LaForce and Johns 2010). Nevertheless, providing the CH 4 bank exists, the CH 4 bank saturation is independent of the initial gas saturation.

Summary and Conclusions
A numerical simulator was developed for 1D, compressible, two-phase, three-component, radially symmetric flow using the method of lines (MOL) and a 3rd-order accurate spatial discretization using weighted essentially non-oscillatory (WENO) scheme. The MOL implementation enabled application of the state-of-the-art MATLAB ODE solver, ODE15s, for time integration. Pressure and overall component mass fraction (z i ) were selected as the primary dependent variables (PDV). The sparsity of the system of equations was taken advantage of by provision of a penta-diagonal Jacobian pattern for the ODE solver. Simulation examples were developed in the context of CO 2 into a reservoir containing a mixture of CH 4 and H 2 O.
Following an assumption of constant equilibrium ratios (K i ) for CO 2 and CH 4 , it was possible to derive a ternary flash calculator by deriving closed-form relationships for exact Illustration of CH 4 bank saturation independence of initial CH 4 mass fraction for the deep reservoir scenario after 1000 days of injection (with parameters as set in Table 1): a and b S 1 I < S 1 G , c and d S 1 I = S 1 G , e and f S 1 I > S 1 G interpolation between equations of state for CO 2 -H 2 O and CH 4 -H 2 O binary mixtures. This helped ensure improved stability and mass conservation and led to reduced computational requirements associated with multi-dimensional lookup tables. The numerical code was successfully tested and verified for a range of scenarios by comparison with an existing analytical solution.