Analysis of Asphaltene Instability Using Diffusive ... - ACS Publications

Subscriber access provided by UB + Fachbibliothek Chemie | (FU-Bibliothekssystem)

Article

Analysis of Asphaltene Instability Using Diffusive and Thermodynamic Models during Gas Charges into Oil Reservoirs Julian Youxiang Zuo, Shu Pan, Kang Wang, Oliver C. Mullins, Hadrien Dumont, Li Chen, Vinay Mishra, and Jesus A Canas Energy Fuels, Just Accepted Manuscript • DOI: 10.1021/acs.energyfuels.6b03305 • Publication Date (Web): 23 Feb 2017 Downloaded from http://pubs.acs.org on February 24, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Energy & Fuels is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Analysis of Asphaltene Instability Using Diffusive and Thermodynamic Models during Gas Charges into Oil Reservoirs Julian Y. Zuo1*, Shu Pan2, Kang Wang3, Oliver C. Mullins4, Hadrien Dumont5, Li Chen5, Vinay Mishra5, Jesus Canas5 1

HFE, Schlumberger, 150 Gillingham Lane, Sugar Land, TX 77478, USA SRC, Schlumberger, 14910 Airline Road, Rosharon, TX 77583-1590, USA 3 BGC, Schlumberger, Chuang Xin Building Beijing Tsinghua Science Park, Beijing100084, China 4 Schlumberger-Doll Research, Cambridge, MA 02139, USA 5 Schlumberger, 1325, S. Dairy Ashford, Houston, TX 77077, USA * Corresponding author: [email protected] 2

Abstract Reservoir fluid analysis can be used to better understand a hydrocarbon reservoir. Reservoir fluids are not always equilibrated in particular because of dynamic processes such as late gas charges occurring in the recent geologic past. It is of great importance to delineate how a late gas charge affects fluid distributions and when asphaltenes become unstable in oil. In this paper, we modify the simplified moving boundary 1-D diffusive model of Zuo et al. (2016) for gas charges into oil reservoirs considering the Maxwell-Stefan diffusivity and thermodynamic nonideality. Live oil is simply treated as three pseudocomponents: gas, maltene and asphaltene. The parameters of the three components in the Flory-Huggins regular solution model are calibrated by the Peng-Robinson EOS. The modified diffusive model is then used to simulate gas charges into oil reservoirs. The simulation results have indicated that fluid density inversion is created by competing effects of asphaltene diffusion and expulsion owing to late gas charges. The density inversion induces gravity currents which cause asphaltene migration from near the top to the base of the oil column fairly rapidly in geological time. In addition, the same thermodynamic model as in the diffusive model is also used to generate vapor-liquid equilibrium (VLE), liquid-liquid equilibrium (LLE) and vaporliquid-liquid equilibrium (VLLE) spinodal and binodal phase boundaries. Furthermore, the fluid compositional variation paths simulated by the modified diffusive model and the calculated phase boundaries are plotted in the same phase diagrams for comparison. The results show that asphaltenes can precipitate locally near the gas/oil contact (GOC) for fluids with high asphaltene content because rapid gas addition to an oil reservoir at the GOC easily makes fluids near the GOC metastable and/or unstable (high gas solution oil decreases its solvency capacity to dissolve asphaltenes). Moreover, with asphaltene-enriched advective flow, asphaltene buildup at the base (oil/water contact, OWC) can give rise to asphaltene phase instability and tar mat formation when the asphaltene content surpasses the maximum asphaltene solubility allowable in the oil. This is consistent with the frequent field observations of asphaltene deposition at the upstructure with shale breaks and the OWC. Additionally, the sensitivity of the thermodynamic model parameters to phase envelopes is analyzed. The sensitivity analysis indicates that decreasing gas and maltene solubility parameters, temperature, and maltene molar mass, and increasing asphaltene solubility parameters and asphaltene molar mass result in increasing asphaltene phase instability. Keywords: Diffusive model, Gas charge, Density inversion, Asphaltene instability, Spinodal Curve, Binodal Curve Introduction Reservoir fluid analysis can be used to better understand a hydrocarbon reservoir. A reservoir oil can be simply treated as dissolved gas, liquid and dissolved solid (asphaltenes). Downhole fluid analysis (DFA) measures reservoir fluid properties downhole in real time [1], which include compositions of CO2, C1, C2, C3, C4, C5, C6+, gas/oil ratio (GOR), density, viscosity, optical density (linearly related to asphaltene content), temperature, pressure, and so on. Under the framework of the Yen-Mullins model [2], the Flory-Huggins-Zuo equation of state (FHZ EOS) [3, 4] was developed for estimating equilibrium asphaltene gradients in oil columns. By integrating DFA, the Yen-Mullins model, and the FHZ EOS, we have successfully described equilibrated asphaltene gradients and predicted reservoir connectivity. That is, the equilibration of asphaltenes requires significant mass transfer in reservoirs which can only occur if the reservoir has good connectivity. [5] Fig. 1 summaries some typical field case studies from light to heavy oil. The sizes of heavy ends change from a true molecular solution of perylene-like resin molecules [6], a true molecular solution of asphaltene (or asphaltene-like) molecules [7], asphaltene nanoaggregates [8, 9] to asphaltene 1

ACS Paragon Plus Environment

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

clusters [10]. The equilibrium distributions of heavy ends in formation oils have suggested the reservoir connectivity which has been proven by the subsequent production data. [5-10]

Fig.1 Gradients of resins and asphaltenes in five different oil reservoirs (a-e from left). (a) volatile oil with nearly no asphaltenes [6]; (b) condensate oil with asphaltene or asphaltene-like molecules[7], (c) intermediate-GOR oil with asphaltene nanoaggregates [8]; (d) low-GOR black oil with asphaltene nanoaggregates [9]; (e) mobile heavy oil with asphaltene clusters [10].

The case studies in Fig. 1 show how important it is to delineate the equilibrium distributions of heavy ends in reservoir fluids in light of the DFA measurements using the FHZ EOS and the Yen-Mullins model. In addition, this same workflow can also be used to identify non-equilibrium asphaltene distributions of reservoir fluids. Reservoir fluids are not always equilibrated. Dynamic processes occurring in reservoirs can give rise to non-equilibrium fluid distributions over geologic time, as observed in many fields. [11 - 14] One dynamic process is a late gas charge into an oil reservoir. In a typical scenario of basin subsidence, the kerogen in source rock generates oil in the oil temperature window and then generates gas in the gas window at higher temperatures while subsiding. The generated oil migrated into the reservoir and subsequently the generated gas entered the oil reservoir over geological time. As such, this process gives rise to a late gas charge into an oil reservoir. In addition, a late oil charge into a gas reservoir would yield a similar initial condition of a gas cap over an undersaturated oil. Fig. 2 depicts a field case of a gas charge into an oil reservoir. [11] The reservoir fluids depicted in Fig. 2 are “density stacked” without mixing [15]. A gas cap is formed above an undersaturated oil column. Gas then diffuses down into the oil column, which results in an increase in both solution gas (or gas/oil ratio, GOR) and saturation pressure in the top half of the oil column as shown in Fig 2C. The enrichment of solution gas resulted in asphaltene expulsion, which is confirmed by the fact that the oil samples from the top of the column are nearly colorless as shown in Fig. 2A. The expelled asphaltenes migrate to the base of the column as shown Fig. 2B. In this tilted sheet reservoir with a relatively small dip angle, most of the distance from the top to bottom is horizontal instead of vertical. Thus, the change in true vertical depth (TVD) takes into consideration only a small part of the total diffusion distance. The large horizontal distance slows down the diffusion process to geologic timescales. If this gas charge process continues, more asphaltene moves down to the base of the reservoir and a tar mat can form. [12]

