- Email: [email protected]

Contents lists available at ScienceDirect

Fuel journal homepage: www.elsevier.com/locate/fuel

Full Length Article

A multicomponent diffusion model for gas charges into oil reservoirs Shu Pan a,⇑, Julian Y. Zuo b, Kang Wang c, Yi Chen b, Oliver C. Mullins d a

DBR Technology Center, Schlumberger, 9450 17 Avenue, Edmonton, AB T6N 1M9, Canada Houston Formation Evaluation Center, Schlumberger, 150 Gillingham Lane, Sugar Land, TX 77478, USA c Beijing Geoscience Center, Schlumberger, Innovation Plaza, Tsinghua Science Park, Beijing 100084, China d Schlumberger-Doll Research Center, Schlumberger, Cambridge, MA 02139, USA b

a r t i c l e

i n f o

Article history: Received 25 December 2015 Received in revised form 30 March 2016 Accepted 11 April 2016 Available online 19 April 2016 Keywords: Multicomponent diffusion Gas charge Asphaltenes Density inversion

a b s t r a c t Asphaltene expulsion appears to be an important mechanism in the formation of a heavy oil layer and/or a tar mat at the base of an oil column during gas charges. To better understand this dynamic process, a multicomponent diffusion model has been developed to simulate the diffusion during gas charge processes in oil reservoirs. In this model, a complex moving boundary diffusion problem has been converted into a fixed boundary problem with simple boundary conditions. The constitutive equations linking the diffusion fluxes to the concentration gradient are given under the framework of the non-equilibrium thermodynamics. All the required thermodynamic properties are obtained from the Flory–Huggins regular solution model. This model has been used to study the gas charge processes in a representative ternary hydrocarbon reservoir system through two case studies. The studies have demonstrated the occurrences of asphaltene expulsion and two different types of density inversion during gas charges. A type I density inversion originates from a synergy of the asphaltene expulsion and an impermeable bottom boundary. It results in a concentration profile with a big, flat asphaltene hump near the reservoir bottom. The predominant cause of a type II density inversion is asphaltene expulsion, and it can induce buoyancy convection, which may accelerate the asphaltene migration to the base of the column. When operative, buoyancy convection is far more efficient than diffusion for asphaltene redistribution. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Oilfield reservoirs often have a layer of heavy oil and/or a tar mat at the base of an oil column. Heavy oil and tar are both characterized by elevated concentrations of asphaltenes, and the viscosity of crude oils and tar is roughly exponentially dependent on the asphaltene content. The impact of heavy oil and tar in oil reservoirs is significant. High viscosities limit the flow rate of crude oils. And tar can act as a pressure and flow seal of the aquifer preventing pressure maintenance of the reservoir by the aquifer during production. When sealing tar mats exist, wells for water injection into the reservoir must be placed within the oil column as opposed to in the aquifer, which is customary. This reduces overall oil recovery. In spite of the considerable importance of heavy oil and tar in reservoirs, the origin of this heavy oil and tar is often not understood or is even misunderstood. Without a proper understanding of the origin of heavy oil and tar, optimization of field development

⇑ Corresponding author. E-mail address: [email protected] (S. Pan). http://dx.doi.org/10.1016/j.fuel.2016.04.055 0016-2361/Ó 2016 Elsevier Ltd. All rights reserved.

planning is greatly hindered. In some work, the biodegradation was suggested as the origin of heavy oil and tar. Biodegradation of reservoir crude oils can result in consumption of certain oil components such as n-alkanes, thereby increasing the asphaltene fraction [1]. Nevertheless, biodegradation may not be the universal mechanism. The extent of biodegradation can be uncovered by analytical tools like gas chromatography due to predictable changes in oil composition. It has been revealed that many heavy oils and tars are not biodegraded. Moreover, the increase in asphaltene concentration in biodegradation is at most a factor of 3 [2], and is borne out in reservoir case studies [3]. However, reservoir field studies often show gradients of asphaltenes that are larger than a factor of 10 [4,5]. When biodegradation is not operative, for instance for reservoirs that are and have been too hot for microbial activity, the occurrence of heavy oil and tar are often considered to be from early oil, the first and least mature charge from the source rock. However, this explanation does not always apply. In some cases, there is no thermal maturity variation in oil columns with 10 variations in asphaltenes [6]. Moreover, the mechanism of early oil does not explain well how a solid tar (of neutral buoyancy) is expelled from nanodarcy source rock, then migrates kilometer

S. Pan et al. / Fuel 180 (2016) 384–395

after kilometer and finally becomes stuck in large fractures (or high porosity formations). In addition, early oil charges to the top of the reservoir, floating on water, with subsequent lighter charge displacing the heavier earliest charge pushing it down [7]. To account for tar in some reservoirs, this would require the current immovable tar at the base of the reservoir to be pushed down and laterally over 30 km without leaving any trace in the reservoir [4]. Clearly, other explanations are needed for asphaltene accumulations at the base of oil columns in reservoirs. One former impediment to the understanding of asphaltene accumulations at the base of oil columns had been the lack of any thermodynamic model to account for asphaltene gradients in reservoirs. The origin of this limitation had been the lack of understanding of asphaltene molecular and colloidal sizes. In recent years, the asphaltene molecular and colloidal structures in laboratory solvents and reservoir crude oils have been resolved and codified in the Yen–Mullins model [8,9]. With this nanostructure model of asphaltenes, a modified regular solution theory, the Flory–Huggins–Zuo equation of state (FHZ EOS) was developed and successfully accounts for equilibrated asphaltene gradients in reservoir crude oils [10,11]. At this juncture, time-dependent processes in reservoir crude oils are only partially understood from a first principles foundation. Nevertheless, oilfield studies guide thinking on dominant dynamic processes affecting fluid gradients. First-principles modeling must be developed to account for these observations. There is a simple and commonly occurring process in many reservoirs that can lead to asphaltene accumulation within a reservoir. In normal basin subsidence, the source rock becomes hotter and enters the oil window thereby filling reservoirs with oil. With further subsidence, the temperature increases yielding more extensive cracking of the hydrocarbons; thus, the source yields lighter hydrocarbons and there is subsequent reservoir filling by these lighter fluids. This normal basin subsidence process can thus result in a late gas charge into an oil reservoir. A chemistry axiom is ‘‘like dissolves like”; as a dark brown solid, asphaltenes are not chemically similar to (colorless) gas; asphaltenes and gas are not mutually soluble. Consequently, a late gas charge into an oil reservoir can yield an increase of solution gas with

385

