A Moncorgé

H A Tchelepi

P Jenny

Sequential Fully Implicit Formulation for Compositional Simulation using Natural Variables

Keywords: nonlinear dynamics, numerical flow simulation, sequential implicit, operator splitting, coupled flow and transport, multiscale methods, flow in porous media, compositional formulation, multiphase flow, multi-component transport, compositional reservoir simulation

The Sequential Fully Implicit (SFI) method was proposed (Jenny et al., JCP 2006), in the context of a Multiscale Finite Volume (MSFV) formulation, to simulate coupled immiscible multiphase fluid flow in porous media. Later, Lee et al. (Comp. Geosci. 2008) extended the SFI formulation to the black-oil model, whereby the gas component is allowed to dissolve in the oil phase. Most recently, the SFI approach was extended to fully compositional isothermal displacements by Moncorgé et al., (JCP 2017). SFI schemes solve the fully coupled system in two steps: (1) Construct and solve the pressure equation (flow problem). (2) Solve the coupled species transport equations for the phase saturations and phase compositions. In SFI, each outer iteration involves this two-step sequence. Experience indicates that complex interphase mass transfer behaviors often lead to large numbers of SFI outer iterations compared with the Fully Implicit (FI) method. Here, we demonstrate that the convergence difficulties are directly related to the treatment of the coupling between the flow and transport problems, and we propose a new SFI variant based on a nonlinear overall-volume balance equation. The first step consists of forming and solving a nonlinear pressure equation, which is a weighted sum of all the component mass conservation equations. A Newton-based scheme is used to iterate out all the pressure dependent nonlinearities in both the accumulation and flux terms of the overall-volume balance equation. The resulting pressure field is used to compute the Darcy phase velocities and the total-velocity. The second step of the new SFI scheme entails introducing the overall-mass density as a degree-of-

freedom, and solving the full set of component conservation equations cast in the natural-variables form (i.e., saturations and phase compositions). During the second step, the pressure and the total-velocity fields are fixed. The SFI scheme with a nonlinear pressure extends the SFI approach of Jenny et al. (JCP 2006) to multi-component compositional processes with interphase mass transfer. The proposed compositional SFI approach employs an overall balance for the pressure equation; however, unlike existing volume-balance Sequential Implicit (SI) schemes (Acs et 2014), which use overall compositions, this SFI formulation is well suited for the natural variables (saturations and phase compositions). We analyze the 'splitting errors' associated with the compositional SFI scheme, and we show how to control these errors in order to converge to the same solution as the Fully Implicit (FI) method. We then demonstrate that the compositional SFI has convergence properties that are very comparable to those of the FI approach. This robust sequential-implicit solution scheme allows for designing numerical methods and linear solvers that are optimized for the sub-problems of flow and transport. The SFI scheme with a nonlinear pressure formulation is well suited for multiscale formulations, and it promises to replace the widely used FI approach for compositional reservoir simulation.

Introduction

The Sequential Fully Implicit (SFI) method was first proposed to model multiphase fluid flow without mass exchange in the context of the Multiscale Finite Volume (MSFV) method [START_REF] Jenny | Adaptive fully implicit multi-scale finite-volume method for multi-phase flow and transport in heterogeneous porous media[END_REF]. The SFI approach was later extended to the black-oil formulation (three pseudo components with gas dissolved in the oil phase) [START_REF] Lee | Multiscale finite-volume formulation for multiphase flow in porous media: black oil formulation of compressible, three-phase flow with gravity[END_REF][START_REF] Møyner | A multiscale restriction-smoothed basis method for compressible black-oil models[END_REF], and recently for compositional multi-component flow by Hajibeygi & Tchelepi [11], Moncorgé et al. [START_REF] Moncorgé | New sequential scheme for mixed-implicit compositional flow simulation[END_REF][START_REF] Moncorgé | Modified sequential fully implicit scheme for compositional flow simulation[END_REF] and Møyner et al. [START_REF] Møyner | A multiscale restriction-smoothed basis method for compositional models[END_REF][START_REF] Møyner | A multiscale restriction-smoothed basis method for high contrast porous media represented on unstructured grids[END_REF]. In these previous works, the SFI algorithm consists of solving two different implicit systems in sequence: (1) a pressure equation, and (2) a saturation/composition system. During the solution of the saturation/composition system, the pressure is fixed. For n c components, a coupled set of n c -1 component conservation equations is solved in the saturation/composition step. Our objective is to design an SFI scheme with a nonlinear convergence rate that is comparable to that of the Fully Implicit (FI) formulation across the full range of compositional modeling of practical interest. This objective is motivated by several factors. First, multiscale formulations rely on sequential treatment of the flow (near-elliptic) and transport (hyperbolic) problems; thus, in order for multiscale methods to fully replace the existing single-level formulations, the SFI approach must have nonlinear convergence rates that are comparable, or even superior, to Fully Implicit (FI) solution schemes. Second, sequential treatment of flow and transport makes it possible to employ advanced discetization schemes and scalable solution methods to the different equations in the coupled system.

Moncorgé et al. [START_REF] Moncorgé | Modified sequential fully implicit scheme for compositional flow simulation[END_REF] enriched the pressure equation with Fully-Implicit (FI) coupling in regions that experience phase appearance to ensure convergence of the full system to any given tolerance. This approach is efficient as long as the fraction of cells that have both liquid and vapor phases is relatively small. However, as the number of cells that experience phase appearance/disappearance grows, this approach becomes expensive. In Lee at al. [START_REF] Lee | Multiscale finite-volume formulation for multiphase flow in porous media: black oil formulation of compressible, three-phase flow with gravity[END_REF] and Møyner et al. [START_REF] Møyner | A multiscale restriction-smoothed basis method for compressible black-oil models[END_REF][START_REF] Møyner | A multiscale restriction-smoothed basis method for compositional models[END_REF][START_REF] Møyner | A multiscale restriction-smoothed basis method for high contrast porous media represented on unstructured grids[END_REF], the conservation equation of one component is removed when solving the saturation/composition system. Hajibeygi & Tchelepi [START_REF] Hajibeygi | Compositional multiscale finitevolume formulation[END_REF] tried to remove the material-balance error by freezing the overall-mass density in the accumulation term and the total-mass flux between the computational cells. Their algorithm is convergent; however, in cases with strong coupling, the scheme needs large numbers of outer iterations. Acs et al. [START_REF] Acs | General purpose compositional model[END_REF], Watts [START_REF] Watts | A compositional formulation of the pressure and saturation equations[END_REF], Trangenstein & Bell [START_REF] Trangenstein | Mathematical structure of the blackoil model for petroleum reservoir simulation[END_REF][START_REF] Trangenstein | Mathematical structure of compositional reservoir simulation[END_REF], Dicks [START_REF] Dicks | Higher Order Godunov Black-Oil Simulations for Compressible Flow in Porous Media[END_REF], Pau et al. [START_REF] Pau | An adaptive mesh refinement algorithm for compressible twophase flow in porous media[END_REF], Faigle et al. [START_REF] Faigle | Multi-physics modeling of non-isothermal compositional flow on adaptive grids[END_REF] and Doster et al. [START_REF] Doster | A robust implicit pressure explicit mass method for multi-phase multi-component flow including capillary pressure and buoyancy[END_REF] developed Sequential Implicit (SI) formulations for compositional multi-component displacements. In these SI formulations, extensive degrees-of-freedom (variables), such as the overall number of moles, or mass, of a component divided by the pore volume, are used.

The developers of SI formulations noted that it is not possible to satisfy all the governing nonlinear equations exactly, and that some inconsistencies must be tolerated. The formulations of Acs However, these volume splitting errors are very local and have rarely large effect on the overall flow. On the contrary, the total-velocity splitting error, has never been documented before and has a much larger support. We show that we need to control this second splitting error to recover the fully-implicit solution. [START_REF] Coats | A note on IMPES and some IMPES-based simulation models[END_REF] proposed an IMPES (IMplicit Pressure Explicit Saturations) pressure equation built by algebraic manipulations of the usual material balance equations. Møyner et al. [START_REF] Møyner | A multiscale restriction-smoothed basis method for high contrast porous media represented on unstructured grids[END_REF][START_REF] Møyner | Comparison of splitting schemes for compositional multiscale solvers[END_REF] use this pressure for their sequential implicit method for compositional flow. The derivatives of the overall compressibility wrt pressure and the derivatives of the partial molar volumes wrt pressure are not accounted for in the linearized 'volume balance' equation and the derivatives of the partial molar volumes wrt pressure are not accounted for in the Coats pressure equation. Both equations experience nonlinear convergence difficulties, especially when the coupling between the multiphase flow and the multi-component transport is strong. Here, we propose a nonlinear 'volume balance' equation that accounts for all the necessary derivatives wrt pressure. For immiscible and incompressible multiphase flow, the nonlinear pressure equation that we derive here simplifies to the pressure equation of [START_REF] Young | Continuous compositional volume-balance equations[END_REF]. It is also identical to the pressure equation for immiscible multiphase flow of Jenny et al. [START_REF] Jenny | Adaptive fully implicit multi-scale finite-volume method for multi-phase flow and transport in heterogeneous porous media[END_REF] and to the pressure equation of Lee et al. [START_REF] Lee | Multiscale finite-volume formulation for multiphase flow in porous media: black oil formulation of compressible, three-phase flow with gravity[END_REF] and Møyner & Lie [START_REF] Møyner | A multiscale restriction-smoothed basis method for compressible black-oil models[END_REF] for the black-oil formulation.

This new SFI scheme with a nonlinear pressure equation can be seen as an adaptation of existing SI methods based on overall compositions to one suited for the natural-variables formulation. The splitting of the coupled system into two nonlinear problems, namely, a nonlinear pressure equation and a nonlinear multi-component transport problem, leads to volume-balance and total-velocity errors. We analyze these 'splitting errors' , and we show how to control them. We then demonstrate that the parabolic pressure equation and the hyperbolic transport systems are easier to solve separately using widely available linear solvers compared with the specialized solvers (e.g., CPR-AMG [START_REF] Cao | Parallel scalable unstructured CPR-type linear solver for reservoir simulation[END_REF]) needed to deal with the full system with mixed parabolichyperbolic behaviors.