2


Page 2 of 19

Page 3 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Fig.2. Gas charge into an oil reservoir [11]. A-Left) Illustration of the increase in solution gas, especially near the gas oil contact, due to gas diffusion into oil. A-Right) Series of dead oil samples from a single reservoir deep water Gulf of Mexico reflecting this process. B) Huge asphaltene gradient (from DFA and lab measurements) showing the expulsion of asphaltenes from the top of the column and accumulation of asphaltene at the base. C) GOR and saturation pressure show huge nonequilibrium gradients associated with gas diffusion into primarily the top of the oil column. Both diffusion and advection are important mechanisms in asphaltene transport. Seifert et al. observed asphaltene transport in very large reservoirs where diffusive processes would be far too slow to account for field measurements of asphaltene accumulations. [16] Because the thermal gradients are too small to give advective overturn in reservoir fluids, the advection here is driven by asphaltene expulsion from the top of the oil column. Asphaltene expulsion due to the gas charge yields density inversion, which, in turn, gives rise to gravity currents (advection). Nevertheless, it is required to show that such a density inversion results from a rigorous treatment. A similar analysis would apply if a light oil charges into a heavier oil reservoir. Such an occurrence can yield some asphaltene instability and migration, but can leave a black oil in the crest of the field. In addition, as mentioned previously, late gas or light hydrocarbon charges into oil reservoirs often result in an increase in GOR and thus a decrease in maximum asphaltene solubility allowable in the oil. When asphaltene content exceeds its solubility in oil, a separate asphaltene phase may occur. The accumulation of asphaltenes leads to the formation of tar mats, which has significantly impact on reservoir quality and field development planning. Therefore, it is of great importance to delineate how gas charges affect fluid distributions and when asphaltenes are unstable in oil. Zuo et al. (2016) developed a 1-D diffusive model with a moving boundary problem to simulate gas charges into oil reservoirs. [17] In their model, a presumption of constant total molarity in a general diffusion equation is made to obtain the simplified diffusive model. Pan et al. (2016) then introduced a new space variable, and rigorously converted the complicated moving boundary problem into a fixed boundary problem with simple boundary conditions. [18] However, this phase variable transformation cannot be easily extended to solve 2-D/3-D problems, therefore it is limited only for 1-D studies. In the first part of present work, a modification is proposed to improve the simplified diffusive model and make it suitable for 2-D/3-D problems in the future studies. [17] On the other hand, phase envelopes are widely used to guide reservoir fluid productions. Phase instability boundaries include spinodal and binodal curves. There are many publications that describe how to obtain phase envelopes and equilibrium calculations in the open literature. [19-23] However, it is still challenging to reliably generate phase boundaries for highly asymmetric gas + solvent (maltene) + asphaltene mixtures. Hence, in the second part of this paper, a computational thermodynamic algorithm is proposed for generating vapor-liquid, liquidliquid and vapor-liquid-liquid spinodal and binodal boundaries for gas + maltene + asphaltene mixtures. With this 3


Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 19

approach, we can compare diffusive compositional variation paths during gas charges into oil reservoirs with the multiphase boundaries to find out where asphaltenes are unstable and better understand the effects of gas charges on asphaltene phase instability. Diffusive Model for Gas Charges into Oil Reservoirs Fig. 3 shows a schematic diagram of a gas charge into an oil reservoir. The reservoir has a gas cap above the oil column and an aquifer below the oil column. To simplify the problem, one-dimension (1-D) diffusion system is taken into account. The reservoir is simply treated as either inclined or vertical column with fixed cross-sectional area. The origin is located at the oil/water contact (OWC). The reservoir hydrocarbon fluid is treated as three pseudocomponents: gas (1, mainly methane or pure methane), asphaltene (2, solute), and maltene (3, solvent which is the rest of components excluding gas and asphaltene in oil). Thermodynamic nonideality of the mixture is described by the Flory-Huggins regular solution model [24] and to be consistent, the same thermodynamic model is used for both diffusion and phase instability. The aquifer is presumed to be incompressible and impermeable for diffusion of hydrocarbon components. It is also presumed that the gas cap is pure gas and asphaltene component cannot move to the gas cap because gas does not have solvency capability to dissolve asphaltene component. At the gas/oil contact (GOC), gas and oil are presumed always in equilibrium. Initial conditions of the oil column can be homogeneous or heterogeneous. Because gas addition to the oil column results in a decrease in gas cap volume and an increase in oil column volume, a moving boundary diffusion problem is involved. To avoid complicating the problem, we ignore sources (or sinks), buoyancy or forced advections, chemical reactions, gravitational and thermal fluxes.

Fig. 3 Schematic diagram of 1-D diffusion for a gas charge into an oil column reservoir. Homogeneous initial conditions are assumed (left); at any time t, gas diffuses down, GOC moves up, and asphaltene migrates down (right). Where t, θ, h, and z are the time, dip angle, true vertical depth (height), and 1-D diffusion path, respectively.

Based on the aforementioned presumptions, a general diffusion governing equation for an N-component mixture (N1 independent variables) is expressed as

∂C = ∇ ⋅ (C t [ D ]∇x ) = ∇ ⋅ (C t [ B ]-1 [ Γ ]∇x ) ∂t

(1)

where C and x are the vectors of the molarity (molar concentration) and the mole fraction, t and Ct denote the time and the total molarity. ∇⋅ is the divergence operator. The diffusion equation can also be used for 1-D, 2-D and 3-D. The matrix [B] of the drag effects is given by the following expression of the Maxwell–Stefan diffusivities (Dij for the i-j pair diffusivity): [25] Bii =

xi + DiN

N

xk

∑D k =1 k ≠i

ik

,

 1 1  Bij (i≠ j ) = xi  − , D   iN Dij 

i, j = 1,2,..., N − 1

(2)

The thermodynamic nonideality of the mixture (matrix [Γ ]) is computed by using the derivatives of the activity coefficients which will be described later [25]:

4


Page 5 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

 ∂ ln γ i Γij = δ ij + x i   ∂x j 

    T ,P , x

i , j = 1, 2,..., N − 1

(3)

k≠ j

where γi and δij are the activity coefficient of component i and the Kronecker delta function for the i-j component pair. For ideal mixtures the activity (or fugacity) coefficients are unity, γi = 1, [Γ Γ] = [I]. Then Eq. (1) is reduced to Fick’s law but diffusion coefficients are given by [B]-1. In our previous work [17], Ct is presumed to be unchanged in Eq. (1) so that it can be canceled out on both sides of the equation. Hence, the following simplified 1-D diffusion equation is obtained:

∂x ∂  ∂x  ∂  ∂x  = [ D]  = [ B]-1[ Γ ]  ∂t ∂z  ∂z  ∂z  ∂z 

(4)

However, in present work, we remove such a presumption to make our model more rigorous. Therefore, the 1-D diffusion equation is given by:

∂C ∂  ∂x  ∂  ∂x  =  Ct [ D]  =  Ct [ B]-1[ Γ ]  ∂t ∂z  ∂z  ∂z  ∂z 

(5)

We convert molarity in Eq. (5) to mole fraction by using the following expression:

xi = C i v =

Ci vv N = vN

Ci vv N N −1

v−

∑ x (v k

k

− vN )

k =1

Ci v N

=

N −1

1−

∑ C (v k

k

i = 1,2,..., N − 1

− vN )

(6)

k =1