expulsion of asphaltenes from the oil [12]. A well-known case of asphaltene expulsion with increased solution gas is shown in Fig. 1 [4,13]. Fig. 1 shows data corresponding to a late gas charge into an oil reservoir. The reservoir filling process occurred without mixing, as expected [7]. Methane then diffused down into the oil column increasing solution gas up to 8000 scf/bbl. This high solution of gas has led to an expulsion of asphaltenes. These asphaltenes migrated to the base of the reservoir, as shown in Fig. 1b. Fig. 2 shows two time steps of this sequence of gas diffusion into an oil column from a late gas charge. Adjacent fault blocks depicted in Fig. 2 are in different stages of the same process, a late gas charge into an oil reservoir. Fig. 2a shows a cartoon of the reservoir immediately after charge. Fig. 2b shows the gas diffusion into the oil and the migration of asphaltenes down in the oil column. Both the gas/oil ratio (GOR) gradient and the asphaltene gradient are not equilibrated due to this flux. Fig. 2c shows the adjacent fault block where the solution gas and asphaltene migration reached the base of the reservoir. Asphaltene content of core extracts show the asphaltenes formed a tar mat at the base of the oil column (Fig. 2c). In one reservoir, the formation of heavy oil and asphaltenes from asphaltene instability at the crest is up to 30 km away from the crest [14]. Thus, the accumulation process necessarily may involve convection because diffusion over these distances would take a trillion years. Many studies have been done to confirm the expected uniform asphaltene chemistry throughout the reservoir expected for convective currents [15]. Moreover, other studies have ruled out other possible mechanisms to account for the massive asphaltene accumulation in this reservoir [6]. The new understanding of timedependent processes of reservoir fluids in geologic time is being treated within the new discipline ‘‘reservoir fluid geodynamics” [4,16] and is enabled by the new asphaltene thermodynamics of the FHZ EOS and the Yen–Mullins model along with the advent of downhole fluid analysis [17] to characterize fluid gradients. The question remains how the asphaltenes can migrate great distances in reservoirs. This migration of asphaltenes in reservoirs must be treated using a first-principles foundation. Clearly, asphaltenes flocs cannot migrate through reservoirs. Even

Fig. 1. Oil properties variation height in the oil column for a single reservoir, (a) 24 bottles of dead crude oils (with solution gas flashed) from a single reservoir stacked by depth, (b) the asphaltene distribution measured downhole and in the laboratory, and (c) gas–oil ratio (GOR) in scf/bbl and the saturation pressure gradient.

386

S. Pan et al. / Fuel 180 (2016) 384–395

Fig. 2. A late gas charge into an oil reservoir. (a) The fluid distribution immediately after the late gas charge; (b) large, disequilibrium gradients of both GOR and asphaltenes associated with gas diffusion into the oil column and expulsion of asphaltenes, and (c) the end of this gas diffusion and asphaltene downward migration with formation of a tar mat at the base of the oil column (high asphaltene content in core extracts).

asphaltenes flocs as small as 1 lm fall until landing on a surface; thermal energy (at reservoir temperatures) is much too small to lift these flocs [18]. However, the addition of gas to oil does not immediately flocculate asphaltenes. Even when slightly destabilized, asphaltenes can remain nanocolloidally suspended in crude oils [8]. Consequently, with gas diffusion from above, asphaltene nanocolloidal species can diffuse downward away from the diffusive gas front [5,6,14]. In this manner, asphaltenes can accumulate at the base of the diffusive gas front. The enrichment of asphaltenes resulting from asphaltene expulsion for regions of high solution gas can yield a hump in vertical density profile [6], referred to as density inversion. Even small magnitudes of density inversions may yield convection in reservoir crude oils that are fast in geologic time [6]. Previous modeling of this asphaltene diffusion downward with advancing solution gas downward showed that indeed asphaltenes can pile up ahead of the solution gas front [6]. However, that initial rudimentary modeling did not account for the increased liquid phase volume with increasing methane addition via diffusion. Thus, this initial modeling did not provide an accurate density of the reservoir fluids with commencement of gas diffusion. Here, a multicomponent model is developed to account for gas diffusion into oil and for asphaltene expulsion from regions of high solution gas. From the modeling perspective, a gas charge into an oil column can be described as an unsteady state multicomponent diffusion process with moving boundaries. Three elements are crucial to model such a process: a reliable thermodynamic representation of asphaltenes contained in hydrocarbon fluid, physically sound description of multicomponent diffusion, and appropriate mathematical and numerical techniques to solve an unsteady state moving boundary diffusion problem. These key elements are briefly reviewed as follows. Asphaltenes are defined as the heaviest group of components in crude oil that are soluble in aromatics (toluene) and insoluble in normal alkanes such as n-heptane. It is difficult to deal with phase behavior for asphaltene-containing mixtures using cubic EOS because of many unknown asphaltenes properties. As mentioned previously, the Flory–Huggins–Zuo equation of state (FHZ EOS) has successfully been developed [10,11,13,19,20] to account for asphaltene equilibrium distributions and delineate reservoir connectivity in oil reservoirs coupled with the Yen–Mullins model and downhole fluid analysis for light oil, black oil, and heavy oil

[21]. If the gravity term is ignored, the FHZ EOS is reduced to the Flory–Huggins regular solution model, which has been widely used for asphaltene phase separation calculations [22–24]. Therefore, the Flory–Huggins regular solution model is used to describe non-ideality. Over the years, different types of diffusion theories have been published to describe the diffusion in multicomponent systems. Two of the most important multicomponent diffusion models are the generalized Fick’s diffusion and the generalized Maxwell– Stefan theory. As discussed in [25,26], both formulations have been generalized and to some extend redefined from their original forms under the framework of Onsager’s non-equilibrium thermodynamics theory, which makes these two formulations equivalent to each other. The generalized Fick’s law explicitly accounts for the linear relationship between diffusion fluxes and the corresponding concentration gradients by introducing a diffusion coefficient matrix that is dependent on reference frames. This formulation facilitates the direct measurement of the Fick’s diffusion coefficients. On the other hand, the generalized Maxwell–Stefan theory uses velocity differences and therefore the generalized Maxwell–Stefan diffusion coefficients do not rely on the selection of reference frames. It was found the diffusion coefficients from the generalized Maxwell–Stefan formulations appear less concentration dependent than those from the generalized Fick’s formulation. Due to this fact, most of modern engineering modeling for multicomponent diffusion in liquid mixtures start from an estimation of the generalized Maxwell–Stefan diffusion coefficients, which are then converted to appropriate formulations in their applications. This approach is also adopted in the present work. The analytical or numerical approaches for unsteady state moving boundary diffusion problems have been extensively reviewed and discussed by Crank in his monograph [27]. However, many of these methods do not fit well to the present study. Some of them assume constant coefficients to simplify and linearize the problem. The others are only suitable for binary diffusion problems and ambiguous for multicomponent systems where inappropriate selection of reference may result in complicated boundary conditions. Following a suggestion from Crank, Duda and Vrentas [28] have reformulated a 1D transient diffusion problem with a moving boundary into a simple form of a fixed boundary problem by utilizing different length and concentration variables. The fixed boundary problem can then be solved by any appropriate

387

S. Pan et al. / Fuel 180 (2016) 384–395

analytical or numerical method. Duda and Vrentas’ formulation has been successfully applied to the study of swelling and diffusion-induced convection in polymer sorption processes [29]. Although originally suggested for binary mixtures, this formulation has a generalized form to handle variable coefficients and multicomponent systems. In the present work, a theoretical multicomponent diffusion model is proposed to simulate gas charge processes in an oil reservoir and to investigate the density inversion induced by the gas charge. In next section, the problem will be described in detail, along with the formulation of a 1D unsteady state model for such a problem. The problem statement and model formulation are followed by the results and discussion from two cases studies with two types of density inversion, named type I and II, in the simulation results. After that, conclusions and recommendations are given.