The paper is organized as follows. First, we introduce a natural-variables formulation, which entails forming a nonlinear volume-balance equation followed by solving the full set of transport equations. Then, we elaborate the details of constructing the volume-balance equation. We also analyze the splitting errors and how to control them. We finally demonstrate the performance of the compositional SFI scheme for very challenging test cases.

Sequential Fully Implicit Scheme for Compositional Flow

In this section, we describe the compositional formulation based on naturalvariables for both Fully Implicit (FI) and Sequential Fully Implicit (SFI) nonlinear solution algorithms. We present a variant of the natural-variables formulation, whereby the system is augmented by an additional equation and an additional variable. Following that, the new SFI scheme for compositional flow is described.

Compositional Formulation Based on Extended Natural-Variables

The generalized compositional formulation, which accounts for multicomponent, three-phase flow with interphase mass transfer, is the target of this work. We assume that we have n c components, i.e., n h hydrocarbon components and the water component, and we employ the natural-variables formulation [START_REF] Coats | An equation of state compositional model[END_REF]. The overall density, ρ t , is used as an additional global variable. The number of unknowns in the case of three-phase flow is 2n h + 7. They are: gas-phase pressure, P g , oil-phase pressure, P o , water-phase pressure, P w , total mole density ρ t , the saturations of each phase, S g , S o , S w , the component mole (mass) fractions in the gas, y c , and the liquid, x c . To close the system, 2n h + 7 equations are required.

Mass conservation equations

The conservation of a hydrocarbon component, c, and of the water component, w, can be written as:

∂ ∂t [φ (y c ρ g S g + x c ρ o S o )] + ∇ • [y c ρ g u g + x c ρ o u o ] = y c ρ g q g + x c ρ o q o (1) and ∂ ∂t [φρ w S w ] + ∇ • [ρ w u w ] = ρ w q w , (2)
where φ is the porosity, ρ g , ρ o and ρ w are the fluid-phase mole densities, and q g , q o and q w source terms. The velocity of each phase p ∈ {g, o, w} is given by Darcy's law:

u p = - k rp µ p K ∇P p -ρ p g∇D = -λ p K ∇P p -ρ p g∇D , (3)
where K is the rock permeability, k rp the relative permeability of each phase, µ p the phase viscosities, λ p = kr p µp the phase mobilities, ρ p the fluid-phase mass densities, g the gravitational acceleration, D the depth. The total-velocity

u t = - p λ p K ∇P p -ρ p g∇D (4)
is the sum of the phase velocities. The phase velocities

u p = λ p λ t u t + q λ p λ q λ t K ∇P q -ρ q g∇D -∇P p -ρ p g∇D = λ p λ t u t + Ψ p (5)
can be expressed as a function of the total-velocity, where λ t = p λ p is the total-mobility and Ψ p = q λpλq λt K ∇P q -ρ q g∇D -∇P p -ρ p g∇D represents the capillary pressure and gravity contributions for each phase. Writing the phase velocities in terms of the total-velocity [START_REF] Coats | A note on IMPES and some IMPES-based simulation models[END_REF], and substituting into equations (1) and (2), rewriting the source terms, and adding the totalmole density ρ t as an additional variable in the accumulation term leads to:

∂ ∂t φ y c ρ g S g + x c ρ o S o p ρ p S p ρ t + ∇ • y c ρ g λ g λ t u t + Ψ g + x c ρ o λ o λ t u t + Ψ o = y c ρ g λ g λ t + x c ρ o λ o λ t q t (6) and ∂ ∂t φ ρ w S w p ρ p S p ρ t + ∇ • ρ w λ w λ t u t + Ψ w = ρ w λ w λ t q t , (7)
where q t = p q p . Equations ( 6) and ( 7) are the conservation equations in 'transport form'.

Constraint equations

The n h + 5 constraint equations involve variables in the control-volume (cell) under consideration. We have:

n h c=1 y c = 1, (8)
n h c=1 x c = 1 (9)
and

S g + S o + S w = 1. ( 10
)
We also have the capillary-pressure versus saturation relationships:

P c GO (S g ) = P g -P o (11) and P c OW (S w ) = P o -P w , (12)
where P c GO and P c OW are the capillary pressures between gas and oil phases and the oil and aqueous phases, respectively. The remaining constraints represent thermodynamic phase equilibrium for each hydrocarbon component (n h equations). These local equilibrium constraints are applied only when both hydrocarbon phases (oil and gas) are present in the control volume.

Existing formulations can be based on black-oil, or K-value correlations [START_REF] Crookston | A numerical simulation model for thermal recovery processes[END_REF], as well as the Peng-Robinson [START_REF] Peng | A new two-constant equation of state[END_REF], Redlich-Kwong [START_REF] Redlich | On the Thermodynamics of Solutions. V. An Equation of State. Fugacities of Gaseous Solutions[END_REF] and Soave-Redlich-Kwong [START_REF] Soave | Equilibrium constants from a modified Redlich-Kwong equation of state[END_REF] cubic equations of state (EOS) models. They can all be written as [START_REF] Martin | Cubic equations of state-which? Industrial & Engineering Chemistry Fundamentals[END_REF]: fc,g (P g , y 1 , ..., y n h ) = fc,o (P g , x 1 , ..., x n h ). ( 13)

Thermodynamic volume equation

The mole-fraction of each phase is denoted β p (P, z) with p ∈ {g, o, w}. The overall mole-fraction of a component, z := (z α ) α=c,w can be expressed in terms of the natural variables as:

z c = y c ρ g S g + x c ρ o S o p ρ p S p (14)
and

z w = ρ w S w p ρ p S p . ( 15
)
The phase-split algorithm is used to solve the system of 2n h + 3 equations and 2n h + 3 variables y c , x c and β p with the pressure, P , and the overall mole-fractions, z, fixed.

z c = y c β g + x c β o (16) z w = β w (17) fc,g (P g , y 1 , ..., y n h ) = fc,o (P g , x 1 , ..., x n h ) ( 18
) p β p = 1.0 (19)
and

c y c = 1.0. ( 20
)
Since p β p = 1.0 and α=c,w z α = 1.0, the last equation is equivalent to c x c = 1.0. To compute the thermodynamic values of the saturations and the mole-fractions, we first compute the overall mole-fractions with ( 14) and ( 15), we then solve the phase-split system to get the thermodynamic values of the mole-fractions y c and x c . Finally, we use the β p to compute the normalized thermodynamic values of the saturations (normalized values needed for the relative permeabilities and the capillary pressures functions):

S thermo p = βp ρp q βq ρq . (21)
The phase-split computations follow the work in [START_REF] Michelsen | Thermodynamic Models: Fundamentals and Computational Aspects[END_REF].

In addition to the n c = n h + 1 conservation equations and the n h + 5 constraint equations, we employ an additional constraint equation and an additional variable. We choose the total-mole density, ρ t , as the additional variable for two reasons. First, ρ t is an intensive variable, but unlike saturation or a mole-fraction, it is not a fraction; second, the phase-split computations are independent of ρ t . Let us define

S thermo p = ρ t β p ρ p (22)
as the 'thermodynamic saturation' of phase p. It is the volume of phase p divided by the pore-volume. At convergence, the sum of the thermodynamic volumes of the fluid phases must equal the pore volume. The usual phase saturation, S p , is referred to as the 'flow saturation'. With the definition of the thermodynamic saturations, the following relationship is always satisfied:

p ρ p S thermo p = p ρ p ρ t β p ρ p = ρ t p β p = ρ t . (23)
The additional constraint equation is the 'thermodynamic volume equation', namely,

p S thermo p = ρ t p β p ρ p = 1.0. ( 24
)
If the thermodynamic volume equation is satisfied, then the thermodynamic saturations and the flow saturations are equal, and the following relationship holds:

ρ t := p ρ p S thermo p = p ρ p S p . (25)

Primary and secondary equations and variables

We have a nonlinear system of 2n h + 7 equations and 2n h + 7 unknowns. The n c = n h + 1 conservation equations are used as the primary equations. In addition to the 'thermodynamic volume equation', we have n h + 5 local constraint equations, which make up the secondary equations. The secondary equations are local constraints and only involve the variables of the control-volume (cell) under consideration. When we linearize the system of equations, the secondary equations are used to eliminate algebraically the n h + 5 variables. After the linear solve that computes the update of the primary variables, the update of the secondary variables are computed from the linearization of the secondary equations.

In our formulation, the gas-phase pressure, P g , is always a primary variable. The remaining n h primary variables depend on the phase state of the control volume. For oil-water cells, the primary variables are: S w and x 1 . . . x n h -1 ; for gas-water cells, the primary variable set is: S w and y 1 . . . y n h -1 ; for three-phase cells, the variables are: S g , S w and y 2 . . . y n h -1 .

The sum of all the component conservation equations leads to:

∂ ∂t φ p ρ p S p + ∇ • p ρ p u p = p ρ p q p . (26)
Substitution of Eq. ( 25) into the accumulation term and Eq. ( 5) in the flux term leads to

∂ ∂t [φρ t ] + ∇ • p ρ p λ p λ t u t + Ψ p = p ρ p λ p λ t q t . (27)
These equations represent the total-mole (mass) conservation equation in the standard form (Eq. ( 26)) and the transport form (Eq. ( 27)). We use one of these equations to replace one of the component conservation equations. Here, we replace the water mole (mass) conservation equation with the total-mole (mass) conservation equation. Even though the rock-fluids system is compressible, the total-mole conservation equation has near-elliptic behavior. The component conservation equations, on the other hand, have near-hyperbolic behavior.

In the following matrix notation, "P " refers to the gas-phase pressure, "R" refers to the total mole density variable, "T " to the total-mole conservation equation, and "C" refers to both the n h primary component variables and the hydrocarbon conservation equations. R C and R T are then the residuals of the hydrocarbon conservation equations and the total-mole conservation equation, respectively. A CC , A CP , A T C and A T P are then the derivatives of the hydrocarbon conservation equations versus the primary component variables, the derivatives of the hydrocarbon conservation equations versus the gas-phase pressure, the derivatives of the total-mole conservation equation versus the primary component variables and the derivatives of the total-mole conservation equation versus the gas-phase pressure, respectively.