where the total molar volume (v) is calculated by the summation of the mole fraction (xk) multiplied by the constant component partial molar volume (vk): N

N −1

k =1

k =1

v = ∑ xk vk = ∑ xk (vk − vN ) + vN

(7)

In this work, we take into account the oil column only. Because the GOC is moving up during the gas charge, we are dealing with a moving boundary problem. For simplicity, initial conditions (t = 0) are set to be homogeneous compositions in the one dimension oil column, i.e.,

Ci (0, z ) = Ci0 ,

i = 1,2 ,..., N

(8)

The boundary conditions are given as follows. At the base of the oil column (z = 0, OWC), impermeable diffusion boundary conditions are applied

∂Ci (t ,0) = 0, ∂z

i = 1,2 ,..., N − 1

(9)

As mentioned previously, a reservoir fluid is grouped into three pseudocomponents: gas (1) + asphaltene (2) + maltene (3). Because gas component is always in equilibrium with oil at the GOC and no asphaltene component moves up from the oil column to the gas cap, we have

∂C 2 (t , GOC ) =0 ∂z where superscript sat refers to the saturated (equilibrium) condition at the GOC. C1 = C1sat ,

and

5


(10)

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 19

Because gas addition into the oil column gives rise to swelling of the oil column, the GOC changes can be estimated by [17] t

z = z0 +

∫C 0

J1GOC dτ − C1GOC −G

(11)

GOC −O 1

where z0 is the initial depth of the GOC. J1 and C1 are the gas flux and molar concentration at the GOC. G and O denote the gas and oil sides of the GOC. The 1-D diffusion model is used for simulation of single phase diffusion problems. Although initial and boundary conditions can be set differently according to actual situations, Eqs. (8) – (10) are applied to this paper as examples. Combined these equations above with the initial and boundary conditions, composition variations with time in the oil column can be calculated numerically. The calculation is conducted by solving the partial differential equations (PDE) mentioned above using a commercial PDE solver with 200 and more spatial and temporal grids.

Phase Boundary Calculations Gas charges into oil reservoirs may result in asphaltene phase instability because gas addition decreases oil solubility parameters and lowers solvency capacity of oil to dissolve asphaltenes. Asphaltene flocculation complicates reservoir evaluation and field development planning and also may cause serious production problems. To ensure asphaltenes to stay in oil stably, it is of great importance to estimate phase envelopes (boundaries) of the mixture in question so that we know what conditions give rise to asphaltene phase instability. In thermodynamics, phase instability boundaries include spinodal and binodal curves. The spinodal condition is the limit of local stability with respect to small fluctuations, which is clearly defined by the condition that the second derivative of Gibbs free energy is zero. The binodal condition is defined as the limit of global minimum energy equilibrium state of the system. The spinodal instability region is less than the binodal instability region. The binodal instability region subtracted by the spinodal instability region is the metastable region. The binodal and spinodal meet at the critical (plait) point. In the following sections, we describe the methods to generate both spinodal and binodal loci. Spinodal Curve Calculations According to thermodynamics, the spinodal criterion for an N-component system at specified temperature and pressure is that the determinant of the Jacobian matrix of molar Gibbs free energy is equal to zero[21], i.e.,

∂2 g ∂x12 ∂2 g det ∂x ∂x 2 1 ... ∂2 g ∂x N −1∂x1

∂2 g ∂x1∂x2 ∂2 g ∂x22 ... ∂2 g ∂x N −1∂x 2

∂2g ∂x1∂x N −1 ∂2g ... ∂x2 ∂x N −1 = 0 ... ... ∂2g ... ∂x N2 −1 ...

(12)

where g is the molar Gibbs free energy that is a function of temperature (T), pressure (P) and mole fraction (x). The molar Gibbs free energy is related to the chemical potential (µ), Eq. (12) is then rewritten as

∂µ1 ∂x1 ∂µ 2 det ∂x 1 ... ∂µ N −1 ∂x1

∂µ1 ∂x 2 ∂µ 2 ∂x 2 ... ∂µ N −1 ∂x 2

∂µ1 ∂x N −1 ∂µ 2 ... ∂x N −1 = 0 ... ... ∂µ N −1 ... ∂x N −1 ...

(13)

The chemical potential is further calculated by an activity coefficient model 6


Page 7 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

µi = µiref + RT ln xiγ i

(14)

where superscript ref is the reference state and γi denotes activity coefficient of component i. Thus, Eq. (13) can be expressed as

det [Γ' ] = 0

(15)

where the matrix [Γ′ ] is associated with [Γ ] in diffusion Eq. (1) or Eq. (3) by Γij' =

 ∂ ln γ i + xi  ∂x j

δ ij

 Γ  = ij  xi T , P , x

i, j = 1,2,..., N − 1

(16)

k≠ j

Therefore, one can determine spinodal instability directly by evaluating the determinant value of matrix [Γ′ ]. If det|[Γ′ ]| < 0, the mixture is unstable (spinodal decomposition); otherwise it is stable (or metastable). If an EOS model is used, the activity coefficients are replaced by the fugacity coefficients. For a ternary mixture, there are three unknowns (x1, x2 and x3) at specified temperature and pressure but there are only two equations which include Eq. (15) and x1 + x2 + x3 = 1. Therefore, if the mole fraction of one component is specified, then the other two mole fractions can be determined from these two equations. At a given temperature and pressure, a series of mole fractions of the first component is specified, the mole fractions of the second and third components are calculated using the Newton iteration method with proper initial guesses until Eq. (15) and the summation equation are finally satisfied. It should be noted that the spinodal calculation is based on a two-phase assumption.

Binodal Curve Calculations We treat the asphaltene phase as a liquid phase. Binodal (coexistence equilibrium) curves are determined in terms of phase equilibrium criterion. According to thermodynamics, vapor-liquid-liquid three phase equilibrium (VLLE) criterion is given by [24] f i V = f i L1 = f i L 2 ,

i = 1,2,..., N

(17)

where fi is the fugacity of component i and superscripts V, L1 and L2 denote the vapor, liquid 1 and liquid 2 phases. For simplicity, if the nonideality of the vapor phase and two liquid phases is described by the activity coefficient model (the model parameters are calculated by a cubic equation of state at reservoir conditions), the phase equilibrium criterion is rewritten as y i γ iV = xiL1γ iL1 = xiL 2 γ iL 2 ,

i = 1,2,..., N

(18)

where yi, and xi are the mole fractions of component i in the vapor and liquid phases. It should be noted that the activity coefficient model instead of an EOS is used for the vapor phase. This is because we use the same model for all the phases to avoid complications and furthermore the model parameters are calibrated by the cubic EOS at reservoir conditions so that the activity coefficient model is suited for the vapor and liquid phases. The reference state for all the components is either hypothetical or real pure liquid at the specified temperature and pressure. The methodology can be easily extended to other thermodynamic frameworks such as the EOS approach and the mixed EOS-activity coefficient model approach. For specified reservoir fluid compositions (z), the tangent plane distance algorithm is used for phase stability check [19] and the flash algorithm is employed for vapor-liquid equilibrium (VLE), liquid-liquid equilibrium (LLE) and VLLE calculations [19]. To improve the convergence issues in the vicinity of the critical point, Michelsen’s phase envelope algorithm [20] is extended to calculate the liquid-liquid coexistence curves and plait (critical) points. For LLE, we have N + 1 equations which are given by

7


Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 19

g i (α , β ) = ln xiL1 + ln γ iL1 − ln xiL 2 − ln γ iL 2 = ln K i + ln γ iL1 − ln γ iL 2 = 0 ,

i = 1,2,..., N

(19)

where Ki is the liquid-liquid equilibrium constant which is dependent on temperature, pressure and compositions, estimated by the activity coefficient model. The (N+1)th equation is the Rachford-Rice equation given by − 1)z i =0 ) i =1 i =1 i −1 where F is the mole fraction of the L1 phase. Adding one more specification equation, we have g N +1 (α , β ) =