2. Problem statement and model formulation To capture the most essential physical picture, the present study has been focused on an idealized gas charge to a homogeneous oil column. Following Stainforth’s petroleum reservoir filling theory [7], the process can be described as follows. Upon the arrival at the reservoir from the source rock, the gas charge quickly flows upward to the top of an existing homogeneous reservoir oil column and forms a gas cap at the crest. The presumption is that this gas charge is rapid and occurs before any significant mass exchange with the oil. One can assume the compositions of hydrocarbons from either the oil column or the charge gas remain unchanged at this starting point of time. The gas–oil contact (GOC) is the only location where a thermodynamic equilibrium compositional gradient is held. The nonequilibrium distribution of chemical potential drives diffusional migrations in the hydrocarbons, which leads to mass exchange between the gas cap and oil column as well as the redistribution of species along the oil column. At the GOC, the primary diffusional migration is the diffusion of gas components into the oil column. It is accompanied by the diffusion of a small amount of intermediate light components from the oil column to the gas cap. As a result, the gas cap is gradually absorbed by the swelling oil column. The GOC moves upward toward the top seal until either a global thermodynamic equilibrium is achieved or the gas cap is fully absorbed by the oil column. During the absorption of the gas cap, the aquifer at the bottom of the hydrocarbon column maintains a static pressure, thus a mechanical equilibrium in the system, and the location of oil–water contact (OWC) adjusts with the change of volume in the hydrocarbon column (oil column + gas cap). The volume and location change in this absorption process is illustrated in Fig. 3a. A closed hydrocarbon system is assumed in this 1D illustration. A fixed reservoir pressure and temperature are presumed. The only fixed reference frame is the top seal rock. Given a constant temperature and pressure in the system, the hydrocarbons migrate in two major mechanisms: diffusion driven by chemical potential gradients and convection induced by the local density variations due to the redistribution of local species concentration. With the introduction of some assumptions, the aforementioned process can be further simplified to a system shown in Fig. 3b. In this system, only an oil column is taken into account, and the bottom of oil column (OWC) is a fixed reference boundary impermeable to all three components. These assumptions are given as below: The pressure and temperature remain constant in the system. At the GOC, mass transfer only occurs from gas cap to the oil column, whereas the mass transfer from the oil column to the gas cap is neglected.

Buoyancy convections are absent in the system. In the real world, buoyancy convections are induced by density inversions. This effect is neglected in the present study. Indeed, the density inversions found herein can be used in future work to model convection. After all the above assumptions and simplifications, a multicomponent diffusion model has been developed for the gas charge processes to an oil column. It is assumed that there are N hydrocarbon components in the system. The components 1 to N c are present in both the charge gas and the oil column, whereas the components N c þ 1 to N exist only in the oil column. The derivation of the model starts from the mass conservation equation of each component i, as given in Eq. (1):

@ðci Þ @ðci ui Þ þ ¼ 0 ði ¼ 1; . . . ; NÞ @t @z

ð1Þ

Here ci is the molarity of component i; ui is the migration velocity of the component i; t is time, and z is the space coordinate originated from the bottom of the oil column. The initial condition can be expressed by

ci ¼ ci0 ;

t ¼ 0 ði ¼ 1; . . . ; NÞ

ð2Þ

where ci0 is the initial molarity of component i in the existing oil column, a constant for a homogeneous oil column. At the lower boundary, i.e., the bottom of the oil column, nonpermeable boundary conditions are applied to all components by

ui ¼ 0:0;

@ z ¼ 0 ði ¼ 1; . . . ; NÞ

ð3Þ

The upper boundary, i.e., the GOC at top of the oil column, is a moving boundary. It is permeable for the gas components in the gas charge, whereas it is impermeable for the other components. The components in the gas charge are at the saturation concentrations at the upper boundary. Assuming the location of this moving upper boundary is Z T ðtÞ, one could specify the corresponding boundary conditions for each component as

ci ðt; ZðtÞÞ ¼ csat i ðT; P; c 1 ; . . . ; c N Þ;

@ z ¼ Z T ðtÞ ði ¼ 1; . . . ; Nc Þ ð4aÞ

ui ðt; ZðtÞÞ ¼

dZðtÞ ; dt

@ z ¼ Z T ðtÞ ði ¼ Nc þ 1; . . . ; NÞ

ð4bÞ

As suggested by Duda and Vrentas [28], Eqs. (1)–(4) can be reformulated into a simple form of a fixed boundary problem. The reformulation is achieved by two key steps. The first step is to express the diffusion by the molar diffusion fluxes with respect to the velocity of a reference components which only exist in the oil column. Here the last component N is specified as the reference component. The definitions of diffusion fluxes are given below:

J i ¼ ci ðui uN Þ ði ¼ 1; . . . ; N 1Þ

ð5aÞ

JN ¼ 0

ð5bÞ

Eq. (1) can then be rewritten as

@ðci Þ @ðJ þ ci uN Þ ¼ i @t @z @ðcN Þ @ðcN uN Þ ¼ @t @z

ði ¼ 1; . . . ; N 1Þ

ð6aÞ ð6bÞ

The second step is to apply a coordinate system transformation technique by introducing a new space variable n:

Z

n¼

z

cN ðt; zÞdz

ð7Þ

0

For any given function of t and z; f ðt; zÞ, one can obtain the following relationship:

388

S. Pan et al. / Fuel 180 (2016) 384–395

Topseal: fixed reference GOC

GOC z

OWC

z=0.0

OWC: fixed reference

(a)

(b) Fig. 3. Schematic illustration of hydrocarbon column change due to diffusional migration.

@f ðt; nÞ @f ðt; zÞ @nðt; zÞ 1 @f ðt; zÞ ¼ ¼ @n @z @z cN @z t t t t

ð8aÞ

In conjunction with Eqs. (1), (3), and (4) for component N, another important relationship can be derived as follows:

@f ðt; nÞ @f ðt; zÞ @f ðt; zÞ @zðt; nÞ ¼ þ @t @t @z @t n z t n @nðt;zÞ @t @f ðt; zÞ @f ðt; zÞ z ¼ @nðt;zÞ @t @z z t @z

Z

ð8bÞ

ð9bÞ

With the relationship given in Eqs. (9a) and (9b), one can transform Eqs. (6) from equations of t and z to those of t and n.

ð10aÞ

ð10bÞ

Eq: ð10aÞ cc2i Eq: ð10bÞ, and N

applying some new concentration variables Si ¼ ccNi , one would obtain

@ t ¼ 0 ði ¼ 1; . . . ; N 1Þ

@ n ¼ 0 ði ¼ 1; . . . ; N 1Þ

Z T ðtÞ

cN ðt; zÞdz ¼ cN0 L

ð12Þ

ð13Þ

ð14Þ

Thus, at the upper boundary, we have n ¼ cN0 L ¼ const. The upper boundary now becomes a fixed boundary with boundary conditions:

Si ¼ csat i =c N ;

@f ðt; zÞ @f ðt; nÞ @f ðt; nÞ ¼ c N uN @t @t @n z n t

1 cN

ci0 ; cN0

0