Sequential Fully Implicit Scheme for Compositional Flow

We use a finite-volume method with single-point upstream weighting for the spatial discretization and a first-order implicit (backward Euler) scheme for the integration in time.

Fully implicit method

The Fully Implicit (FI) method with the natural-variables formulation [START_REF] Coats | An equation of state compositional model[END_REF] uses the usual form of the conservation equations ( 1) and ( 26) without the thermodynamic volume equation; that is, the variable ρ t is not used. Upon convergence, if ρ t is computed according to equation ( 24) , then we have S thermo p = S p for each phase p ∈ {g, o, w}. Algorithm 1 describes the FI method, in which denotes the convergence tolerance. compute R C and R T from equations ( 1) and ( 26) 1) and ( 26) 1) and ( 26) end Algorithm 1: FI method

while max(|R C |, |R T |) > do compute A CC , A CP , A T C , A T P from equations (
compute δC δP = - A CC A CP A T C A T P -1 R C R T C = C + δC P = P + δP compute phase-split only to detect phase appearance recompute R C and R T from equations (

Nonlinear volume-balance equation

The temporal discretization of the material balance equations ( 1) and ( 2) can be written for α = c, w as:

φ p ρ p S p z α n+1 -φ p ρ p S p z α n = ∆t f n+1 α (28)
The contributions from the fluxes and the wells for a hydrocarbon component c and the water component are:

f n+1 c = -∇ • [y c ρ g u g + x c ρ o u o ] + y c ρ g q g + x c ρ o q o (29) and f n+1 w = -∇ • [ρ w u w ] + ρ w q w ( 30
)
The superscripts n and n + 1 refer to the previous and current timesteps. The Sequential Fully Implicit (SFI) nonlinear volume-balance equation for immiscible flow [START_REF] Jenny | Adaptive fully implicit multi-scale finite-volume method for multi-phase flow and transport in heterogeneous porous media[END_REF] was obtained as follows: each mole (mass) conservation equation was multiplied by the inverse of the molar (mass) density of the phase, and then summed together. When the flow is immiscible, the inverse of the molar density of the phase is the molar volume of that phase. The equivalent of the molar volume of each phase for compositional flow is the partial molar volume of each component. We define the partial molar volume of component α as:

V Tα = ∂V T ∂N α ( 31
)
where V T is the total fluid volume and N α is the number of moles of component α = c, w. In the case of immiscible multiphase flow, where each component can exist in only one phase, the following relationships are satisfied:

c y c V Tc = 1 ρ g , (32)
c x c V Tc = 1 ρ o ( 33
)
and

V Tw = 1 ρ w . ( 34
)
As stated by Watts [START_REF] Watts | A compositional formulation of the pressure and saturation equations[END_REF], the partial molar volumes, the overall mole-fractions and the overall specific volume 1 p ρpSp satisfy the following relationship:

α=c,w z α V Tα = 1 p ρ p S p (35)
In this paper, the pressure equation is obtained as the weighted sum of the nonlinear conservation equations [START_REF] Redlich | On the Thermodynamics of Solutions. V. An Equation of State. Fugacities of Gaseous Solutions[END_REF]; the weights are the partial molar volumes. The result can be written as

α=c,w V n+1 Tα φ p ρ p S p z α n+1 = α=c,w V n+1 Tα φ p ρ p S p z α n + ∆t α=c,w V n+1 Tα f n+1 α (36)
Using relationship [START_REF] Zaydullin | Nonlinear formulation based on an equation-of-state free method for compositional flow simulation[END_REF], it follows that α=c,w

V n+1 Tα φ p ρ p S p z α n+1 = φ n+1 , (37)
and the new nonlinear pressure equation can be written as:

φ n+1 = α=c,w V n+1 Tα φ p ρ p S p z α n + ∆t α=c,w V n+1 Tα f n+1 α (38)
Equation ( 38) is the nonlinear overall volume-balance. This equation describes the pressure field that equilibrates the system if the flux terms are taken as implicit functions of the pressure, but the overall mole-fractions are from the previous composition update. In our sequential-implicit solution strategy, we use the Newton method to iterate out the pressure nonlinearities in the overall-volume balance; thus, the derivatives of the partial molar volumes with respect to pressure are required. During the pressure solution stage, the overall mole-fraction of each component (in the control volume) is kept fixed. The constraint equations ( 8) to ( 13) are used to get the derivatives of y c and x c with respect to P , and equation ( 21) is used to get the derivatives of the saturations with respect to P . After each Newton iteration (pressure solution), a phase-split computation is performed for each cell. Then, the saturations and mole-fractions for each phase are reset to their thermodynamic equilibrium values. Recent research show that the cost of the phase-split calculations can be reduced to a small fraction of the total simulation cost [START_REF] Michelsen | Thermodynamic Models: Fundamentals and Computational Aspects[END_REF][START_REF] Voskov | Compositional space parameterization: Multicontact miscible displacements and extension to multiple phases[END_REF][START_REF] Khait | Operator-based linearization for general purpose reservoir simulation[END_REF]. Equation ( 38) is identical to the pressure equation of the SFI method for the dead-oil (immiscible multiphase flow with no mass transfer), where the phase densities are used for the decoupling. It is also identical to the decoupling used by [START_REF] Lee | Multiscale finite-volume formulation for multiphase flow in porous media: black oil formulation of compressible, three-phase flow with gravity[END_REF][START_REF] Møyner | A multiscale restriction-smoothed basis method for compressible black-oil models[END_REF] for the black-oil formulation (i.e., gas component can dissolve in the oil phase).

Compositional system

The second step of the sequential implicit method consists of freezing the pressure and total-velocity fields and advecting the components. For this purpose, we use the transport form of the conservation equations (i.e., Eqs. ( 6) and ( 27)).

To deal with the transport, two classes of methods have been proposed. The first class uses variables that represent the full fluid content. This is the case for Acs et al. [START_REF] Acs | General purpose compositional model[END_REF], Watts [START_REF] Watts | A compositional formulation of the pressure and saturation equations[END_REF], Trangenstein & Bell [START_REF] Trangenstein | Mathematical structure of the blackoil model for petroleum reservoir simulation[END_REF][START_REF] Trangenstein | Mathematical structure of compositional reservoir simulation[END_REF], Dicks [START_REF] Dicks | Higher Order Godunov Black-Oil Simulations for Compressible Flow in Porous Media[END_REF], Pau et al. [START_REF] Pau | An adaptive mesh refinement algorithm for compressible twophase flow in porous media[END_REF], Faigle et al. [START_REF] Faigle | Multi-physics modeling of non-isothermal compositional flow on adaptive grids[END_REF] and Doster et al. [START_REF] Doster | A robust implicit pressure explicit mass method for multi-phase multi-component flow including capillary pressure and buoyancy[END_REF]. The variables used are: the total number of moles (mass) of each component, or the total number of mole (mass) of each component divided by the pore volume. This is equivalent to updating the variables ξ c and ξ w defined as ξ c = φ (ycρgSg+xcρoSo) p ρpSp ρ t and ξ w = φ ρwSw p ρpSp ρ t and recomputing all the other variables by phase-split computations.

Since these formulations are based on total-mole (total-mass) variables, it is not possible to employ a natural-variables (i.e., saturations-based) formulation to skip the phase-split computations. In the natural-variables formulation, the local equilibrium constraints Eqs. ( 13) are solved simultaneously with the conservation equations ( 6) and ( 27) by the Newton method, and the phase-split computations are only used to determine if a fluid phase has appeared. After convergence of the coupled system of component conservation equations, the moles (mass) of each component are conserved; however, the thermodynamic volume equation [START_REF] Møyner | A multiscale restriction-smoothed basis method for compositional models[END_REF] may not be satisfied exactly. Since this error is local, Acs et al. [START_REF] Acs | General purpose compositional model[END_REF], Watts [START_REF] Watts | A compositional formulation of the pressure and saturation equations[END_REF], Pau et al. [START_REF] Pau | An adaptive mesh refinement algorithm for compressible twophase flow in porous media[END_REF], Faigle et al. [START_REF] Faigle | Multi-physics modeling of non-isothermal compositional flow on adaptive grids[END_REF] and Doster et al. [START_REF] Doster | A robust implicit pressure explicit mass method for multi-phase multi-component flow including capillary pressure and buoyancy[END_REF] proposed to accept the timestep, but to introduce a local relaxation term in order to keep the error bounded in time. Dicks [START_REF] Dicks | Higher Order Godunov Black-Oil Simulations for Compressible Flow in Porous Media[END_REF] studied the thermodynamic volume error and introduced a smoothing algorithm that locally perturbs the material balance. Trangenstein & Bell [START_REF] Trangenstein | Mathematical structure of the blackoil model for petroleum reservoir simulation[END_REF][START_REF] Trangenstein | Mathematical structure of compositional reservoir simulation[END_REF] mentioned that the pressure equation is a linearization of the thermodynamic volume error, and that it could be possible to add an outer-loop on the pressure and the compositional system. However, in their applications, they accept the timestep after one outer-loop. No study exists of the convergence rate of the outer iterations.

The second class of methods [START_REF] Jenny | Adaptive fully implicit multi-scale finite-volume method for multi-phase flow and transport in heterogeneous porous media[END_REF][START_REF] Lee | Multiscale finite-volume formulation for multiphase flow in porous media: black oil formulation of compressible, three-phase flow with gravity[END_REF][START_REF] Hajibeygi | Compositional multiscale finitevolume formulation[END_REF][START_REF] Møyner | A multiscale restriction-smoothed basis method for compressible black-oil models[END_REF][START_REF] Moncorgé | Modified sequential fully implicit scheme for compositional flow simulation[END_REF][START_REF] Møyner | A multiscale restriction-smoothed basis method for compositional models[END_REF][START_REF] Møyner | A multiscale restriction-smoothed basis method for high contrast porous media represented on unstructured grids[END_REF] uses saturations and mole-or mass-fractions. Since the composition variables are only fractions, it is necessary to remove one conservation equation during the solution of the coupled transport. It then means that the pressure equation is not considered a volume-balance, but rather a total-mole or total-mass conservation law. In immiscible situations, with the exception of small time scale transient situations, flow and transport usually can be considered decoupled. But for compositional systems, flow and transport usually cannot be decoupled, and it is thus necessary to perform several outer iterations to reduce to an acceptable tolerance. Recently, this second class of methods has been abandoned, that is, the first class of methods is preferred [START_REF] Moncorgé | Stable sequential fully implicit scheme for compositional simulation[END_REF][START_REF] Møyner | Comparison of splitting schemes for compositional multiscale solvers[END_REF].

Here, we propose to converge the system as in the first class of methods with our new extended natural-variables formulation. The pressure and the total-velocity are frozen. The total mole density is used as an independent variable, and we converge the transport system with the n h + 1 conservation equations ( 6) and [START_REF] Peng | A new two-constant equation of state[END_REF]. The variables are the n h primary component variables and the additional total mole density variable ρ t . The advantage of this formulation -compared with using the number of moles, or mass, per component -is that for the compositional update, we retain the benefits of the natural-variables formulation. Namely, we better control the nonlinearities with the saturations as variables and we compute the phase-split calculations only when a new phase is detected. At convergence, all the components are conserved to the prescribed tolerance; however, the thermo-dynamic volume equation may not be satisfied. We define the total-velocity error as the difference between the total-velocity based on the latest variables after the compositional update and the total-velocity obtained from the pressure solution and frozen during the compositional update. We show that the total-velocity error, which had not been identified before, is present in the invaded regions, and that it is the primary source of nonlinear convergence problems associated with the splitting scheme. This total-velocity error is particularly important in regions experiencing appearance/disappearance of the gas phase during the timestep. The thermodynamic volume residual is defined as

R thermo = p S thermo p -1.0 (39)
and the relative total-velocity residual as

R ut = 1 u t (u f t -u t ), (40)
where u f t is the total-velocity, which is kept fixed during the composition update, and u t is the total-velocity computed by equation ( 4) with the latest updated variables. Upon solving the nonlinear pressure equation, two strategies are possible. The first strategy consists of continuing with outer iterations until the dimensionless residuals |R thermo | ∞ and |R ut | ∞ (infinity norm) fall below a tight tolerance. For strongly coupled flow and transport, even a tolerance of 0.01 may take many outer iterations, or no convergence may be achieved at all. This strategy is not practical, as it requires too many iterations in order to be competitive with methods like mSFI [START_REF] Moncorgé | Modified sequential fully implicit scheme for compositional flow simulation[END_REF]. The second strategy consists of relaxing the tolerances of |R thermo | ∞ and |R ut | ∞ in order to achieve convergence with less outer iterations. Another measure of the total-velocity error is based on the divergence of the total-velocity. The dimensionless divergence of the total-velocity is written as

∇ • u t = j (Q t ) j,i max( (V P ) i ∆t , max j (| (Q t ) j,i |)) (41)
where j are the cells surrounding cell i. (Q t ) j,i are the total volumetric rate from cell j to cell i (positive if entering cell i, negative if leaving cell i), (V P ) i the pore volume of cell i and ∆t the discretized timestep. The divergence of the total-velocity error is defined by the dimensionless residual

R ∇•ut = j Q f t j,i -j (Q t ) j,i max( (V P ) i ∆t , max j (| Q f t j,i |), max j (| (Q t ) j,i |)) (42)
where

Q f t j,i
are the total volumetric rates from cell j to cell i that have been frozen during the composition update.

For challenging test cases, we observe that tolerances of 0.2 for |R thermo | ∞ , 0.4 for |R ut | ∞ and 0.1 for |R ∇•ut | ∞ are good enough for strongly coupled flowtransport problems. From this point on, we can accept the timestep with the converged accumulation terms computed with the quantities ξ c and ξ w . The accepted volume balance error: p S thermo p -S p can be added to the pressure system in order to not propagate this error. However, the volumebalance equation (38) already prevents the propagation of the volume-balance error that potentially exists at time t n by equilibrating volumes at time t n+1 . Adding the volume-balance error does not change the results, but may perturb the solution in some situations. We do not recommend its use.

The algorithm of this second strategy is described in Algorithm 2.

while max(|R C | ∞ , |R T | ∞ ) > or |R thermo | ∞ > thermo or |R ut | ∞ > ut or |R ∇•ut | ∞ > ∇•
while max(|δC| ∞ , |δρ t | ∞ ) > C do
compute Jacobian and residual of equations ( 6) and ( 27) (all conservation equations) wrt variables C and ρ t ; constraints equations ( 8) to ( 13) are used to convert the derivatives wrt secondary variables in function of derivatives wrt primary variables C solve and update C = C + δC and ρ t = ρ t + δρ t compute phase-split only to detect phase appearance compute R C and R T from equations ( 6) and ( 27) end compute R thermo , R ut and R ∇•ut unfreeze total-velocity u t end Exact conservative quantities are ξ c and ξ w Algorithm 2: Compositional SFI method In the case of immiscible multiphase flow with one component per phase, our solution and the solution of Jenny et al. [START_REF] Jenny | Adaptive fully implicit multi-scale finite-volume method for multi-phase flow and transport in heterogeneous porous media[END_REF] are identical. This is not the case for black-oil formulations, where the solution of Lee et al. [START_REF] Lee | Multiscale finite-volume formulation for multiphase flow in porous media: black oil formulation of compressible, three-phase flow with gravity[END_REF] and Møyner & Lie [START_REF] Møyner | A multiscale restriction-smoothed basis method for compressible black-oil models[END_REF] is relaxing the material balance of one component.

Comparison with other compositional pressure equations

Similar to the nonlinear volume-balance equation, we use the gas-phase pressure P and the overall mole-fractions z α for each component α = c, w. The overall mole-fractions z α are fixed, and the constraint equations ( 8) to ( 13) are used to get the derivatives of y c and x c wrt P . Equation ( 21) is used to get the derivatives of the saturations with respect to P . Sequential formulations for compositional flows where derived in a semi-implicit context by Acs et al. [START_REF] Acs | General purpose compositional model[END_REF], Watts [START_REF] Watts | A compositional formulation of the pressure and saturation equations[END_REF], Trangenstein & Bell [START_REF] Trangenstein | Mathematical structure of the blackoil model for petroleum reservoir simulation[END_REF][START_REF] Trangenstein | Mathematical structure of compositional reservoir simulation[END_REF], Dicks [START_REF] Dicks | Higher Order Godunov Black-Oil Simulations for Compressible Flow in Porous Media[END_REF], Pau et al. [START_REF] Pau | An adaptive mesh refinement algorithm for compressible twophase flow in porous media[END_REF], Faigle et al. [START_REF] Faigle | Multi-physics modeling of non-isothermal compositional flow on adaptive grids[END_REF] and Doster et al. [START_REF] Doster | A robust implicit pressure explicit mass method for multi-phase multi-component flow including capillary pressure and buoyancy[END_REF]. In these previous methods, the pressure equation is derived by linearizing the difference between the sum of the thermodynamic volumes of each phase and the pore volume. The resulting pressure equation accounts for the change of volume of the fluids due to the advection of the compositions. We refer to [START_REF] Watts | A compositional formulation of the pressure and saturation equations[END_REF] for the derivation of the 'linearized volume-balance pressure equation' that can be written in semi-discretized in time form

dφ dP + φ p ρ p S p ∂ p ρ p S p ∂P (P n+1 -P n ) = ∆t α=c,w V Tα f n+1 α . (43)
In a semi-implicit situation, we would first converge

dφ dP + φ p ρ p S p ∂ p ρ p S p ∂P n (P n+1 -P n ) = ∆t α=c,w V n Tα f n+1 α ( 44
)
and then switch to the nonlinear composition update. In a sequential implicit situation we would converge the equation

dφ dP + φ p ρ p S p ∂ p ρ p S p ∂P ν (P ν+1 -P n ) = ∆t α=c,w V ν Tα f ν+1 α (45)
until |P ν+1 -P ν | is small, where the superscripts ν and ν + 1 denote the old and new pressure iteration levels. Coats [START_REF] Coats | A note on IMPES and some IMPES-based simulation models[END_REF] built a pressure equation by algebraic manipulations of the usual material balance equations, and Møyner et al. [START_REF] Møyner | A multiscale restriction-smoothed basis method for high contrast porous media represented on unstructured grids[END_REF][START_REF] Møyner | Comparison of splitting schemes for compositional multiscale solvers[END_REF] use this pressure for their sequential implicit method for compositional flow. Coats' pressure equation in a sequential implicit setting is equivalent to converging the equation α=c,w

V ν Tα φ p ρ p S p z α ν+1 = α=c,w V ν Tα φ p ρ p S p z α n + ∆t α=c,w V ν Tα f ν+1 α ( 46
)
until |P ν+1 -P ν | is small, where the superscripts ν and ν + 1 denote the old and new pressure iteration levels, and z α in the left-hand-side term is taken at the latest composition iteration level. The overall compressibility of the system can dramatically change during the transition from a state with only liquid phases present to a state where a gas phase is present. Moreover, V Tα can be very strong functions of pressure as seen in the relationships (32), ( 33) and [START_REF] Young | Continuous compositional volume-balance equations[END_REF]. It is then possible to encounter cases where the pressure equations ( 45) and ( 46) do not converge. It is then necessary to use both the derivatives of the compressibility and the partial molar volumes. The linearized equation ( 45) and the nonlinear volume-balance equation (38) are only equivalent for incompressible immiscible systems.. The Coats pressure equation ( 46) and the nonlinear volume-balance equation (38) are different because V Tα depend on pressure, and the left-hand-side of (46) is only equals to φ n+1 at convergence. However, when the Coats' pressure equation converges, the solution is the same as the nonlinear volume-balance equation.

Results and Discussion

The nonlinear strategies used for the numerical experiments with the FI method, the volume balance system, and the compositional SFI method ensure that the variations of pressure, saturations and mole-fractions over a Newton iteration get damped to respectively 1000.0 Psi, 0.2 and 0.2, and that the saturations and mole-fractions are kept between 0.0 and 1.0. If the value of the total mole density ρ t after a Newton iteration is negative, it is set to zero, and if the variation of ρ t over a Newton iteration is larger than the maximum allowed value, i.e., if |δρ t | ∞ > max p ρ p , the value of ρ t is set to ρ t = p ρ p S p .

1D test cases

In this section we test our compositional SFI algorithm using 1D models; first for dead-oil and then for live-oil fluids.

Dead-oil

First, the compositional SFI algorithm is tested for cases with no mass transfer between the fluid phases. We consider injection of water into a reservoir saturated with dead-oil. With a CFL number of 140 [START_REF] Coats | IMPES stability: Selection of stable timesteps[END_REF], the FI method converges in 16 Newton iterations for the first timestep. The SFI method was stopped after one outer iteration of Algorithm 2. The first outer iteration of the SFI method consists of three pressure iterations and 16 compositional iterations. The results are shown is Figure 1a. The saturation fronts obtained with the FI and SFI methods are in very good agreement, but the pressure profiles are quite different. For the SFI method, the material balance equations are satisfied, but the splitting errors

|R thermo | ∞ , |R ut | ∞ and |R ∇•ut | ∞ are,
respectively, 2e-3, 0.55 and 0.04. Figure 2 shows the profiles of the three splitting errors. R thermo is very small, R ut has large values across the 1D domain, and R ∇•ut is around 1% except for some localized areas. The total-velocity errors have an impact on the absolute value of the pressure; however, since the system is immiscible and the divergence of the totalvelocity is correctly reconstructed, the transport is already correct during the first outer iteration. Two outer SFI iterations (a total of 6 pressure iterations and 18 compositional iterations) are required to fully converge the system with |R thermo | ∞ , |R ut | ∞ and |R ∇•ut | ∞ set to 1e-5, 1e-2 and 6e-3, respectively. These results are shown in Figure 1b, and one can observe that the profiles obtained with the FI and SFI methods are identical.

In the second test case, injection of dry-gas into a dead-oil reservoir is considered. The gas, which is less dense, less viscous, and more compressible than the oil phase, is injected with a CFL number of 780. The ratio of the oil to gas mass densities is 735; the ratio of the oil to gas viscosities is 8.1, and the ratio of the gas to oil compressibilities is 0.005. For the first timestep, the FI method converges in 42 Newton iterations. The SFI method was stopped after one outer iteration. The results are shown in Figure 3a are required to converge the SFI method. The results are shown in Figure 3b, and it can be observed that the profiles obtained with the FI and SFI methods are identical.

As a last dead-oil case, gas injection with gravity is considered. Gas is injected into a reservoir initially filled with pure oil, and a CFL number of 65 was applied. The FI method converged in 16 Newton iterations for the first timestep. The SFI method was stopped after one outer iteration. The results are shown in Figure 4a

Depletion of a live-oil reservoir

Now that the algorithm is validated for dead-oil cases, we consider a depletion case with the SPE 5 fluid [START_REF] Killough | Fifth comparative solution project: Evaluation of miscible flood simulators[END_REF]. The reservoir is originally saturated with pure oil, and we drop the pressure at the right-side. In one timestep with a CFL number of 56, gas appears on the right side of the reservoir. The FI method requires seven iterations to converge. One outer iteration of the SFI method consisting of six pressure and three compositional iterations was employed (|R thermo | ∞ , |R ut | ∞ and |R ∇•ut | ∞ are 0.36, 0.07 and 0.05). The pressure and gas saturation profiles for the FI and SFI methods are shown in Figure 5. The SFI method converges to the correct solution in one outer iteration even though the thermodynamic volume error of 0.36. However, R thermo is large for a single cell where the gas phases appears. In two outer iterations (consisting of nine pressure and four compositional iterations), |R thermo | ∞ , |R ut | ∞ and |R ∇•ut | ∞ decrease to, respectively, 1e-3, 2e-3 and 2e-3.

Gas injection in a live-oil reservoir

This test case starts with the initial state shown in Figure 1b. Live-gas is injected into the reservoir. The SPE 5 fluid [START_REF] Killough | Fifth comparative solution project: Evaluation of miscible flood simulators[END_REF] is used, whereas the light component is injected. Figure 6a shows the pressure and saturation profiles after one timestep with a CFL number of 60 obtained using the FI and the SFI methods. The FI method required 22 iterations to converge. One outer iteration of the SFI method consisting of three pressure and 11 compositional iterations was employed (|R thermo | ∞ , |R ut | ∞ and |R ∇•ut | ∞ are 4.95, 0.89 and 0.57). The SFI solutions with one outer iteration are different from the FI results. For complete convergence, 27 iterations (three outer, 9 pressure and 18 compositional iterations) were required; the profiles are shown in Figure 6b.

The next test case describes injection of wet-gas into a reservoir initially filled with live-oil. Figure 7 shows pressure and gas saturation profiles after one timestep with a CFL number of 64 using the SFI method; once after one outer iteration consisting of three pressure and 11 compositional iterations and once after full convergence (requiring in total three outer iterations consisting of 9 pressure and 20 compositional iterations). The FI method requires 37 iterations to converge. After one outer iteration of the SFI method the spiltting-errors

|R thermo | ∞ , |R ut | ∞ and |R ∇•ut | ∞ are 10
.2, 0.90 and 0.52, respectively.

The next test case starts with the initial state of Figure 8a and live-gas is injected into the reservoir. Figure 8b shows pressure and gas saturation profiles after one timestep with a CFL number of 76 for the FI method and the SFI method stopped after one outer iteration. The FI method requires 10 iterations to converge. One outer iteration of the SFI method requires three pressure and 9 compositional iterations. The saturation profiles look identical, but the pressure profiles are different. Figure 9 shows the thermodynamic volume error R thermo , the total-velocity error R ut and the divergence of the total-velocity error R ∇•ut profiles for the SFI method. |R thermo | has values up to 8%, |R ut | has values up to 34% and |R ∇•ut | has values up to 14%. On the contrary to the total-velocity error R ut that has a support on the whole 1D domain, R thermo and R ∇•ut have only non-zero values in localized areas. This is an indication that when the compressibility and thermodynamic effects are not too strong we can reconstruct the correct transport in one outer iteration. Table 1 shows the evolution for the SFI method of the thermodynamic volume error, the total-velocity error, the divergence of the total-velocity error as well as the divergence of the total-velocity in function of the number of outer iterations. We observe that the splitting errors

|R thermo | ∞ , |R ut | ∞ and |R ∇•ut | ∞ decrease
with the outer iterations. At some point, the divergence of the total-velocity, ∇ • u t , converges to the true value, Figure 6: Gas injection after water injection: pressure and saturations profiles versus cell number for the FI and the SFI methods.

Table 1: Gas injection: evolution of the thermodynamic volume error, the total-velocity error, the divergence of the total-velocity error and the divergence of the total-velocity in function of the number of outer iterations. since the flow is compressible.

Starting from the final state of Figure 7, water is injected for one timestep with a CFL number of 190. Figure 10 shows the pressure and saturations profiles after this timestep. The FI method required 18 iterations to converge. One outer iteration of the SFI method required three pressure and 16 compositional iterations. After this outer iteration, the maximum values of |R thermo | ∞ , |R ut | ∞ and |R ∇•ut | ∞ were 0.34 (in cell 39), 0.89 (in cell 1) and 0.05 (in cell 1), respectively. The pressure and gas-saturation profiles obtained after one iteration with the SFI and FI methods are different. To achieve convergence, SFI scheme required two outer iterations consisting of a total of 6 pressure and 18 compositional iterations.

Similar test cases have been simulated with the "system 1" fluid from [START_REF] Zaydullin | Nonlinear formulation based on an equation-of-state free method for compositional flow simulation[END_REF] and the decane/CO 2 /methane fluid from [START_REF] Møyner | A multiscale restriction-smoothed basis method for compositional models[END_REF], and the same observations apply.

2D test cases

In this section, we present 2D test cases with gas and water injection.

2D Case-1: Gas and Water Injection

The first test case is from Moncorgé et al. [START_REF] Moncorgé | Modified sequential fully implicit scheme for compositional flow simulation[END_REF]. The light component is injected into multi-component oil at the lower-left corner of the model; the producer, for which the bottom-hole pressure (BHP) is set to 4000 Psi, is located at the upper-right corner of the model. See figure 11a. We study the behaviors for one particular timestep. Figure 11 shows the solution of the SFI method after one outer iteration. The mass balance equations are converged, but the residuals R thermo and R ut are 'relaxed'. The CFL numbers for explicit treatment of saturation (CFL S ) and compositions (CFL X ) [START_REF] Coats | IMPES stability: Selection of stable timesteps[END_REF] are 33 and 224 for the timestep, respectively. The FI method converges in 6 iterations. For the SFI method, we summarize in table 2 the evolution of the errors in the thermodynamic volume, total-velocity, and divergence of the total-velocity, as functions of the number of outer iterations. We observe that the splitting errors decrease with the outer iterations. Since the flow is compressible, the divergence of the total-velocity ∇ • u t converges to the correct nonzero value.

Table 2: Evolution of the thermodynamic volume error, the total-velocity error, the divergence of the total-velocity error and the divergence of the total-velocity as functions of the number of outer iterations for 2D Case-1. After one outer iteration with two pressure and six compositional solves, we observe highly localized errors at the front between the gas and the oil phases in the thermodynamic volume equation with values that amount up to 7%. The cells with errors of more than 1% for the volume are plotted in Figure 11b. Total-velocity errors as large as 29% are observed, and they have a large support as seen in Figure 11c. The cells with errors of more than 5% in the total-velocity are plotted in Figure 11d. Figure 11e shows that errors greater than 5% in the total-velocity divergence are localized. Nevertheless, the largest error is about 44%. For this timestep, we required three outer iterations (a total of six pressure and 10 compositional solves) to converge to the FI solution.

We now study the behaviors for the full simulation. The first period involves injection of the light-component at a rate of six MSCF/Day; the second period entails water injection at a rate of four STB/Day, and the third period involves light-component injection at a rate of six MSCF/Day. This injection strategy leads to complex behaviors with gas fronts chased by water fronts that force the gas to dissolve back into the oil phase. FI requires 399 and SFI 401 timesteps to reach convergence. Figures 12a, 12b and12c show the gas saturation profiles at the end of each of the three periods. For the first period, the average CFL numbers for the explicit treatment of saturation (CFL S ) and for the explicit treatment of compositions (CFL X ) are 40 and 71, respectively, and the corresponding maximum values are 204 and 228. For the second period, the average values of CFL S and CFL X are 39 and 15, respectively, and their maximum values are 500 and 136. For the third period, the average values of CFL S and CFL X are 182 and 191, respectively, and the maximum values are 635 and 490. Figure 13 shows the cumulative number of Newton iterations for the FI method and the cumulative number of pressure and composition iterations for the SFI method, as well as the percentage of gas in the model for both the FI and SFI methods. For this test case, the SFI method requires 26% more pressure iterations and 55% more composition iterations than the FI method Newton iterations. The percentage of gas is exactly the same for both models. Figure 14 shows the cumulative number of CPR(AMG,ILU0)-FGMRES iterations required by the FI method and by the pressure solve of the SFI method and the cumulative number of ILU0-GMRES iterations required by the SFI method for the compositional system. Both the full system of the FI method and the pressure systems of the SFI method require around 5 CPR iterations to converge, but the compositional systems of the SFI method only require 2.2 ILU0-GMRES iterations to converge. As a result, for this test case the SFI method takes 33% more linear pressure iterations than the FI method, but only 65% of the FI linear composition iterations. This test case is quite challenging for SFI methods, as the coupling between flow and transport is strong.

2D Case-2: Gas and Water Injection

The second test case is challenging for the FI method, but easier for the SFI method. Water is injected at the top-left corner at a rate of four STB/Day and gas at the bottom-left corner at a rate of one MSCF/Day. Production occurs at the top-right and bottom-right corners, where the BHP is set to 4000 Psi. To increase the velocity of the water phase, the porosity in the 20 top layers (out of 60 layers) is reduced by a factor of 100. Note that such contrasts are realistic and represent low porosity corridors. Both the FI Figure 13: Cumulative Newton iterations for the FI method and cumulative pressure and composition iterations for the SFI method for 2D Case-1. The percentage of gas in the model for the FI and SFI methods is also shown. and SFI methods require 19 timesteps to converge. Figures 15a and15b show the water and gas saturation maps at the end of the simulation. A large part of the reservoir has been flooded by water and experiences very fast velocities; gas is only present in the lower part of the reservoir. While there is no gas-phase in the reservoir initially, at the end of the simulation almost 8% of the cells have gas. The average CFL number per timestep is 480 and the corresponding maximum value is 1140. Figure 16 shows the cumulative number of Newton iterations for the FI method and the cumulative number of pressure and composition iterations for the SFI method as well as the percentage of gas in the model for both the FI and SFI methods. For this case, the SFI method requires 36% fewer pressure iterations and 17% more composition iterations than the FI method Newton iterations. The percentage of gas for the models is exactly the same. Figure 17, shows the cumulative number of CPR(AMG,ILU0)-FGMRES iterations for the FI method and the pressure solves of the SFI method and the cumulative number of ILU0-GMRES iterations for the compositional system of the SFI method. The full system of the FI method requires in average around 8 CPR iterations to converge, while the pressure systems of the SFI method requires around 6 CPR iterations to converge. The compositional systems of the SFI method requires in average 5.3 ILU0-GMRES iterations to converge. As a result, the SFI method for this test case takes 49% of the linear pressure iterations and 78% of the linear composition iterations of the FI method. The strong oil-water flow is challenging for the FI method and very well handled by the SFI method, and the mild compositional effects are handled well by both methods. This is a configuration that is closer to real reservoir simulation scenarios.

3D test cases

The 3D model here represents an anticlinal formation. The fluid is modeled as a two-component compositional black-oil formulation. There are four water injectors and six producers. Figure 18 shows the model and the gas saturation at the end of the simulation. Figures 19 and20 show the injection and production profiles for each phase as well as the average reservoir pressure. We observe that the FI and SFI produce exactly the same results. The FI method converges in 304 timesteps and the SFI method in 361 timesteps. Figure 21 shows the cumulative Newton iterations for the FI method, the cumulative pressure and composition iterations for the SFI method as well as the percentage of gas in the model for both the FI and SFI methods. We see that the percentage of gas in the model changes quite rapidly, and that Figure 17: Cumulative CPR iterations for the systems of the FI method and the pressure of the SFI method and cumulative ILU0-GMRES iterations for the composition system of the SFI method for 2D Case-2. both methods capture exactly the same physics. These rapid changes are due to the complexity of the medium with strong contrasts between the geological facies. The SFI method requires 20% fewer pressure iterations than the FI method requires Newton iterations and 15% more composition iterations than the FI method requires Newton iterations. Figure 22, shows the cumulative number of CPR(AMG,ILU0)-FGMRES iterations for the FI method and the pressure solves of the SFI method and the cumulative number of ILU0-GMRES iterations for the compositional system of the SFI method. The CPR(AMG,ILU0)-FGMRES system of the FI method requires in average 4.2 iterations, while the pressure system of the SFI method requires in average 3.1 linear iterations and the composition system of the SFI method requires in average 2 linear iterations per system. As a result, the number of linear iterations from the pressure systems of the SFI method amount to 58% of the number of CPR iterations for the coupled system of the FI method, and the number of linear iterations from the composition systems of the SFI method amount to 55% of the number of CPR iterations for the coupled system of the FI method.

Conclusions

We proposed a new compositional Sequential Fully Implicit (SFI) solution scheme with a nonlinear pressure equation for compositional flow simulation. The compositional SFI scheme is an extended natural-variables formulation. The pressure equation used here is obtained as a weighted sum of the nonlinear component conservation equations, and it can be seen as a nonlinear overall-volume balance. This pressure equation is different from the 'usual' pressure equations used in sequential formulations, which are obtained as a linear combination of the linearized component conservation equations or by linear combination of the component conservation equations with constant coefficients. The nonlinear volume-balance equation simplifies to the usual volume-balance equations used by SFI methods for dead-oil and black-oil systems. There are two splitting errors associated with the compositional SFI scheme: (1) volume balance errors, and (2) total-velocity errors. Our analysis indicates that the volume-balance errors are highly localized around fronts where a phase change takes place; however, the errors in the totalvelocity field can be quite large with large spatial support that spans the invaded regions of the reservoir. The total-velocity error associated with the sequential-implicit split is the primary cause of nonlinear convergence problems for SFI schemes including ours. Here, we demonstrate that the splitting errors can be reduced to arbitrarily small values when a nonlinear pressure equation is used. Moreover, we show how to effectively control the splitting errors in order to achieve convergence to any desired tolerance. We tested the compositional SFI approach for challenging test cases in 1D, 2D and 3D. For strongly coupled problems, the SFI algorithm requires a few more iterations than the FI method, but for weakly coupled cases, the sequential algorithm requires fewer iterations than the FI method. Overall, the compositional SFI scheme enjoys convergence properties that are very competitive with the fully-implicit method, even when the coupling between the flow and the transport is quite strong. The SFI formulation has several advantages, however. With the sequential formulation, the parabolic and hyperbolic operators of compositional multi-component flow and transport are now identified and decoupled. As a result, it is now possible to design specific numerical methods optimized for each sub-problem. That includes the use of specific linear solvers for the parabolic and hyperbolic operators, as well as opening the door for advanced (high-order) spatial and temporal discretization schemes. Finally, this new SFI scheme is well suited for multi-scale compositional formulations that rely on sequential coupling of the flow and transport problems. • µ p phase dynamic viscosity in cP

• k rp relative permeability for phase p

Appendix A: Description of the Models

The 1D model used in this study has 220 cells of dimensions of dx = 20 ft, dy = 10 ft, dz = 2 ft. The porosity is φ = 0.2 and the permeability K = 100 mD. The 2D model has 220×60 cells with porosity and permeabilities from the top layer of SPE comparative solution project [START_REF] Christie | Tenth SPE comparative solution project: A comparison of upscaling techniques[END_REF]. The rock compressibility for the 1D and 2D cases is 1.78e-5 1/Psi. The relative permeabilities are quadratic for both the oil and the gas phases. The capillary pressures are zero. The fluid is taken from the SPE 5 comparative solution project [START_REF] Killough | Fifth comparative solution project: Evaluation of miscible flood simulators[END_REF]. The test cases have initial conditions with P = 4000.0 Psi, T = 160 • F and S o = 1.0; table 3 provides the initial composition values. The connection Table 3: Initial oil composition for SPE 5 fluid.

  et al., Watts, Trangenstein & Bell, Dicks, Pau et al., Faigle et al. and Doster et al. consider the pressure equation as a 'volume balance' equation; in the second step of solving the transport, they keep the conservation equations of all the components. Upon convergence, the mass conservation equation of each of the components is satisfied subject to the desired tolerance; however, some discrepancies in the overall-volume balance and the total-velocity persist. Only the overallvolume balance splitting error has been defined in these previous works. Attempts to reduce this volume splitting error have been done in Acs et al., Trangenstein & Bell, Pau et al., Faigle et al. and Doster et al. by using a local relaxation term to keep the error bounded in time or by local smoothing by Dicks. They all result in local changes relaxing the mass balance equations.

  . The figure shows that the pressure and saturation profiles obtained with the FI and SFI methods are quite different. The first outer iteration of the SFI method consists of eight pressure iterations and seven compositional iterations (|R thermo | ∞ , |R ut | ∞ and |R ∇•ut | ∞ are 0.40, 0.81 and 0.12). Three outer SFI iterations consisting of 18 pressure and 12 compositional iterations

  (a) SFI profiles after one outer iteration and converged FI profiles. (b) Converged SFI and FI profiles.

Figure 1 :

 1 Figure 1: Water injection: pressure and water saturation profiles versus cell number for the FI and the SFI methods. 22

Figure 2 :

 2 Figure 2: Water injection: thermodynamic volume error R thermo , totalvelocity error R ut and divergence of the total-velocity error R ∇•ut for the SFI method after one outer iteration.

  , which depicts pressure and saturation profiles as functions of depth. The pressure and saturation profiles obtained with the FI and the SFI methods are quite different. The first outer iteration of the SFI method consists of 14 pressure iterations and eight compositional iterations (|R thermo | ∞ , |R ut | ∞ and |R ∇•ut | ∞ are 0.25, 0.73 and 0.16). Three outer SFI iterations consisting of 27 pressure and 13 compositional iterations are required to converge the SFI method. For this test case with initially pure oil in the reservoir, it is not possible to reach convergence if the usual linearized pressure equation (45) or the Coats pressure equation (46) are employed. The results are shown in Figure 4b, and it can be observed that the profiles obtained with the FI and SFI methods are identical.

Figure 3 :

 3 Figure 3: Injection of dry-gas: pressure and gas saturation profiles versus cell number for the FI and the SFI methods. 24

Figure 4 :

 4 Figure 4: Injection of dry-gas with gravity: pressure and gas saturation profiles versus cell number for the FI and the SFI methods. 25

Figure 5 :

 5 Figure 5: Depletion of live oil: pressure and gas saturation profiles versus cell number for the FI and the SFI methods.

  (a) Initial SFI and FI profiles.(b) SFI profiles after one outer iteration and converged FI profiles for a timestep with a CFL number of 76.

Figure 8 :

 8 Figure 8: Gas injection: pressure and gas saturation profiles versus cell number for the FI and the SFI methods. 31

Figure 9 :

 9 Figure 9: Gas injection: thermodynamic volume error R thermo , total-velocity error R ut and divergence of the total-velocity error R ∇•ut for the SFI method after one outer iteration.

Figure 10 :

 10 Figure 10: Water injection: pressure and saturations profiles versus cell number for the FI and the SFI methods. 33

  (a) S g . (b) Volume error above 1%. Error ranges from -7% to +7%. (c) Total-velocity error reaches up to 30%. (d) Velocity error exceeding 5%.(e) Divergence of velocity error exceeding 5%. Error ranges from -39% to +44%.

Figure 11 :

 11 Figure 11: Maps for 2D Case-1 after application of one outer iteration of the sequential algorithm. 35

  (a) After gas injection at the end of the first period. (b) After water injection at the end of the second period.(c) After gas injection at the end of the third period.

Figure 12 :

 12 Figure 12: Gas saturation for 2D Case-1.

Figure 14 :

 14 Figure 14: Cumulative CPR iterations for the systems of the FI method and the pressure of the SFI method and cumulative ILU0-GMRES iterations for the composition systems of the SFI method for 2D Case-1.

  (a) Water saturation at the end of the simulation. (b) Gas saturation at the end of the simulation.

Figure 15 :

 15 Figure 15: Saturation maps for 2D Case-2.

Figure 16 :

 16 Figure16: Cumulative Newton iterations for the FI method and cumulative pressure and composition iterations for the SFI method for 2D Case-2. The percentage of gas in the model for the FI and SFI methods is also shown.

Figure 18 :

 18 Figure 18: Gas saturation map for the anticlinal model. The water injectors are the blue arrows and the producers are the green arrows. Plotted with ResInsight [25].

Figure 19 :

 19 Figure 19: Liquid rates (water injection/production and oil production) profiles for the FI and SFI methods for the anticlinal model.

Figure 20 :

 20 Figure 20: Average reservoir pressure and gas production rate profiles for the FI and SFI methods for the anticlinal model.

Figure 21 :

 21 Figure 21: Cumulative Newton iterations for the FI method and cumulative pressure and composition iterations for the SFI method for the anticlinal model. The percentage of gas in the model for FI and SFI methods are also shown.

Figure 22 :

 22 Figure 22: Cumulative CPR iterations for the FI method and cumulative linear iterations for the pressure and composition systems for the SFI method for the anticlinal model.

•

  φ porosity of the rock • K permeability in mD • D depth in ft • g gravitational acceleration in ft 2 -Psi/lb • n h number of hydrocarbon components • M W c molar fraction of component c • P p pressure of phase p in Psi • P c GO (S g ) = P g -P o capillary pressure in Psi • P c OW (S w ) = P o -P w capillary pressure in Psi • S p flow saturation of phase p • S thermo p thermodynamic saturation of phase p • S thermo p normalized thermodynamic saturation of phase p • x c mole-fraction of component c in the oil phase • y c mole-fraction of component c in the gas phase • z α total mole-fraction of component α = c, w • β p mole-fraction of phase p • ρ p phase mole density in lbmol/ft 3 • ρ t total mole density in lbmol/ft 3 • ρ p phase mass density in lb/ft 3 , ρ g = ( c y c M W c )ρ g , ρ o = ( c x c M W c )ρ o and ρ w = M W w ρ w .

  0.03 0.07 0.2 0.15 0.05 factors for all the cell-well connections are fixed at 0.3 Rbbl.cP/day-Psi. The fluid is produced at 1800 Psi during the depletion.

  al., SPEJ 1985; Watts, SPEJ 1986; Trangenstein & Bell, SIAM 1989; Pau et al., Comp. Geosci. 2012; Faigle et al., Comp. Meth. App. Mech. Eng. 2014; and Doster et al., CRC

  Acs et al., Watts, Trangenstein & Bell, Pau et al., Faigle et al. Doster et al. as well as the pressure equations of Coats and Young

  ut do freeze z α=c,w while |δP | ∞ > P do compute derivative wrt P and residual of volume-balance equation (38) with derivatives of y c and x c given by constraints (8) to (13) and derivatives of S p by (21) solve and update P = P + δP compute phase-split and set y c , x c and S p to thermodynamic values with (16) to (21) recompute residual of volume-balance equation (38) end unfreeze z α=c,w compute and freeze total-velocity u t

(a) SFI profiles after one outer iteration and converged FI profiles. (b) Converged SFI and FI profiles.

(a) SFI profiles after one outer iteration and converged FI profiles. (b) Converged SFI and FI profiles.Figure 7: Gas injection: pressure and gas saturation profiles versus cell number for the FI and the SFI methods.

• fc,p fugacity of hydrocarbon component c in phase p in Psi • u p velocity of phase p in ft/day • u t total-velocity in ft/day • u f t frozen total-velocity during the composition update in ft/day

• λ p mobility for phase p in cP -1

• λ t total mobility in cP -1

• q p and q t well contributions (source or sink) in lbmol/day

• R thermo thermodynamic volume splitting error

• R ut normalized total-velocity splitting error

• R ∇•ut normalized divergence of total-velocity splitting error • V P discretized pore volume in ft 3

• Q t total volumetric rate in ft 3 /day

• V Tα partial molar volume of component α = c, w in ft 3 /lbmol

• ξ α density of component α = c, w in lbmol/ft 3

Acknowledgments

The authors would like to thank TOTAL management for permission to publish this work and O. Møyner for constructive discussions.


E Baudin

M E Hayden

G Tastevin

P J Nacher

NMR Time Reversal Experiments in Highly Polarised Liquid 3 He-4 He Mixtures

Keywords: numbers: 76.60.Jx, 76.60.Lz, 82.56-b, 82.56.Jn, 67.65+z, 75.40-Gb

Long-range magnetic interactions in highly magnetised liquids (laser-polarised 3 He-4 He dilute mixtures at 1 K in our experiment) introduce a significant non-linear and non-local contribution to the evolution of nuclear magnetisation that leads to instabilities during free precession. We recently demonstrated that a multi-echo NMR sequence, based on the magic sandwich pulse scheme developed for solid-state NMR, can be used to stabilise the magnetisation against the effect of distant dipolar fields. Here, we report investigations of echo attenuation in an applied field gradient that show the potential of this NMR sequence for spin diffusion measurements at high magnetisation densities.

INTRODUCTION

Long-range dipolar interactions between the local magnetisation and that of the remainder of a liquid sample play a significant role in determining nuclear spin dynamics at high magnetisation densities. They introduce a nonlocal contribution to the evolution of magnetisation and give rise to complex non-linear behaviour. Spectacularly abrupt decays of NMR free induction signals have been observed at large tip angles in several laser-polarised liquid systems [1], and have been explained by the development of unstable inhomogeneous magnetisation patterns that grow exponentially in time [2]. The onset of precession instabilities is a generic feature of highly magnetised systems, and is related to turbulent spin motion [2]. It is potentially relevant to understanding the dynamics of polarised quantum fluids, and is becoming increasingly important in high-field 1 H-based NMR spectrometry where chaotic spatiotemporal behaviour occurs under the joint action of distant dipolar fields and strong radiation damping [3].

We recently demonstrated that the magic sandwich (MS), a sophisticated pulse sequence designed to refocus the magnetisation lost through broadening by local dipolar fields in solids, can be successfully used to fight the deleterious effect of distant dipolar fields in highly magnetised liquids and to investigate the spatiotemporal development of instabilities [4]. We also demonstrated that it is possible to achieve dynamic stabilisation of transverse precession in laser-polarised He mixtures using a series of MS cycles to periodically refocus the dipolar interactions. In this paper, we report preliminary results that establish the potential of multiple-echo MS-based NMR sequences for spin diffusion measurements at high magnetisation densities, and confirm the fact that conventional spin echo techniques completely fail to achieve this goal [5].

EXPERIMENTAL SETUP

Polarised samples are prepared by condensation of optically pumped 3 He gas into superfluid 4 He in a 1 metre-long Pyrex glass cell. 3 He gas is continuously injected (at mbar pressures with a flow rate of order 0.1 µmol/s) into a room temperature volume where laser optical pumping is performed, as described elsewhere [6]. It then flows down a narrow tube into a 0.44 cm 3 spheroidal volume partly filled with liquid 4 He at temperature T ≃ 1.15 K. Cs coatings are used to avoid wall relaxation, and bulk dipole-dipole relaxation times reach several hours in our dilute samples ( 3 He molar fractions X = 1-5%) [6]. By adjusting the 3 He molar fraction and nuclear polarisation (that can be as high as 40%), we control the magnetisation density M in the sample. In the following, M is conveniently reported in terms of a dipolar frequency F dip = γB dip /2π associated with the characteristic dipolar field B dip = µ 0 M , where µ 0 is the permeability of free space and γ the 3 He gyromagnetic ratio (γ/2π = 32.4 MHz/T). NMR excitation by 90 • rf pulses induces an irreversible loss of this out-of-equilibrium magnetisation, so that a new batch of polarised gas has to be condensed for each experiment.

The cryogenic apparatus includes a commercial aluminum and fiberglass dewar and a homemade non-magnetic insert. The lower end of the cell lies inside a stainless-steel vacuum can, and is thermally anchored to a copper 1 K-pot via bundles of enameled copper wires that do not allow for eddy current loops.

Low frequency NMR (74 kHz) is performed in a 2.3 mT magnetic field (shimmed to 20 ppm over the sample) using crossed coils. The receive coil is capacitively tuned with a very low Q factor (≃ 2) to avoid radiation damping. The transmit coil is actively shielded to avoid Joule heating by eddy currents induced in the surrounding metallic parts during intense cw rf excitation [7]. With the shield coil in place, the rf field outside of the transmit coil is dramatically reduced and, e.g., the measured heat load on the 1 K-pot is more than two orders of magnitude lower. Shielding also minimises field distortions and improves rf homogeneity (±200 ppm over the sample). In-house PC-controlled analog hardware allows for flexible control of hard pulse sequences.

ECHO TRAINS AND DIFFUSION-INDUCED ATTENUATION

Spin diffusion coefficients are usually measured by monitoring the attenuation of echo amplitudes in the presence of a uniform applied gradient. Conventional spin echo techniques are inefficient in highly magnetised liquids. The (quadratic) effect of dipolar interactions actually remains unchanged throughout the sample when a 180 • rf pulse is applied to rotate the local magnetisation. The decay of the average transverse magnetisation is thus expected to remain dominated by dipolar-induced instabilities at high magnetisation densities. Indeed, early attempts to generate single Hahn spin echoes and Carr-Purcell spin echo trains in highly polarised 3 He- 4 He dilute mixtures have led to the observation of multiple echoes [8] and of non-exponential decay [5], respectively.

Figure 1 (traces a and c), spin echoes have the shapes that are expected given the sample geometry and the gradient strengths. Mono-exponential decays of the echo amplitudes are observed. At high magnetisations (traces b and d), the echo shapes are strongly altered and the much faster decay is no longer exponential, as previously reported [5,10]. Since 180 • pulses only remove the linear phase dispersion introduced by the applied gradient, dipolar instabilities quickly take over and induce a very strong signal decay. At low G, minor refocussing is achieved by the 180 • pulses and the echo envelope is reminiscent of a typical dipolar-induced signal decay. At high G, the large gradient-induced frequency spread over the sample (∼ 10 × F dip ) dominates the very short measured FID lifetime. However, irreversible signal loss occurs during subsequent free evolution due to the development of unstable magnetisation patterns.

Dipolar-induced phase dispersion can be refocussed by continuous application of appropriately phased intense rf pulses that mix the longitudinal and transverse components of the magnetisation in such a way that it undergoes a time-reversed evolution, at half the pace of free evolution [11]. Successful implementation and use of the original magic sandwich pulse scheme in the context of liquid NMR has been recently reported [4]. We use a slightly modified MS cycle, 90 • y (360 • x ,-360 • x ) n 90 • y , that incorporates an additional 180 • rotation to refocus the gradient-induced linear phase dispersion. Periodic refocussing is achieved using a 'repeated magic sandwich' (RMS) sequence 90 • -τ -(MS -2τ ) n , where the duration of the MS cycle is 4τ [4,12]. The RMS sequence is thus analogous to a conventional spin echo sequence where the 180 • pulses have been replaced by MS cycles. It leads to the formation of so-called magic echoes with a time period equal to 6τ .

Figure 2 shows magic echo trains obtained under conditions similar to those of 2: Echo trains obtained using a RMS sequence with 6τ =145 ms in mixtures with X=2.4-3.0% at T =1.16 K, for various gradients and magnetisations. Upper plots: FID and first 3 echoes for G=0.89 µT/cm (left) and G=3.34 µT/cm (right), at low magnetisations (a: F dip =0.8 Hz, c: F dip =0.9 Hz) and at high magnetisations (b: F dip =10 Hz, d: F dip =13 Hz). Note that the low and high F dip traces in the upper-left-hand plot are nearly identical. No signal is acquired during the MS cycles. Lower plot: Semi-log plot of the echo amplitudes for these echo trains (filled/open symbols: low/high F dip ; circles/triangles: weak/strong G). For low F dip , decay rates Γ=0.099 and 0.49 s -1 for a and c, respectively. For high F dip , Γ=0.098 and 0.52 s -1 for b and d, respectively.

shapes and exhibit monoexponential decay at a rate that does not depend on F dip . NMR signals are efficiently recovered at all magnetisation densities since frequent time reversal heals spatial variations of magnetisation that develop during periods of free evolution.

For both CPMG and RMS echo trains, experimental decay rates Γ are obtained from an exponential fit to the decay of the squared amplitude of the Fourier transform of the band-pass filtered NMR signals. This method preserves the accuracy of the analysis down to low signal-to-noise ratios.

DISCUSSION

Our results confirm that CPMG sequences cannot be used to probe spin diffusion in highly magnetised liquids, due to the dramatic dipolar-induced signal losses (Fig. 1, traces b andd). In the stable linear regime (at low F dip ), the diffusion-induced decay rate for spin echo amplitudes is [13]:

Γ CPMG = D(2τ γG) 2 /12, (1)
where D is the spin diffusion coefficient in the mixture. The quadratic dependence of Γ CPMG on τ and G has been checked in a previous series of experiments, in which it was demonstrated that D can be accurately measured with this method in our apparatus [5]. Here, the two decay rates Γ extracted from traces a and c in Fig. 1 yield D = 1.54 × 10 -3 cm 2 /s at X = 2.8%, in good agreement with published values [5,14].

The RMS sequence appears to be a lot more robust than the CPMG sequence against dipolar effects. Over the range of experimental conditions we have explored, the shape and attenuation of magic echoes are insensitive to F dip . The expected diffusion-induced decay rate for spin echo amplitudes is [12]:

Γ RMS = 7 27 D(6τ γG) 2 /12, (2)
which is less than Γ CPMG because the dephasing action of G is arrested during each MS cycle [12]. Assuming that the decay rates measured for the strongest G only result from diffusion-induced attenuation (i.e., Γ = Γ RMS ), traces c and d in Fig. 2 combined with Eq. 2 yield D = 2.4 × 10 -3 cm 2 /s at X = 2.5%, which is 40% higher than expected. Decay rates at zero applied gradient (obtained from direct measurements at G=0 or from extrapolation of RMS data at various G) are also higher than expected from CPMG measurements (0.13-0.24 s -1 instead of 0.033 s -1 ). The physical origin of this additional damping has not yet been identified. We believe that rf heating of the sample (of order 1 mK/s during RMS sequences) has negligible influence on signal decay. In our experiment, evaporation hardly modifies the 3 He molar fraction (a 140 ppm/s relative decrease of X is estimated from thermodynamic data [15]). Additional signal losses incurred during RMS sequences, that are observed to be F dip -independent, may be caused by minor rf amplitude and timing imperfections that will soon be addressed. The present work clearly demonstrates that robust RMS sequences have the potential to enable one to perform accurate diffusion measurements in highly magnetised liquids. More generally, the efficient MS sequence provides a valuable time-reversal tool to investigate the complex non-linear NMR dynamics that can be encountered in polarised quantum fluids such as Bose-Einstein condensates, superfluid 3 He, degenerate 3 He-4 He mixtures, and 2D hydrogen gases.

FIG. 1 :

 1 FIG.1: Echo trains obtained using a CPMG sequence with 2τ =151 ms in mixtures with X=2.4-2.9% at T =1.16 K, for various gradients and magnetisations. Upper plots: FID and first 3 echoes for G=0.47 µT/cm (left) and G=3.34 µT/cm (right), at low magnetisations (a: F dip =0.4 Hz, c: F dip =0.5 Hz) or high magnetisations (b: F dip =11.5 Hz, d: F dip =25 Hz). 180 • rf pulses are indicated by vertical lines. Lower plot: Semi-log plot of echo amplitudes versus time, for these echo trains (filled/open symbols: low/high F dip ; circles/triangles: weak/strong G). For low F dip , decay rates Γ=0.054 and 1.42 s -1 are extracted from echo trains a and c, respectively.

Fig. 1 .

 1 FIG.2: Echo trains obtained using a RMS sequence with 6τ =145 ms in mixtures with X=2.4-3.0% at T =1.16 K, for various gradients and magnetisations. Upper plots: FID and first 3 echoes for G=0.89 µT/cm (left) and G=3.34 µT/cm (right), at low magnetisations (a: F dip =0.8 Hz, c: F dip =0.9 Hz) and at high magnetisations (b: F dip =10 Hz, d: F dip =13 Hz). Note that the low and high F dip traces in the upper-left-hand plot are nearly identical. No signal is acquired during the MS cycles. Lower plot: Semi-log plot of the echo amplitudes for these echo trains (filled/open symbols: low/high F dip ; circles/triangles: weak/strong G). For low F dip , decay rates Γ=0.099 and 0.49 s -1 for a and c, respectively. For high F dip , Γ=0.098 and 0.52 s -1 for b and d, respectively.