∑ (x N

L1 i

− xiL 2 ) =

(K

N

∑ 1 + F (K i

g N +2 (α, β ) = α k − S = 0

(20)

(21)

The vector of specified variables are given by β T = (F , T , P )

(22)

The vector of dependent variables are

α T = (ln K 1 , ln K 2 ,..., ln K N , ln z j , ln z m )

(23)

Let zm be a specified variable and zj one of the primary variables, the remaining compositions are calculated by

zi = 1 − zm − z j ,

i ≠ j, m

(24)

Because the Newton-Raphson method is used to solve the set of nonlinear equations, good initial guesses are required. Therefore, the same extrapolation method is utilized as Michelsen. [20] A polynomial extrapolation is used to obtain a good initial estimation for kth iteration

α k (S ) = α k 0 + α k 1S + α k 2 S 2 + α k 3 S 3

(25)

where akj, j = 0,1, ... ,3, are the coefficients of the polynomial regressed from the information of the last two points. This extrapolation method needs two points to initialize. The LLE flash algorithm is used to generate the first point. A slight composition change is made for component m to obtain the second point. Eq. (25) is then used afterward. At plait (critical) points, Ki = 1 for i = 1, 2, …, N.

Activity Coefficient Model In general, there two approaches to model asphaltene phase separation behavior: EOS and activity coefficient. The PC SAFT EOS and the Flory-Huggins type models have broadly been used in modeling of phase behavior of asphaltene precipitation in the oil and gas industry. [26-31] In this work, the activity coefficient is calculated by the Flory-Huggins regular solution model [24]

ln γ i = ln

vi v v +1− i + i v v RT

∑∑ φ φ [(δ N

N

j

j =1

k

− δ j ) + 2lij δ i δ j − 0.5(δ j − δ k ) − l jk δ j δ k 2

i

2

]

(26)

k

where R, T, and δ is the gas constant, the temperature, and the solubility parameter and ljk is the binary interaction parameter between components j and k. For pure component j, ljj = 0. If ljk = 0, the Eq. (26) is reduced to vi v v (δ − δ ) +1− i + i i v v RT

2

ln γ i = ln

(27)

8


Page 9 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

N

δ = ∑ φk δ k

(28)

i =1

φi =

xi vi

(29)

N

∑ xk vk i =1

The molar volume (v) is estimated by Eq. (7).

Results and Discussions Basic Properties of Initial Reservoir Oil Before the aforementioned 1-D diffusive and thermodynamic models are applied to simulate gas charges into oil reservoirs, basic input properties of reservoir fluids are required. For simplicity, a live oil with GOR of ~50 sm3/sm3 and API gravity of ~37 °API is used as the initial reservoir fluid. The reservoir temperature and pressure are 339 K and 30 MPa, respectively. The reservoir fluid is characterized by the method of Zuo and Zhang [32] and then grouped into three pseudocomponents: gas + asphaltene + maltene. In this work, we presume gas is pure methane and the rest components excluding the asphaltene component are lumped as the maltene component. The experimental pressure-volume-temperature (PVT) data are matched by tuning the Peng-Robinson equation of state (PR EOS) parameters [33]. The density and partial molar volume of gas and maltene components are calculated by the PR EOS [33] at the specified reservoir conditions after EOS model tuning. The saturated gas solubility in the oil at the reservoir condition is also calculated by the PR EOS, which is equal to 0.575 in mole fraction. The solubility parameters of gas and maltene components are calculated by [34]

δ i = 17.347 ρ i + 2.904

(30)

where ρi is the density of component i in g/cm3 and δi is the solubility parameter of component i in MPa0.5. The solubility parameter of asphaltene component in MPa0.5 is calculated by [5]

δ = δ 0 (1 − 1.07 × 10 −2 (T − 298.15))

(31)

where δ0 is 21.85 MPa0.5 for asphaltenes and T is the temperature in K. Asphaltenes are present in the oil as nanoaggregates for black oil. [5] Based on the Yen-Mullins model [2], the size of asphaltene nanoaggregates is 2 nm, thus corresponding to molar volume of 2523 cm3/mol and molar mass of 3027 g/mol if density is 1.2 g/cm3.[5] It is assumed that the basic parameters of each component are constant at the specified reservoir pressure and temperature as given in Table 1. The initial oil column has homogeneous mole fractions for all the components as shown in column 1 of Table 1 and the gas cap only has pure methane at the beginning. At the GOC, no asphaltene flux is presumed because pure methane in the gas cap does not have solvency capacity to dissolve any asphaltenes and the saturated methane mole fraction is 0.575.

Table 1 Basic Properties of Initial Reservoir Oil at 339 K and 30 MPa Initial Molecular Specific Component Partial molar volume composition weight gravity mol% g/mol cm3/mol

Solubility parameter MPa0.5

Methane (gas)

33.16

16

0.314

51

8.35

Maltene

66.66

191

0.84

227

17.48

Asphaltene

0.28

3027

1.2

2523

20.89

9


Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 19

Diffusion and Density Inversion Methane migration diffusion coefficients in different real source rocks were measured by Krooss, which is ~ 1.62×10-11 m2/s. [35] This value is taken as diffusion coefficient for methane–maltene pair: D13= D31 = 1.62×10-11 m2/s. In addition, methane diffusion coefficient of (~1.3×10-10 m2/s) in Athabasca bitumen at 50 – 75 °C and 8 MPa was measured by Sheikha et al. [36] Tortuosity of porous media often reduces diffusion coefficients. In typical sandstones, tenfold reduction of diffusion coefficients is observed compared to those in free volume. [37] Therefore, diffusion coefficient for methane–asphaltene pair is: D12= D21 = 0.9×10-11 m2/s in porous media. Because there is no data available for maltene–asphaltene pair in the open literature, diffusion coefficient for maltene-asphaltene pair is set to: D23= D32 = 0.5×10-11 m2/s in porous media, which is smaller than that for methane–asphaltene pair. The diffusion model is at first used to simulate a gas charge process into an initially 100 m thick vertical oil column. We adjust l13 = l31 = -0.01 in the activity coefficient model to match the gas solubility in oil at the GOC, which is estimated by the PR EOS (0.575 in mole fraction). We assign l12 = l21 = 0.01, and l23 = l32 = -0.015 (more discussions are given later). It should be noted that the results simulated by the modified diffusive model are just slightly different from those by the simplified diffusive model [17]. For example, the GOC moves up ~ 2.2 m over one million years (MY) by the modified model whereas ~2 m by the simplified diffusive model. In addition, because Zuo et al. gave the detailed results and sensitivity studies in reference [17] for the simplified diffusive model, we only briefly summarize main results obtained by the modified diffusive model herein. The simulated profiles of methane mole fraction, asphaltene mole fraction and fluid density with initial asphaltene mole fractions of 0.0028 and 0.0024 (corresponding to 6.0 and 5.2 wt%) are shown in Fig. 4.

Fig. 4. Simulation results for a 100 m vertical reservoir. (a) Simulated variation of methane mole fraction with initial x2=0.0028; (b) Simulated variation of asphaltene mole fraction with initial x2=0.0028; (c) Simulated variation of fluid density with initial x2=0.0028; (d) Simulated variation of methane mole fraction with initial x2=0.0024; (e) Simulated variation of asphaltene mole fraction with initial x2=0.0024; (f) Simulated variation of fluid density with initial x2=0.0024;Obviously, density inversion is generated by the diffusion model in a small range of initial asphaltene content. It is observed from Fig. 4 (a) and (d) that gas does not reach the base of the oil column over one million years because methane mole fraction at the base is still unchanged. Asphaltene accumulation occurs because the opposite effects of asphaltene diffusion from high to low concentration and asphaltene expulsion due to the gas charge. That is, the gas charge increases solution gas in oil and lowers its solvency capacity to dissolve asphaltenes, thus 10


Page 11 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

expelling asphaltenes out from the oil. Consequently, asphaltenes build up somewhere behind the gas diffusion. This asphaltene accumulation starts near the GOC and then moves downward with time as shown in Fig. 4 (b) and (e). Density inversion is generated due to the gas charge into the oil column as shown in Fig. 4 (c) and (f). The extent of density inversion depends on initial asphaltene concentrations. The higher asphaltene content the larger density inversion. If initial asphaltene mole fraction is set to 0.0024 (5.2 wt%), the maximum density inversion is < 1 kg/m3 whereas if it is set to 0.0028 (6.0 wt%), the maximum density inversion is ~3 kg/m3. No density inversion can be created when initial asphaltene mole fraction is less than 0.0022 (4.8 wt%). Because we assume constant partial molar volumes for each component, the density inversion is obviously caused by composition differences, in particular asphaltene accumulation. On the other hand, increasing initial asphaltene content further results in asphaltene phase instability. Large density inversion often occurs near the asphaltene phase instability boundary, in particular in the metastable region in between spinodal and binodal boundaries, which will be discussed in more detail thereafter. At the GOC, fluid density is much lower than the initial fluid density because higher equilibrium gas concentration lowers fluid density. At a given time, fluid density increases rapidly (almost linearly) along the depth and reaches a maximum point where fluid density is slightly (∆ρ ≈ 3.0 kg/m3) higher than the initial fluid density. Then fluid density gradually decreases below the initial density value and subsequently increases to the initial density value. Evidently, fluid density inversion is created. The density inversion generated by the gas charge can induce density driven advection (gravity currents), thus moving asphaltenes down to the base of the oil column. Time scale of gravity currents induced by density inversion can be compared to diffusion briefly below. Buoyancy advection velocity induced by density inversion can be estimated by [38, 39]

u=

∆ρgk sin θ Φµ

(32)

For example, if a reservoir has permeability k = 300 mD, porosity Φ = 0.2, viscosity µ = 1 cP, dip angle θ = 5°, and density inversion ∆ρ = 3 kg/m3 the estimated gravity current velocity is 37 km per million years, which is much faster than diffusion (~100 m per million years). Therefore, the created density inversion can migrate asphaltenes to the base of the oil column over geological time. It should be noted that asphaltene phase instability can occur during the gas charge process. In this work, asphaltene phase instability is thermodynamically analyzed by generating spinodal and binodal boundaries for the same mixture using the same thermodynamic model as in the diffusion simulation. Then we compare the simulated compositional variation paths during the gas charge process and find out when and where asphaltene phase instability occurs.

Phase instability analysis First of all, the phase diagrams of the three constitutive binary mixtures for the ternary system mentioned above are calculated using the aforementioned method. The calculated binodal and spinodal curves are shown in Fig. 5. In general, the mixtures stay in a single stable phase outside the binodal curves. These binodal curves are also called the coexistence curves. The mixtures are separated into two distinct (vapor-liquid or liquid-liquid) phases inside the spinodal curves, which is referred to as spinodal decompositions where a new phase is formed spontaneously. The metastable region is in between the binodal and spinodal curves. It should be noted that the spinodal and binodal curves depend on the thermodynamic model and its parameters but they are independent of the aforementioned diffusive model. Furthermore, the diffusive model relies not only on the thermodynamic model but also on diffusive model parameters. The diffusive model gives compositional variation paths (composition profiles with depth and time) during the gas charge into the oil column. Combining those together, we can know when and where compositional paths cross the phase instability boundaries. To check the sensitivity of phase diagrams to the binary interaction parameters (lij) in the activity coefficient model for the three constitutive binary mixtures, the phase diagrams with and without binary interaction parameters are compared in Fig. 5. It can be seen that positive binary interaction parameters (lij) enlarge metastable and unstable regions whereas negative values decrease metastable and unstable regions.

11


Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 19

In addition, according to Prausnitz et al.[24], the binary interaction parameter (lij) in the regular solution model is essentially a correction to the geometric mean rule for the cohesive energy density (or to the solubility parameter). The geometric mean rule assumes to be strictly valid for non-polar molecules according to the London theory. Thus, this lij is similar to the well-known kij in cubic equations of state. Unfortunately, in most cases, the lij parameter cannot be correlated with physical properties of the compounds in mixtures. However, some rough approximations have been proposed for aromatic-saturated hydrocarbon mixtures, where lij’s are in a range of -0.03 to 0.015, as given by Prausnitz et al.[24] As mentioned previously, l13 = l31 = -0.01 was obtained by matching the saturated gas solubility in oil at the GOC (x1 = 0.575) obtained from the PR EOS. This l13 value falls into the range given by Prausnitz et al. [24]. It is interesting that maltene and asphaltene components can be completely miscible at any given compositions as shown in Fig. 5(c) when temperature above 320 K by changing l23 = l32 from 0 to -0.015. It can be seen that the immiscible region has very high asphaltene contents (x2 > 0.4 when T slightly greater than 320 K with l23 = 0, corresponding to > 90 wt% asphaltene). Oil with such high asphaltene concentration (> 90 wt%) is not mobile at all owing to high viscosity. Asphaltene content from core extraction in the tar mat zone is typically around 60 wt%. [16] The immiscible region is not anticipated for the maltene-asphaltene binary mixture without gas in reservoirs. Therefore, it is plausible to choose a negative l23 = l32 so that the immiscible region convers a more realistic range. We use all these as a guide to assign l23 = l32 = -0.015 and l12 = l21 = 0.01, which are all within the range given by Prausnitz et al [24].

Fig. 5 Binodal and spinodal curves for three binary mixtures: (a) methane + maltene, (b) methane + asphaltene, and (c) maltene and asphaltene.

For the ternary (gas + asphaltene + maltene) mixture, the calculated ternary diagram at 339 K and 30 MPa is given in Fig. 6 using the same binary interaction parameters (l12 = l21 = 0.01, l13 = l31 = -0.01, and l23 = l32 = -0.015) as the previous diffusion simulation. In Fig. 6(a), the thick black dashed curve is the spinodal curve. Binodal (solid) curves include LLE curves on the top left (green-L1 and red-L2) and VLE curves close to the origin (blue-V, almost invisible) and on the bottom right (purple-L2). Tie-lines are also drawn with thin dashed lines (red-LLE, purpleVLE and light blue-VLLE). There is a three phase VLLE region on the bottom left inside the three open circles. In this region any compositions are separated into the three VLLE phases and the three phase compositions are located at the three open circles. According to the Gibbs phase rule, the degree of freedom for this ternary system is given by: F = C – P + 2 = 2 where C = 3 (3 components), P (3 phases). If temperature and pressure are specified, then, F = 0 (no freedom). That explains why the three phase equilibrium compositions are fixed at the reservoir temperature and pressure. Again, the stable single phase region is outside the binodal curves; the unstable region is inside the spinodal curves and the metastable region is in between the binodal and spinodal curves. It should be noted that l23 = l32 = -0.015 can avoid the immiscible region for binary maltene-asphaltene mixture at the specified reservoir condition, i.e., no spinodal and binodal curves close to the hypotenuse on the right-hand side, which we are not interested in at reservoir conditions because live fluid compositions cannot be in that area. As mentioned previously, the phase transition loci depend on the thermodynamic model and its parameters only. If the thermodynamic model and its parameters are fixed, the phase transition loci are fixed. However, compositional variation paths simulated by the diffusive model rely not only on the thermodynamic model and its parameters but also on the parameters of the diffusive model such as initial conditions, time and component diffusivity. 12


Page 13 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

The simulated compositional variation paths for the gas charge into the 100 m vertical reservoir over 1 million years in the previous section are also plotted in Fig. 6. It can be seen that the simulated compositional variation paths are very close to the top-left side owing to very low asphaltene mole fractions in the oil column. To see the results clearly, this area is enlarged in Fig. 6(b). The initial reservoir compositions are those top triangles and the GOC is behind the open circle almost on the y-axis. The simulated compositional variation paths of the diffusion with initial asphaltene mole fractions x2 = 0.0013 (2.9 wt%) and x2 = 0.0024 (5.2 wt%) are all in the stable region. However, the simulated compositional variation paths of the diffusion with x2 = 0.0028 (6.0 wt%) are very close to the metastable region whereas those with x2 = 0.0029 (6.3 wt%) cross the binodal and spinodal curves. In the case of high initial asphaltene content (x2 = 0.0029), asphaltenes start to concentrate [as shown Fig. 4(b) and (e)] near the GOC because of the aforementioned opposite effects of asphaltene diffusion and expulsion and then to precipitate there because the mixture is in the metastable region and quickly falls into the unstable region (across the spinodal boundary – dashed black curve). If asphaltenes accumulate in the upstructure of the reservoir like that, local asphaltene phase instability can occur. The field case [40] with a rapid gas charge is in the similar situation where the gas charge increases GOR and decreases oil solvency capacity to dissolve asphaltenes, thereby yielding asphaltenes instability and flocculation locally. It should be noted that the simulated compositional variation paths are almost identical at different times for a specified initial asphaltene content but the asphaltene accumulation (the largest asphaltene content) corresponding to the largest density in Fig. 4 is located at different depths as shown in Fig 7. In this 100-m vertical oil column, the asphaltene accumulation does not reach the base of the oil column, over 1 million years, where the impermeable diffusion boundary conditions are set. In addition, the simulated compositional variation path with the initial x2 = 0.0028 is much closer to the binodal curves compared to that with the initial x2 = 0.0024. Bigger density inversion is obtained with the initial x2 = 0.0028 as shown in Fig. 4. This means that large density inversion is obtained near asphaltene phase instability boundaries.

Fig. 6 Ternary diagram for gas + maltene + asphaltene mixtures at 339 K and l12 = 0.01, l13 = -0.01, and l23 = 0.015

13


Energy & Fuels

0.00

100

95

90 0.50 85

Height, m

0.25

Time, MY

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 19

0.75 80

1.00 0

0.001

0.002

0.003

0.004

0.005

75 0.006

Maximum asphaltene mole fraction Fig. 7 Variation of maximum asphaltene mole fraction with time and height for gas + maltene + asphaltene mixtures at 339 K and l12 = 0.01, l13 = -0.01, and l23 = -0.015 If simulation time increases from 1 million years to 8 million years (or decreasing thickness of the vertical oil column) with initial x2 = 0.0018 (4.0 wt%), the simulated results are also plotted in Fig. 6(b) (dotted purple curve). Because simulation time is much longer (or the thickness of the oil column is much smaller), gas can diffuse down to the base of the oil column as shown in Fig. 8. Gas addition lowers the oil solvency capacity to dissolve asphaltenes and expels asphaltenes out to the base. Therefore, asphaltenes are accumulated at the base of the oil column after 8 million years (or in the 35-m thick oil column after 1 million years). Asphaltenes at the base become unstable (spinodal decomposition) and the PDE solver has failed with x2 = 0.0028. Therefore, we lower the initial asphaltene mole fraction from 0.0028 (6.0 wt%) to 0.0018 (4.0 wt%) or shorten the gas charge time so that the simulation can be completed. It can be seen that compositional variation paths of the diffusion (dotted purple curve) in Fig. 6(b) cross the binodal curve near to the base of the oil column. The mixture is now in the metastable region. The simulated methane, asphaltene and density profiles are compared in Fig. 8 up to 8 million years. Gas reaches the base of the oil column and asphaltenes are concentrated at the base, thus causing asphaltene phase instability and forming a tar mat. Similar situation was observed in a Norwegian field case shows the similar situation. [41] In addition, as mentioned previously, density inversion induced advection significantly moves asphaltenes to the base. If enriched asphaltenes at the base of the oil column exceed the maximum asphaltene solubility in the oil, then a tar mat is formed, as observed in a Saudi field case study. [16]

Fig. 8 Simulated methane, asphaltene and density profiles with x2 = 0.0018 (4.0 wt%). (a) methane profiles; (b) asphaltene profiles; (c) density profiles.

14


Page 15 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

To analyze sensitivity of ternary phase diagrams to the lij’s, the calculated ternary diagram is illustrated in Fig. 9 with l12 = l21 = 0.002, l13 = l31 = l23 = l32 = -0.002. With this set of parameters, the LLE region is significantly increases compared to the ternary diagram in Fig. 6. The VLE region decreases whereas the VLLE region increases. Fig. 9(b) is the enlarged view of the interesting region from the initial reservoir compositions (top triangles) to the GOC (triangle on the y-axis). The simulated compositional variation paths of the diffusion are all in the metastable and unstable regions with the initial asphaltene mole fractions x2 = 0.0028 (6.0 wt% asphaltene, dashed-dotted blue curve) and x2 = 0.0013 (2.9 wt% asphaltene, dashed orange curve) whereas those are almost all in the metastable region with x2 = 0.0003 (0.7 wt% asphaltene, dotted light blue) in the 100-m thick vertical oil column over 1 million years. In the case of intermedium initial asphaltene content (x2 = 0.0013), the compositional variation paths of the diffusion are all in the metastable region and 0.5 and 1 million year curves are overlapped.

Fig. 9 Ternary diagram for gas + maltene + asphaltene mixtures at 339 K and 30 MPa with l12 = 0.002, l13 = 0.002, and l23 = -0.002.

Impacts of parameters on spinodal curves The impact of lij’s on phase diagrams is shown in the previous section. To a certain extent, this is equivalent to changing the other parameters of the thermodynamic model. The sensitivity of spinodal curves to other parameters of the thermodynamic model is also analyzed by changing one parameter at a time and keeping the others unchanged. Fig. 10 summarizes the results. From the previous discussions, we know that compositional variation paths simulated by the diffusive model are very close to top-left of ternary diagrams where spinodal and binodal curves are nearby. If spinodal curves are shifted toward the top-left region, opportunities of asphaltene phase instability increase. This can be used as a guide for selecting thermodynamic model parameters. The origin in Fig. 10 represents 100% gas (methane). Fig. 10(a) shows the impact of gas solubility parameters on spinodal curves. As anticipated, gas with lower solubility parameters slightly increases the unstable region, thus increasing a chance of asphaltene phase instability. The impact of asphaltene solubility parameters is given in Fig. 10(b). The higher asphaltene solubility parameters the bigger unstable region and the higher chance of asphaltene phase instability. Fig. 10(c) depicts the effect of maltene solubility parameters on spinodal curves. Lower maltene solubility parameters increase the unstable region and asphaltene phase instability. All these are easy to understand with the principle of “like dissolves like” in the solubility theory. [42, 43] Substances with similar solubility parameters can be mutually dissolved otherwise they are partially or completely immiscible and separated into different phases. Hence, poor solvents such as gas often lead to asphaltene phase instability like gas injection in enhanced oil recovery. The effect of temperature is shown in Fig. 10(d). An increase in temperature gives rise to a decrease in an unstable region and asphaltene phase instability. Fig. 10(e) shows the impact of asphaltene sizes on spinodal curves. The smaller asphaltene size the shorter and wider unstable region. Increasing molar mass of asphaltenes (nanoaggregates) increases chances of asphaltene phase instability because spinodal curves reach higher maltene and lower methane concentration regions where live oil compositions are often located. Finally, Fig. 10 (f) illustrates the

15


Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 19

impact of maltene molar mass on spinodal curves. Lighter oil slightly shifts the unstable region to the left-hand side, which is closer to live oil compositions and more easily to have asphaltene phase instability issues.

Fig. 10 Impacts of properties on ternary diagram for gas + maltene + asphaltene mixtures for l12 = 0.01, l13 = 0.01, and l23 = -0.01. (a) effect of gas solubility parameter; (b) effect of asphaltene solubility parameter; (c) effect of maltene solubility parameter; (d) effect of temperature; (e) effect of asphaltene molar mass; (f) effect of maltene molar mass.

Conclusions The simplified moving boundary 1-D diffusive model of Zuo et al. [17] is modified for gas charges into oil reservoirs coupled with the Maxwell-Stefan diffusivity. The thermodynamic nonideality is described by the FloryHuggins regular solution model. Reservoir live oil is grouped into three pseudocomponents: gas, asphaltene and maltene. The Peng-Robinson EOS fitted to the experimental PVT data is employed to calibrate the parameters of the three components in the Flory-Huggins regular solution model. The modified 1-D diffusive model is then utilized to simulate the gas charge process into the oil column. The fluid density inversion is created by the modified 1-D diffusive model because of opposite competing influences of asphaltene diffusion and expulsion. That is, the asphaltenes that are expelled from increasing solution gas at the top of the column build up somewhere behind the gas diffusion thereby increasing the density of this oil layer above that of the original oil. This density inversion then gives rise to advective or gravity currents. The gravity currents then migrate asphaltenes from top to base fairly rapidly in geological time. In addition, the same thermodynamic Flory-Huggins regular solution model calibrated by the Peng-Robinson EOS is used to generate multiphase spinodal and binodal loci by extending Michesen’s phase behavior algorithms. Then, the fluid distributions simulated by the modified 1-D diffusive model are compared with the phase diagrams. Asphaltenes may precipitate locally near the GOC for fluids with high asphaltene content such as shale breaks in upstructure where asphaltenes can be concentrated easily. This is because rapid gas addition to oil reservoirs easily 16


Page 17 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

makes the asphaltene enriched fluids near the GOC metastable and/or unstable (high gas solution oil lowers solvency capacity to dissolve asphaltenes). Furthermore, accumulation of asphaltenes at the base (OWC) by diffusion and advection can cause asphaltene phase instability and a tar mat formation when asphaltene content exceeds the maximum asphaltene solubility in the oil. This also leads to an increase in viscosity with depth since viscosity increases exponentially with asphaltene content. It is observed that large density inversion occurs near the asphaltene phase instability boundary, in particular in the metastable region in between spinodal and binodal boundaries. The sensitivity of the model parameters to phase envelopes is analyzed. Chances of asphaltene phase instability increases by decreasing gas and maltene solubility parameters, temperature, and maltene molar mass, and by increasing asphaltene solubility parameters and asphaltene molar mass.

References: [1] O’Keefe, M. Reservoir Fluid Properties Measured Downhole. JPT, 2009, August, 22–25. [2] Mullins, O. C. The modified Yen model. Energy & Fuels 2010, 24(4), 2179−2207. [3] Freed, D. E.; Mullins, O. C.; Zuo, J. Y. Theoretical Treatment of Asphaltenes in the Presence of GOR Gradients. Energy & Fuels 2010, 24(7), 3942–3949. [4] Zuo, J.Y.; Freed, D.; Mullins, O.C.; Zhang, D.; Gisolf, A. Interpretation of DFA Color Gradients in Oil Columns Using the Flory-Huggins Solubility Model. Paper SPE 130305 presented at the CPS/SPE International Oil & Gas Conference and Exhibition in China, Beijing, China, 8–10 June 2010. [5] Zuo, J. Y.; Mullins, O.C.; Freed, D.; Elshahawi, H.; Dong, C.; Seifert, D.J. Advances in the Flory−Huggins−Zuo Equation of State for Asphaltene Gradients and Formation Evaluation. Energy & Fuels 2013, 27, 1322−1335. [6] Juyal, P.; McKenna, A.M.; Yen, A.; Rodgers, R.P.; Reddy, C.M.; Nelson, R.L.; Andrews, A.B.; Atolia, E.; Allenson, S.J.; Mullins, O.C.; Marshall, A.G.; Analysis and identification of biomarkers and origin of blue color in an unusually blue crude oil. Energy & Fuels 2011, 25, 172−182. [7] Gisolf, A.; Dubost, F. X.; Zuo, Y. J.; Williams, S.; Kristoffersen, J.; Achourov, V.; Bisarah, A.; Mullins, O. C. Real Time Integration of Reservoir Modeling and Formation Testing. Paper SPE# 121275. Presented at the EUROPEC/EAGE Annual Conference and Exhibition, Amsterdam, 8–11 June, 2009. [8] Dong, C.; Petro, D.; Pomerantz, A.E.; Nelson, R.K.; Latifzai, A.S.; Nouvelle, X.; Zuo, J.Y.; Reddy, C.M.; Mullins, O.C. New Thermodynamic Modeling of Reservoir Crude Oil. Fuel 2014, 117, 839-850. [9] Betancourt, S. S.; Dubost, F. X.; Mullins, O. C.; Cribbs, M. E.; Creek, J. L.; Mathews, S. G. Predicting Downhole Fluid Analysis Logs to Investigate Reservoir Connectivity. Paper SPE IPTC 11488 presented at IPTC, Dubai, UAE, 4–6 December, 2007. [10] Pastor, W.; Garcia, G.; Zuo, J. Y.; Hulme, R.; Goddyn, X.; Mullins, O. C. Measurement and EOS Modeling of Large Compositional Gradients in Heavy Oils. SPWLA Paper presented at the SPWLA 53rd Annual Symposium in Cartagena, Colombia, June 16-20, 2012. [11] Elshahawi, H.; Hows, M.; Dong, C.; Venkataramanan, L.; Mullins, O.C.; McKinney, D.; Flannery, M.; Hashem, M., Paper SPE 109684 presented at the SPE Annual Technical Conference and Exhibition, Anaheim, California, 11–14 November, 2007. [12] Zuo, J. Y.; Elshahawi, H.; Mullins, O. C.; Dong, C.; Zhang, D.; Jia, N.; Zhao, H. Asphaltene Gradients and Tar Mat Formation in Reservoirs under Active Gas Charging. Fluid Phase Equilibria 2012 315, 91-98. [13] Jackson, R.; Zuo, J.Y.; Agarwal, A.; Herold, B.; Kumar, S.; De Santo, I.; Dumont, H.; Ayan, C.; Mullins, O.C.; Mapping and Modeling Large Viscosity and Asphaltene Variations in a Reservoir Undergoing Active Biodegradation, SPE, 170794, ATCE, 2014. [14] Zuo, J. Y.; Jackson, R.; Agarwal, A.; Herold, B.; Kumar, S.; De Santo, I.; Dumont, H.; Ayan, C.; Beardsell, M.; Mullins, O. C. Diffusion Model Coupled with the Flory-Huggins-Zuo EoS and Yen-Mullins Model Accounts for Large Viscosity and Asphaltene Variations in a Reservoir Undergoing Active Biodegradation. Energy & Fuels 2015, 29, 1447−1460. 17


Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 19

[15] Stainforth, J.G. New Insights into Reservoir Filling and Mixing Processes. In Cubit, J.M.; England, W.A.; Larter, S.; (Eds.), Understanding Petroleum Reservoirs: Toward an Integrated Reservoir Engineering and Geochemical Approach, Geological Society, London, Special Publication, 2004. [16] Mullins, O.C.; Zuo, J.Y.; Seifert, D.; Zeybek, M.; Clusters of Asphaltene Nanoaggregates Observed in Oilfield Reservoirs, Energy & Fuels 2013, 27, 1752–1761. [17] Zuo, J.Y.; Chen, Y.; Pan, S.; Wang, K.; Mullins, O.C. Investigation of Density Inversion Induced by Gas Charges into Oil Reservoirs Using Diffusion Equations. Energy 2016, 100, 199 – 216. [18] Pan, S.; Zuo, J.Y.; Wang, K.; Chen, Y.; Mullins, O.C. A Multicomponent Diffusion Model for Gas Charges into Oil Reservoirs. Fuel 2016, 180, 384 – 395. [19] Michelsen, M. L. The Isothermal Flash Problem. Part I. Stability. Fluid Phase Equilibria 1982, 9(1): 1–19. [20] Michelsen, M. L. Calculation of Phase Envelopes and Critical Points for Multicomponent Mixtures. Fluid Phase Equilibria 1980, 4(1-1): 1–10. [21] Peng, D.Y.; Robinson, D.B. A Rigorous Method for Predicting the Critical Properties of Multicomponent Systems from an Equation of State. AIChE J., 1977, 23, 137-144. [22] Harding, S.T.; Floudas, C.A. Phase Stability with Cubic Equations of State: Global Optimization Approach. AIChE J., 2000, 46, 1422 – 1440. [23] Hoteit H.; Firoozabadi, A. Simple Phase Stability-Testing Algorithm in the Reduction Method. AIChE J. 2006, 52, 2909 – 2902. [24] Prausnitz, J. M. Lichtenthaler, R. N.; Azevedo, E. G. D. Molecular Thermodynamics of Fluid Phase Equilibria, 3rd ed. Prentice Hall PTR, Upper Saddle River, NJ, 1999. [25] Krishna, R.; Wesselingh, J.A. The Maxwell-Stefan Approach to Mass Transfer. Chem. Eng. Sci. 1997, 52, 861– 911. [26] Buckley, J. S.; Wang, X.; Creek, J. L. Solubility of the Least Soluble Asphaltenes. Ch. 16, pp. 401– 438 in Mullins, O. C.; Sheu, E. Y.; Hammami, A.; Marshall, A. G. eds. Asphaltenes, Heavy Oils and Petroleomics, New York, New York, USA: Springer, 2007. [27] Wang, J. X.; Buckley, J. S. A Two-Component Model of the Onset of Asphaltene Flocculation in Crude Oils. Energy & Fuels 2001, 15, 1004–1012. [28] Wang, J. X.; Creek, J. L.; Buckley, J. S. Screening for Potential Asphaltene Problems. Paper SPE# 103137. Presented at the 2006 SPE Annual Technical Conferences and Exhibition, San Antonio, TX, USA, 24–27 September, 2006. [29] Creek, J. L.; Wang, J.; Buckley, J. S. Verification of Asphaltene-Instability-Trend (ASIST) Predictions for Low-Molecular-Weight Alkanes. SPE Production & Operations 2009, 5, 360–367. [30] Vargas, F. M.; Gonzalez, D. L; Creek, J. L.; Wang, J. X.; Buckley, J. Hirasaki, G.J.; Chapman, W.G. Development of a General Method for Modeling Asphaltene Stability. Energy & Fuels 2009, 23, 1147–1154. [31] Panuganti, S.R.; Vargas, F.M; Gonzalez, D.L.; Kurup, A.S; Chapman, W.G. PC-SAFT Characterization of Crude Oils and Modeling of Asphaltene Phase Behavior. Fuel 2012, 93, 658–669. [32] Zuo, J.Y.; Zhang, D. Plus Fraction Characterization and PVT Data Regression for Reservoir Fluids near Critical Conditions. SPE Paper 64520, presented at the SPE Asia Pacific Oil and Gas Conference and Exhibition, Brisbane, Australia, 16-18 October, 2000. [33] Peng, D.-Y.; Robinson, D. B. A New Two-Constant Equation of State. Industrial Engineering Chemistry Fundamentals 1976, 15, 59–64. [34] Zuo, J. Y.; Mullins, O. C.; Freed, D.; and Zhang, D. A Simple Relation between Densities and Solubility Parameters for Live Reservoir Fluids. Journal of Chem. Eng. Data 2010, 55(9), 2964–2969. [35] Krooss, B. M. Experimental Investigation of the Molecular Migration of C1-C6 Hydrocarbons: Kinetics of Hydrocarbon release from Source Rocks. Org. Geochem. 1988, 13(1–3), 513–523.

18


Page 19 of 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

[36] Sheikha, H.; Mehrotra, A. K.; Pooladi-Darvish, M. An Inverse Solution Methodology for Estimating the Diffusion Coefficient, of Gases in Athabasca Bitumen from Pressure-Decay Data. J. Petro. Sci. Eng. 2006, 53, 189–202. [37] Pisani, L. Simple Expression for the Tortuosity of Porous Media. Transp. Porous Med. 2011, 88, 193–203. [38] Chen, Y.; Wang, K.; Zuo, J.Y.; Mullins, O.C. Dynamics of Tar Mat Formation due to Asphaltenes Accumulation under Gas Charge in Reservoirs. OTC-25752-MS, presented at the Offshore Technology Conference, Houston, TX, USA, 4–7 May 2015. [39] Vella, D.; Huppert, H.E. Gravity Currents in a Porous Medium at an Inclined Plane. J. Fluid Mech. 2006, 555, 353–362. [40] Dumont, H.; Mishra, V.; Zuo, J.Y.; Mullins, O.C.; Permeable Tar Mat Formation within the Context of Novel Asphaltene Science, SPE 163292, KIPCE, Kuwait, 2012. [41] Achourov, V.; Pfeiffer, T.; Kollien, T.; Betancourt, S.S.; Zuo, J.Y.; di Primio, R.; Mullins, O.C. Gas Diffusion into Oil, Reservoir Baffling and Tar Mats Analyzed by Downhole Fluid Analysis, Pressure Transients, Core Extracts and DSTs. Petrophysics 2015, 56(4), 346–357. [42] Hansen, C., Hansen Solubility Parameters: A User's Handbook, Second Edition, Boca Raton, Fla: CRC Press, 2007. ISBN 978-0-8493-7248-3 [43] Mullins, O.C.; Pomerantz, A.E.; Andrews, A.B.; Zuo, J.Y. Asphaltenes Explained for the Nonchemist. Petrophysics 2015, 56(3), 266–275.

19


Analysis of Asphaltene Instability Using Diffusive ... - ACS Publications

Recommend Documents