ð9aÞ

Applying the operation

Si ¼

The most important change is applied to the upper boundary. It is converted from a moving boundary to a fixed boundary due to the fact that total mass of component N stays constant in the system. Given the original height of oil column, L, i.e., Z T ð0Þ ¼ L, the total moles of component N at any time t is given as

Eqs. (8a) and (8b) can be rewritten as

@cN @cN @ðcN uN Þ @uN ¼ c N uN cN ¼ c2N @n @t n @n t @n t t

ð11Þ

The initial conditions now become

J i ¼ 0:0;

@f ðt; zÞ @f ðt; nÞ ¼ cN @z @n t t

@ci @ci @ðJ i þ ci uN Þ ¼ cN uN cN @n @t n @n t t @J i @uN ¼ cN c N ci ði ¼ 1; . . . ; N 1Þ @n t @n t

ði ¼ 1; . . . ; N 1Þ

The lower boundary condition can be reduced to

t

R z @cN dz z @f ðt; zÞ @f ðt; zÞ 0 @t ¼ @t @z cN ðt; zÞ z t R z @ðcN uN Þ dz 0 @z @f ðt; zÞ @f ðt; zÞ z ¼ þ @t @z cN ðt; zÞ z t @f ðt; zÞ @f ðt; zÞ þ uN ¼ @t @z z t

@Si @J i ¼ @t n @n t

J i ¼ 0:0;

@ n ¼ cN0 L ði ¼ 1; . . . ; Nc Þ

ð15aÞ

@ n ¼ cN0 L ði ¼ Nc þ 1; . . . ; N 1Þ

ð15bÞ

To this point, the model has been formulated into a concise and tangible partial differential equation (PDE) set in the expression of Eqs. (11)–(13), and (15). The next step is to introduce appropriate constitutive equations to express the fluxes J i in terms of variables Si , and eventually to obtain a closed set of PDEs. The required constitutive equations can be derived from the generalized Maxwell–Stefan theory, which links the chemical potential gradients to relative velocities between components as below: N 1 X xi xi xk ðui uk Þ xi xN ðui uN Þ r z li ¼ Ðik ÐiN RT k¼1

ði ¼ 1; . . . ; N 1Þ ð16Þ

where li is the chemical potential of component i; xi is the mole fraction of component i; R is the ideal gas constant, and Ðik is the generalized Maxwell–Stefan diffusion coefficients between the pair of components i and k.

389

S. Pan et al. / Fuel 180 (2016) 384–395

In conjugation with the definition of J i and Si , Eq. (16) can be expressed as

r z li

0 1 J J Ji N1 Sk i k Si Sk RT @X Si A ¼ þ ct k¼1 ÐiN Ðik 0 1 N1 N RT BX 1 1X Sk C ¼ J JA @ ct k¼1; Ðik k Si k¼1 Ðik i

ð17Þ

k–i

where SN ¼

cN cN

¼ 1 and ct ¼

PN

N X Sk 6 S11 Ð1k 6 k¼1 6 6 7 RT 6 .. .. 6 6 7¼ . 4 5 . ct 6 6 6 rz lN1 4 1

3

ÐN1j

..

1 Ð1N1

.

1 SN1

.. . N X k¼1

Sk ÐN1k

72 3 2 3 7 J1 J1 7 76 . 7 RT 6 . 7 76 . 7 ¼ 7 ½B6 74 . 5 4 .. 5 ct 7 7 J J N1 5 N1

ð18Þ Here the elements for the matrix ½B are given as

½Bij ¼

1 Ðij

when i – j

ð19aÞ

and

½Bii ¼

N 1X Sk Si k¼1 Ðik

ð19bÞ

The gradients of chemical potential in Eq. (18) are expressed in z-space. They can be converted to n-space using Eq. (9a) and then Eq. (18) can be reorganized as

2

3

2

3

2

3

rz l1 rn l1 J1 ct 6 .. 7 6 7 ct cN 1 6 7 .. .. ¼ ½B1 4 ½B 4 . 5¼ 4 5 5 . . RT RT J N1 rz lN1 rn lN1 2 3 P rn S1 c2N Nk¼1 Sk 1 @ l 6 7 .. ¼ ½B 4 5 . RT @S rn SN1

In the present work, each component is assumed to be of constant partial molar volume, v i ; therefore, cN can be calculated by

c N ¼ PN

i¼1

v i Si

ð21Þ

Meanwhile, the Flory–Huggins regular solution model is adopted to describe thermodynamics in the system. The chemical potential can be expressed as [30]

(

xi v i v i v i XX h / / ðdi dj Þ2 þ1 þ V V RT j k j k ) i þ 2lij di dj 0:5ðdj dk Þ2 ljk dj dk

ð23Þ

Z Z T ðtÞ ¼

cN0 L

0

1 dn cN ðt; nÞ

ð24Þ

In the present work, a simplified but representative ternary hydrocarbon system has been studied using the model described above. In this system, all hydrocarbons have been lumped into three pseudocomponents—gas, maltenes and asphaltenes. The gas component primarily consists of methane and represents the species in the gas cap. The gas is assumed to be the only pseudocomponent present in the gas cap. The asphaltenes are simply the asphaltenes contained in hydrocarbons. The maltenes represent all intermediates between the gas and the asphaltenes. For this ternary system, Eqs. (11)–(13), (15), and (20)–(22) can be written as below, where subscript G stands for the component gas, A stands for the component asphaltenes, and M stands for the component maltenes. The PDE set is

@SG @J G ¼ @t n @n t

ð25aÞ

@SA @J A ¼ @t n @n t

ð25bÞ

The initial conditions are

ð20Þ

In Eq. (20), li and cN are functions of Si : They can be determined by any appropriate relationships and thermodynamic models. And @l is the Jacobian matrix of chemical potentials li @S (i ¼ 1; . . . ; N 1Þ with respect to Si (i ¼ 1; . . . ; N 1Þ; i.e., @@Sl ij ¼ @@Slji .

1

1 dn cN ðt; nÞ

As a special case of Eq. (23), the location of the moving upper boundary can be calculated by

3

n

0

i¼1 c i .

2

rz l1

Z z¼

Eq. (17) can be rewritten in a matrix format as 2

The equations above can be solved numerically using standard nonlinear diffusion equation solvers. In this work, PDEPE from Matlab 2012b [31] was used for solutions. It should be noted that the solutions to the equation set are given the concentration variables Si in terms of t and n. One additional relationship must be engaged to convert the solutions from n space back to z space, and that is the space coordinate in the real world:

li ¼ l0i þ RT ln

SG ¼

cG0 ; cM0

@t¼0

ð26aÞ

SA ¼

cA0 ; cM0

@t¼0

ð26bÞ

The boundary conditions are

J G ¼ 0:0;

@n¼0

ð27aÞ

J A ¼ 0:0;

@n¼0

ð27bÞ

SG ¼ kG ;

@ n ¼ cM0 L

ð28aÞ

J A ¼ 0:0;

@ n ¼ cM0 L

ð28bÞ

where kG ¼ csat G =cN , is the solubility constant of gas in maltenes. Meanwhile, the diffusion fluxes are given in terms of Si as

ð22Þ

JG ÐGM SG ðÐAM SG þ ÐGA Þ ÐAM ÐGM SA SG c2 1 þ SA þ SG ¼ M RT ÐAM SG þ ÐGA þ ÐGM SA JA ÐAM ÐGM SA SG ÐAM SA ðÐGA þ ÐGM SA Þ 2 3 @ lG @ lG

r S @S @S n G 4 @ lG @lA 5 ð29Þ A A rn S A @S @S G

where /i is the volume fraction of component i; di is the solubility parameter of component i; lij is the binary interaction parameters between components i and j, and V is the mixture molar volume P which is calculated by V ¼ i xi v i . At this point, the full set of equations for the proposed model have been laid out, including Eqs. (11)–(13), (15), and (20)–(22).

A

where

cM ¼

1

v M þ v A SA þ v G SG

Chemical potential is calculated by Eq. (22).

ð30Þ

390

S. Pan et al. / Fuel 180 (2016) 384–395 Table 2 Generalized Maxwell–Stefan diffusion coefficients (1011 m2/s).

3. Simulation results and discussion In this section, some simulation results from the proposed multicomponent diffusion model will be presented. Since the purpose of the present work is to provide a framework to analyze coupled diffusion of gas and asphaltenes, rather than conduct a comprehensive sensitivity analysis, the simulations were performed for a typical hydrocarbon reservoir system with fluid and reservoir parameters given in Section 3.1. Some of the thermodynamic model parameters were adjusted within reasonable ranges to appreciate the diverse and complex nature of reservoir fluid. These adjustments resulted in various simulation results, which will be shown and discussed in Sections 3.2 and 3.3. Further model sensitivity analysis will appear in future papers.

3.1. Fluid and reservoir parameters for simulations In this study, the formation is assumed to be sandstone of porosity 10% and permeability 300 md. This reservoir contains a typical intermediate black oil with viscosity 1 cP. The thermodynamic parameters of each lumped component—gas, asphaltenes, and maltenes—have been estimated from the original fluid composition and given in Table 1. Three generalized Maxwell–Stefan diffusion coefficients, ÐGM ; ÐAM and ÐGA , are required to model the diffusion in such a ternary mixture. To the best knowledge of authors, these parameters have not been accurately determined, particularly for ÐAM and ÐGA . Some estimations have been made in this work, and the estimated generalized Maxwell–Stefan diffusion coefficients are given in Table 2. The estimations were made following two general guidelines: With the tortuosity in porous media, diffusion coefficients in typical sandstones is 10 times lower than those in free volume [32]. Since the average viscosity and molecular size of three pairs of components in the concentration range of the present study can be generally ranked as gas–maltenes < gas–asphaltenes < maltenes–asphaltenes, it is fair to assume that ÐGM >GA > ÐAM . As such, some typical values of ÐGM ; ÐAM , and ÐGA are given in Table 2. Future work will explore the impact of the ratios and values of these diffusion constants. The initial and boundary conditions of the simulations are summarized in Table 3. The initial height of hydrocarbon L is assumed to be 80 m. The initial molar concentrations in the oil column are 3

3

cG0 ¼ 1911:09 mole=m ; cA0 ¼ 13:38 mole=m ,

and

3

cM0 ¼ 3822:18 mole=m . At the upper boundary, the mole ratio between gas and maltenes is fixed to kG ¼ 1:2 mole=mole. The only missing model parameters are lij , the binary interaction parameters in the Flory–Huggins regular solution model. lij are introduced as the correction to the geometric mean rule for the cohesive energy density. Unfortunately, there has been little success in the efforts to correlate lij even in binary mixtures. Prausnitz et al. [30] have studied 23 pairs of aromatic-saturated

ÐGM

ÐAM

ÐGA

1.6

0.5

0.9

hydrocarbon mixtures and found those lij fall into a range of 0.03 to 0.015. In the following two case studies, the binary interaction parameters lij have been adjusted to reflect the diversity and complexity of hydrocarbon composition, which has been heavily lumped into three components in the present work. 3.2. Case study 1: Gas charge process with type I density inversion Table 4 lists lij for the case study 1. These parameters were adopted with other parameters given in Section 3.1 to simulate the gas charge process in a 1-million-year (MY) time span. In our simplified gas charge process, the component gas continuously enters the oil column while the other two components remain in the oil column. The oil column swells with the addition of new mass from the charge gas. The height of oil column, i.e., the distance between the GOC and OWC, increases with continuing gas charge. As shown in Fig. 4, the height of oil column grows rapidly at the beginning due to higher diffusion driving force. The growth becomes relatively stable with time. The oil column height rises about 5 m after 3 MY. Meanwhile, Fig. 5 shows the change of total mass of gas in the oil column follows the same growth trend as the height of oil column growth due to the assumption of constant partial molar volume. The distribution of mass concentration of gas, asphaltenes, and maltenes are plotted against the normalized height in Figs. 6–8, respectively. Here the normalized height, defined as z=Z T ðtÞ, is introduced so that all profiles can be plotted under the same scale. Fig. 6 shows the concentration of gas remains almost unchanged at the top of oil column and decreases with the depth. With the gas charge going on, more and more gas migrates to the lower part of oil column, and the gas mass concentration increases accordingly. The profile change of asphaltene mass concentration is plotted in Fig. 7. The asphaltene expulsion phenomenon is clearly demonstrated in this figure. Corresponding to the frontal edge of gas migration, a hump is formed in the asphaltene mass concentration profile. The hump is pushed downward along the direction of gas migration due to the chemical incompatibility between the gas and asphaltenes. The hump is slightly spread out when the gas distributes more uniformly in the first 2 MY. Then it gradually builds up near the bottom while the downward migration of asphaltenes pushed by gas migration is stopped by the impermeable bottom boundary. This synergy of asphaltene expulsion and the impermeable bottom boundary leads to a concentration profile with a big, flat asphaltene hump near the reservoir bottom at 3.5 MY, a shape similar to the field observation in Fig. 2c. This similarity suggests the asphaltene expulsion might be relevant to account for observed large asphaltene gradients in the field studies. The simulation could not continue after 3.5 MY as the oil with peak asphaltene concentration reaches the

Table 1 Component thermodynamic parameters. Component

Molecular or nanoaggregate weight (g/mole)

Molar volume (cm3/mole)

Solubility parameter (MPa0.5)

Solubility in maltenes (mole/mole)

Gas Asphaltenes (nanoaggregate) Maltenes

16 3000 191

51.00 2500.00 227.38

8.346 20.895 17.475

1.2 – –

391

S. Pan et al. / Fuel 180 (2016) 384–395 Table 3 Initial and boundary conditions.

The Gas mass concentration profile 1

G0 SG ¼ ccM0 ¼ 0:5; @ t ¼ 0

Initial conditions

A0 ¼ 0:0035; @ t ¼ 0 SA ¼ ccM0

0.9

Table 4 lij for case study 1.

0.8

Normalized height

JG ¼ 0:0; @ n ¼ 0 JA ¼ 0:0; @ n ¼ 0 SG ¼ kG ¼ 1:2 mole=mole; @ n ¼ cM0 L ¼ 3822:18 mole=m3 80 m JA ¼ 0:0; @ n ¼ cM0 L ¼ 3822:18 mole=m3 80 m

Boundary conditions

0.7 t=0.7 MY t=1.4 MY t=2.1 MY t=2.8 MY t=3.5 MY

0.6 0.5 0.4 0.3

lGM

lAM

lGA

0.2

0.0

0.01

0.0

0.1 0 30

35

40

45

50

55

60

65

70

3

Oil column height change during gas charge

Gas mass concentration, kg/m

85 Fig. 6. Case 1, gas mass concentration profile in the oil column.

The Asphaltenes mass concentration profile 1

83

0.9 82

0.8

Normalized height

Oil column height, m

84

81

80

79 0

0.5

1

1.5

2

2.5

3

0.7 0.6 0.5 0.4 0.3

t=0.7 MY t=1.4 MY t=2.1 MY t=2.8 MY t=3.5 MY

3.5 0.2

Time, MY

0.1 Fig. 4. Case 1, oil column height change during gas charge.

0 0

Total gas mass in oil column

20

30

40

50

60

70

Asphaltenes mass concentration, kg/m

4000

80

3

Fig. 7. Case 1, asphaltene mass concentration profile in the oil column.

3800

Total gas mass, kg

10

3600 3400 3200 3000 2800 2600 2400 0

0.5

1

1.5

2

2.5

3

3.5

Time, MY Fig. 5. Case 1, total mass of gas in the oil column.

two-phase region and cannot be dealt with using the current model. Meanwhile, as shown in Fig. 8, depressions form in the maltene mass concentration profile near where those humps in asphaltene profiles are observed.

The overall mass density distribution in the oil column is determined by the summation of local mass concentrations of all three components. Fig. 9a plots the whole density profile of this case study, and the top of density profile has been enlarged and given in Fig. 9b. It is shown that the density increases from the top to the bottom of oil column when the frontal edge of gas migration is far from the bottom of the reservoir. The density variation along the oil column during this period is primarily controlled by the gas migration while the hump in asphaltenes mass concentration has been well compensated by the depression in maltene mass concentration profile. After 2.8 MY, the gas migration front approaches the bottom of the reservoir. The depression in maltene mass concentration profile appears insufficient to compensate the growth of hump in the asphaltene concentration profile. And therefore a density inversion starts to form. This density inversion due to the synergy of both asphaltene expulsion and the impermeable bottom boundary is named a type I density inversion. The convection induced by a type I density inversion may have limited impact on the distribution of species in the reservoir fluid since the inversion only forms near the bottom of reservoir and also the fluid in

392

S. Pan et al. / Fuel 180 (2016) 384–395 Table 5 lij for case study 2.

The Maltenes mass concentration profile 1 0.9

lGM

lAM

lGA

0.015

0.03

0.015

Normalized height

0.8 0.7 t=0.7 MY t=1.4 MY t=2.1 MY t=2.8 MY t=3.5 MY

0.6 0.5 0.4

this region has significantly increased viscosity due to the accumulation of asphaltenes. Therefore, the impact of convection induced by type I density inversion is not further discussed in the present paper.

0.3

3.3. Case study 2: Gas charge process with density inversion

0.2 0.1 0 660

670

680

690

700

710

720

730

Maltenes mass concentration, kg/m3 Fig. 8. Case 1, maltene mass concentration profile in the oil column.

The density profile 1 0.9

(a)

Normalized height

0.8 0.7 t=0.7 MY t=1.4 MY t=2.1 MY t=2.8 MY t=3.5 MY

0.6 0.5 0.4 0.3 0.2 0.1 0 720

740

760

780

Density ρ, kg/m

800

820

3

The enlarged density profile 1 0.9

(b)

Normalized height

0.8 0.7 t=0.7 MY t=1.4 MY t=2.1 MY t=2.8 MY t=3.5 MY

0.6 0.5 0.4 0.3 0.2 0.1 0 795

800

805

Density ρ, kg/m3 Fig. 9. Case 1, density profile in the oil column, (a) the whole profile and (b) enlarged profile with density inversion.

In case study 2, the binary interaction parameters lij were adjusted within the range suggested by Prausnitz et al. [30] as given in Table 5. In reality, this variation of lij could reflect the different characteristics, such as aromaticity and composition distributions, of heavily lumped asphaltene and maltene components. Similar to case study 1, lij from Table 5 and other model parameters given in Section 3.1 have been adopted to simulate the gas charge process for 1 MY. The growth of oil column height and change of total mass of gas are plotted in Figs. 10 and 11. The minor changes to lij have shown only a little impact on the absorption of gas and the curves appears nearly the same as those from the case study 1. Figs. 12–14 plot the distributions of mass concentration of gas, asphaltenes, and maltenes, respectively. The asphaltene expulsion phenomenon is again illustrated in these figures. A significant hump forms in asphaltene mass concentration profile and spreads out gradually while being driven downward by the migration of gas. Compared to the results from the case study 1, the height of humps in asphaltene mass concentration profile has slightly dropped from above 65 kg/m3 to nearly 60 kg/m3; however, the depressions in the maltene mass concentration profile almost disappears. These changes in mass concentration profiles result in a hump in density profile along the oil column at the early stage of charge, well before the migration approaches the bottom of the reservoir, as shown in Fig. 15. This density inversion near the top of the oil column appears to be predominantly due to the asphaltene expulsion with negligible impact from the impermeable bottom boundary. It is therefore named a type II density inversion, to distinguish it from the type I density inversion, which only forms near the bottom of the reservoir due to the synergy of both asphaltene expulsion and impermeable bottom boundary. The overall density profile is provided in Fig. 15a, and the density inversion part is enlarged and given in Fig. 15b. Fig. 15 also shows that the density hump moves downward and spreads with the ongoing gas charging. It appears similar to the migration behavior of the asphaltene mass concentration profile since the latter provides the major contribution to the formation of the density hump. Unlike a type I density inversion, the convection induced by a type II inversion may have a significant impact on the hydrocarbon species distribution in the reservoir, since type II density inversion forms near the top of reservoir in the early stage of gas charge. A further analysis is therefore given below to gain a rough understanding of this impact. Type II density inversion here can be characterized by the difference between maximum density qmax in the oil column and the density at the bottom of oil column qz¼0 :

Dq ¼ qmax qz¼0

ð31Þ

Dq is equivalent to zero without the presence of density inversion. The evolution of Dq in the current simulation is described in a

393

S. Pan et al. / Fuel 180 (2016) 384–395

Oil column height change during gas charge

The Asphaltenes mass concentration profile 1

83

0.9 0.8

Normalized height

Oil column height, m

82.5

82

81.5

81

0.7 t=0.2 MY t=0.4 MY t=0.6 MY t=0.8 MY t=1.0 MY

0.6 0.5 0.4 0.3 0.2

80.5

0.1 80 0

0.2

0.4

0.6

0.8

0

1

0

10

20

30

40

50

60

70

Asphaltenes mass concentration, kg/m3

Time, MY Fig. 10. Case 2, oil column height change during gas charge.

Fig. 13. Case 2, asphaltene mass concentration profile in the oil column.

Total gas charge The Maltenes mass concentration profile

3200

1

3100

0.9 0.8

3000

Normalized height

Total gas charge, kg

3300

2900 2800 2700 2600 2500

0.7 0.6 0.5 0.4

t=0.2 MY t=0.4 MY t=0.6 MY t=0.8 MY t=1.0 MY

0.3 0.2

2400 0

0.2

0.4

0.6

0.8

1

0.1 0 660

Time, MY

670

680

690

700

710

720

730

740

3

Fig. 11. Case 2, total mass of gas in the oil column.

Maltenes mass concentration, kg/m

Fig. 14. Case 2, maltene mass concentration profile in the oil column.

The Gas mass concentration profile 1 0.9

Normalized height

0.8 0.7 t=0.2 MY t=0.4 MY t=0.6 MY t=0.8 MY t=1.0 MY

0.6 0.5 0.4 0.3 0.2 0.1 0 30

35

40

45

50

55

60

65

Gas mass concentration, kg/m3 Fig. 12. Case 2, gas mass concentration profile in the oil column.

70

semi-logarithmic plot given in Fig. 16. As shown in the figure, the density inversion only initiated after nearly 3 105 MY (30 years). The level of inversion continues growing and reaches a maximum of 0.95 kg/m3 at nearly 0.8 MY. As indicated in Section 2, a convection-free world is assumed in the simulation; however, in reality, the pure diffusion processes are inevitably interrupted by the buoyancy convections induced by the density inversion, particularly type II density inversion, which forms near the top of reservoir at the early stage of gas charge. In the remaining of this section, some rough comparisons will be made to understand the impact of the buoyancy convection on the pure diffusion process. The comparisons are made between the characteristic diffusional migration velocity of each component and the characteristic buoyancy convective velocity. These characteristic velocities are defined below. The diffusional migration velocity of each component, ui , can be calculated using Eq. (32) with known spatial and temporal composition variation in the oil column:

394

S. Pan et al. / Fuel 180 (2016) 384–395

The density profile

Ub

(a) Characteristic velocity, m/s

0.9 0.8

Normalized height

The evolution of characteristic velocities

-6

10

1

0.7 t=0.2 MY t=0.4 MY t=0.6 MY t=0.8 MY t=1.0 MY

0.6 0.5 0.4 0.3 0.2

UG -8

UA

10

UM -10

10

-12

10

0.1 -14

10

0 720

740

760

780

800

-4

-6

10

820

-2

10

10

0

10

Time, MY

Density ρ, kg/m3

Fig. 17. Case 2, the evolution of characteristic velocities.

The enlarged density profile ui ðt; zÞ ¼

1 0.9

(b)

Normalized height

0.7

0.5 0.4

t=0.2 MY t=0.4 MY t=0.6 MY t=0.8 MY t=1.0 MY

0.2 0.1 0 795

800

805

Density ρ, kg/m3 Fig. 15. Case 2, density profile in the oil column, (a) the whole profile and (b) enlarged profile with density inversion.

@ðci Þ dz ði ¼ G; A; or MÞ @t

ð32Þ

uG ¼

JG þ uM CG

ð33aÞ

uA ¼

JA þ uM CA

ð33bÞ

U i ðtÞ ¼ maxðjui ðt; 0 < z < Z T ÞjÞ ði ¼ G; A; or MÞ

1

ð34Þ

Here the absolute value is introduced to better compare the characteristic diffusional migration velocity of each component, among which uG and uA are always opposing uG and move toward the bottom of the oil column. On the other hand, the buoyancy convection is driven by the pressure gradient induced by the buoyant forces. The pressure gradient can be approximated as

dP Dqgzmax ¼ Dqg dz zmax

Density inversion

ð35Þ

where zmax is the location of maximum density in the oil column and g is the gravitational acceleration 9.8 m/s2. With the Darcy’s law, the characteristic buoyancy convective velocity U b is then estimated by

0.9

Density inversition Δρ, kg/m 3

0

z

In this work, uM is calculated by Eq. (32) whereas uG and uA are calculated using Eq. (33) to reduce numerical errors. The characteristic diffusional migration velocity of each component, U i , is then defined by the maximum of absolute value of ui in the oil column at any time.

0.3

0.8 0.7

k dP kDqg

¼

U b ¼

/l dz /l

0.6 0.5 0.4 0.3 0.2 0.1 0 -6 10

Z

Meanwhile, uG and uA can also be obtained from the definition of J G and J A in Eq. (5):

0.8

0.6

1 Ci

-2

-4

10

10

Time, MY Fig. 16. Case 2, evolution of density inversion.

0

10

ð36Þ

Here / is the porosity, k is the permeability, and l is the viscosity. Their values for the current study are given in Section 3.1. Fig. 17 shows the evolution of all aforementioned characteristic velocities in 1 MY. It is shown that all three characteristic diffusional migration velocities remain relatively stable and decrease slowly with lower mass transfer driving forces. They hold the relationship of U A > U G > U M since the asphaltenes is expelled by the filling of gas whereas the maltenes is a relatively inertial component in the system. It is also shown that U b quickly grows shortly after the density inversion starts forming. The magnitude of U b catches up with all

S. Pan et al. / Fuel 180 (2016) 384–395

the three characteristic diffusional migration velocities in about 4 105 MY (40 years), only 10 years after the formation of density inversion. In about 5 104 MY (500 years), U b appears stable and more than 100 times stronger than the diffusional migration velocities. It implies in the real world the diffusional gas charging process will be distorted in about 40 years and heavily suppressed by buoyancy convection in 500 years, far before the estimated time when the fluid in the reservoir achieves a thermodynamic equilibrium only through diffusional migration. The buoyancy convection induced by the density inversion, rather than diffusion, might be the dominant mechanism of hydrocarbon species redistribution after gas charge. It should be noted that in geological time, 40 years or even 500 years is far too short even for the assumption of an instantaneous charge; therefore, the reality must be more complicated than the ideal scenario simulated in the current case study. However, the results appear sufficient to illustrate the incipience of convection in the early stage of gas charge and its significant impact on the distribution of hydrocarbon species in reservoir. 4. Conclusions Asphaltene expulsion, the diffusional migration of asphaltenes driven by charge gas, appears to be an important mechanism that leads to asphaltene accumulation near the bottom of a reservoir and therefore the formation of a heavy oil layer and/or a tar mat at the base of an oil column. In this study, a multicomponent diffusion model has been proposed to simulate the diffusional migrations during gas charge processes in oil reservoirs. In the model, the mass diffusion is expressed in terms of the molar diffusion fluxes with the velocity of maltenes as the reference frame. Without any further simplification, these expressions have been adopted with a space coordinate transformation technique to convert this moving boundary diffusion problem into a fixed boundary problem with simple boundary conditions. Under the framework of non-equilibrium thermodynamics, some constitutive equations have been developed through the generalized Maxwell–Stefan theory to link the diffusion fluxes to the concentration gradient and make the equation set closed. The required thermodynamic properties are obtained from the Flory–Huggins regular solution model to account for the distributions of asphaltenes. Although this multicomponent diffusion model has been proposed based on a ternary system, it can be readily applied to systems of any number of components. The model has been used to study the gas charge processes in a typical hydrocarbon reservoir system through two case studies. In both studies, the asphaltene expulsion and density inversion phenomena have been clearly demonstrated. The simulations have also shown some reasonable adjustment in thermodynamic model in the two case studies can lead to two different types of density inversions. A type I density inversion occurs near the bottom of a reservoir at the late stage of gas charge. It is originated from a synergy of the asphaltenes expulsion and an impermeable bottom boundary, resulting in a concentration profile with a big, flat asphaltene hump near the reservoir bottom. A type II density inversion, which is dominantly caused by the asphaltene expulsion, initializes shortly after the beginning of the gas charge process and continues growing if there is no other disturbance. The nondisturbance assumption, however, is not realistic due to the presence of buoyancy convection induced by the type II density inversion near the top of a reservoir. In this work, some characteristic velocities have been defined to understand the potential impact of buoyancy convection. The comparisons have suggested that in the real world the diffusional gas charging process is distorted and heavily suppressed by buoyancy convection, well before the estimated time to achieve a thermodynamic equilibrium only through diffusional migration. It is implied that in the

395

presence of type II density inversion, the buoyancy convection, rather than diffusion, might be the dominant mechanism in hydrocarbon species redistribution during gas charge. References [1] Peters KE, Walters CC, Moldowan JM. The biomarker guide. 2nd ed. Cambridge University Press; 2005. [2] Head IM, Jones DM, Larter SR. Biological activity in the deep subsurface and the origin of heavy oil. Nature 2003;426:344–52. [3] Jackson RR, Zuo JY, Agarwal A, Herold BH, Kumar S, De Santo I, et al. Mapping and modeling large viscosity and asphaltene variations in a reservoir undergoing active biodegradation. In: SPE annual technical conference and exhibition, 27–29 October, Amsterdam, The Netherlands. SPE-170794-MS, Society of Petroleum Engineers; 2014. [4] Mullins OC, Zuo JY, Wang K, Hammond PS, De Santos R, Dumont H, et al. The dynamics of reservoir fluids and their substantial systematic variations. Petrophysics 2014;55:96–112. [5] Achourov V, Pfeiffer T, Kollien T, Betancourt SS, Zuo JY, di Primio R, et al. Gas diffusion into oil, reservoir baffling and tar mats analyzed by downhole fluid analysis, pressure transients, core extracts and DSTs. Petrophysics 2015;56:346–57. [6] Forsythe JC, Pomerantz AE, Seifert DJ, Wang K, Chen Y, Zuo JY, et al. A geological model for the origin of fluid compositional gradients in a large Saudi Arabian oilfield: an investigation by two-dimensional gas chromatography (GC GC) and asphaltene chemistry. Energy Fuels 2015;29:5666–80. [7] Stainforth JG. New insights into reservoir filling and mixing processes. Geol Soc London Special Publ 2004;237:115–32. [8] Mullins OC. The modified Yen model. Energy Fuels 2010;24:2179–207. [9] Pomerantz AE, Mullins OC, Zuo JY, Andrews AB. Asphaltenes explained for the nonchemist. Petrophysics 2015;56:226–75. [10] Freed DE, Mullins OC, Zuo JY. Theoretical treatment of asphaltene gradients in the presence of GOR gradients. Energy Fuels 2010;24:3942–9. [11] Freed DE, Mullins OC, Zuo JY. Heuristics for equilibrium distributions of asphaltenes in the presence of GOR gradients. Energy Fuels 2014;28:4859–69. [12] Mishra VK, Zuo JY, Dumont H, Mullins OC. Permeable tar mat formation described within context of novel asphaltene science. In: SPE Kuwait international petroleum conference and exhibition, 10–12 December, Kuwait City, Kuwait. SPE-163292-MS, Society of Petroleum Engineers; 2012. [13] Zuo JY, Elshahawi H, Mullins OC, Dong C, Zhang D, Jia N, et al. Asphaltene gradients and tar mat formation in reservoirs under active gas charging. Fluid Phase Equilib 2012;315:91–8. [14] Seifert DJ, Zeybek M, Dong C, Zuo JY, Mullins OC. Black oil, heavy oil and tar in one oil column understood by simple asphaltene nanoscience. In: Abu Dhabi international petroleum conference and exhibition, 11–14 November, Abu Dhabi, UAE. SPE-161144-MS, Society of Petroleum Engineers; 2012. [15] Pomerantz AE, Wu QH, Mullins OC, Zare RN. Laser-based mass spectrometric assessment of asphaltene molecular weight, molecular architecture, and nanoaggregate number. Energy Fuels 2015;29:2833–42. [16] Zuo YJ, Chen Y, Pan S, Wang K, Mullins OC. Investigation of density inversion induced by gas charges into oil reservoirs using diffusion equations. Energy 2016;100:199–216. [17] Mullins OC. The physics of reservoir fluids: discovery through downhole fluid analysis. Schlumberger; 2008. [18] Joshi NB, Mullins OC, Jamaluddin A, Creek J, McFadden J. Asphaltene precipitation from live crude oil. Energy Fuels 2001;15:979–86. [19] Zuo JY, Mullins OC, Mishra V, Garcia G, Dong C, Zhang D. Asphaltene grading and tar mats in oil reservoirs. Energy Fuels 2012;26:1670–80. [20] Zuo JY, Mullins OC, Dong C, Zhang D. Modeling of asphaltene grading in oil reservoirs. Nat Resour 2010;1:19–27. [21] Zuo JY, Mullins OC, Freed D, Elshahawi H, Dong C, Seifert DJ. Advances in the Flory–Huggins–Zuo equation of state for asphaltene gradients and formation evaluation. Energy Fuels 2013;27:1722–35. [22] Wang JX, Buckley JS. A two-component solubility model of the onset of asphaltene flocculation in crude oils. Energy Fuels 2001;15:1004–12. [23] Mohammadi AH, Richon D. A monodisperse thermodynamic model for estimating asphaltene precipitation. AIChE J 2007;53:2940–7. [24] Alboudwarej H, Akbarzadeh K, Beck J, Svrcek WY, Yarranton HW. Regular solution model for asphaltene precipitation from bitumens and solvents. AIChE J 2003;49:2948–56. [25] Taylor R, Krishna R. Multicomponent mass transfer. Wiley; 1993. [26] Kjelstrup S, Bedeaux D, Johannessen E, Gross J. Non-equilibrium thermodynamics for engineers. World Scientific; 2010. [27] Crank J. The mathematics of diffusion. Clarendon Press; 1979. [28] Duda JL, Vrentas JS. Mathematical analysis of sorption experiments. AIChE J 1971;17:464–9. [29] Alsoy S, Duda JL. Influence of swelling and diffusion-induced convection on polymer sorption processes. AIChE J 2002;48:1849–55. [30] Prausnitz JM, Lichtenthaler RN, Azevedo EGd. Molecular thermodynamics of fluid-phase equilibria. 3rd ed. Upper Saddle River, N.J.: Prentice-Hall PTR; 1999 [c1999]. [31] MATLAB. The MathWorks, Inc., Natick, Massachusetts, United States; 2012b. [32] Pisani L. Simple expression for the tortuosity of porous media. Transp Porous Med 2011;88:193–203.