Asphaltene phase instability analysis in gas charges into oil reservoirs

A method includes placing a downhole acquisition tool in a wellbore in a geological formation containing a reservoir fluid. The method includes performing downhole fluid analysis using the downhole acquisition tool to determine at least one measurement of the reservoir fluid. The method includes using a processor to estimate at least one fluid component property by using an equation of state based at least in part on the at least one measurement of the reservoir fluid and to simulate a diffusion process using a diffusion model that takes into account the at least one estimated fluid property to generate a composition path. The method includes using a processor to estimate one or more phase envelopes based in part on the at least one fluid property and compare the one or more phase envelopes with the composition path. The method includes outputting a visualization identify potential areas of asphaltene instability.

BACKGROUND

This disclosure relates to determining one or more dynamic processes for a reservoir in a geological formation occurring over geological time.

Reservoir fluid analysis may be used to better understand a hydrocarbon reservoir in a geological formation. Indeed, reservoir fluid analysis may be used to measure and model fluid properties within the reservoir to determine a quantity and/or quality of formation fluids—such as liquid and/or gas hydrocarbons, condensates, drilling muds, and so forth—that may provide much useful information about the reservoir, such as where areas of asphaltene instability occur. This may allow operators to better assess the economic value of the reservoir, obtain reservoir development plans, and identify hydrocarbon production concerns for the reservoir. Numerous possible reservoir models may be used to describe the reservoir. For a given reservoir, however, different possible reservoir models may have varying degrees of accuracy. The accuracy of the reservoir model may impact plans for future well operations, such as enhanced oil recovery, logging operations, and dynamic formation analyses. As such, the more accurate the reservoir model, the greater the likely value of future well operations to the operators producing hydrocarbons from the reservoir.

SUMMARY

This summary is provided to introduce a selection of concepts that are further described below in the detailed description. This summary is not intended to identify key or essential features of the subject matter described herein, nor is it intended to be used as an aid in limiting the scope of the subject matter described herein. Indeed, this disclosure may encompass a variety of aspects that may not be set forth below.

In one example, a method includes placing a downhole acquisition tool in a wellbore in a geological formation, where the wellbore or the geological formation, or both, contain a reservoir fluid. The method includes performing downhole fluid analysis using the downhole acquisition tool in the wellbore to determine at least one measurement associated with the reservoir fluid. The method includes using a processor to estimate at least one fluid component property by using an equation of state based at least in part on the at least one measurement associated with the reservoir fluid and to simulate a diffusion process using a diffusion model that takes into account the at least one estimated fluid property to generate a composition path. The method includes using a processor to estimate one or more phase envelopes based at least in part on the at least one fluid property and compare the one or more phase envelopes with the composition path. The method includes using a processor to output a visualization based at least in part on the comparison of the phase envelopes with the composition path to identify potential areas of asphaltene instability.

In another example, one or more tangible, non-transitory, machine-readable media include instructions to receive at least one measurement representative of a portion of a reservoir fluid as analyzed by a data acquisition tool in a wellbore in a geological formation within a hydrocarbon reservoir. The instructions estimate at least one fluid component property by using an equation of state. The instructions simulate a diffusion process by executing a diffusion model that takes into account the at least one estimated fluid property to generate a composition path. The instructions estimate one or more phase envelopes based at least in part on the at least one fluid property. The instructions compare the one or more phase envelopes with the composition path, and output a visualization based at least in part on the comparison of the phase envelopes with the output to identify potential areas of asphaltene instability.

In another example, a system a downhole acquisition tool housing including a plurality of sensors configured to measure at least one fluid property of a reservoir fluid within a geological formation of a hydrocarbon reservoir; and a data processing system configured to predict areas of potential asphaltene instability from a visualization that depends at least in part on an instability analysis associated with the at least one measured fluid property; wherein the data processing system comprises one or more tangible, non-transitory, machine-readable media comprising instructions. The instructions estimate at least one fluid component property by using a suitable equation of state, and simulate a diffusion process using a diffusion model that takes into account the at least one estimated fluid property to generate a composition path. The instructions estimate one or more phase envelopes based at least in part on the at least one fluid property. The instructions compare the one or more phase envelopes with the composition path. The instructions output a visualization based at least in part on the comparison of the phase envelopes with the output.

DETAILED DESCRIPTION

The present disclosure relates to systems and methods for reservoir characterization, such as performing asphaltene phase analysis to identify areas of instability in reservoirs. Petroleum includes a complex mixture of hydrocarbons of various molecular weights, plus other organic compounds. The molecular composition of petroleum varies widely from formation to formation. The proportion of hydrocarbons in the mixture is highly variable and ranges from as much as 97 percent by weight in the lighter oils to as little as 50 percent in the heavier oils and bitumens. The hydrocarbons in petroleum are mostly alkanes (linear or branched), cycloalkanes, aromatic hydrocarbons, or more complicated chemicals like asphaltenes. The other organic compounds in petroleum may contain carbon dioxide (CO2), nitrogen, hydrogen sulfide (H2S), oxygen, and sulfur, and trace amounts of metals such as iron, nickel, copper, and vanadium.

Acquisition and analysis representative of formation fluids downhole in delayed or real time may be used in reservoir modeling. A reservoir model based on downhole fluid analysis may predict or explain reservoir characteristics such as, but not limited to, connectivity, productivity, lifecycle stages, type and timing of hydrocarbon, hydrocarbon contamination, reservoir fluid dynamics, composition, and phase. Over the life of the reservoir, reservoir fluids such as oil, gas, condensates may behave dynamically in the reservoir. The reservoir fluids may flow (e.g., diffuse) into and out of the reservoir and/or biodegrade. This may result in spatial variations in the reservoir fluids throughout the reservoir, which may appear as fluid gradients in the composition characteristics of the reservoir fluids. For example, a concentration of compositional components of the reservoir fluid (e.g., gas, condensates, asphaltenes, etc.) may or may not vary along a vertical depth of the reservoir.

FIG. 1depicts an example of a hydrocarbon reservoir8that may use a wireline downhole tool10to employ the systems and techniques described herein. The downhole tool10is inserted into a wellbore12located within the hydrocarbon reservoir8. For example, in the illustrated embodiment, the downhole tool10is suspended in the wellbore12from the lower end of a multi-conductor cable14that is spooled on a winch at surface16. The cable14is communicatively coupled to an electronics and processing system18. The downhole tool10includes an elongated body20that houses modules24,26,28,30, and32, that provide various functionalities including fluid sampling, fluid testing, operational control, and communication, among others. For example, the modules24and26may provide additional functionality such as fluid analysis, resistivity measurements, operational control, communications, coring, and/or imaging, among others. As should be noted, other data acquisition tools may be used to employ the systems and techniques described herein. For example, in certain embodiments, the data acquisition tool may be a logging while drilling tool, seismic acquisition tool, or the like.

In the example shown inFIG. 1, the module28is a fluid communication module28that has a selectively extendable probe36and backup pistons38that are arranged on opposite sides of the elongated body20. The extendable probe36selectively seals off or isolates selected portions of a wellbore wall40of the wellbore12to fluidly couple to the adjacent formation42and/or to draw fluid samples (e.g., reservoir fluid48) from the formation42. The probe36may include a single inlet or multiple inlets designed for guarded or focused sampling. The reservoir fluid48may be expelled to the wellbore through a port in the body20or the reservoir fluid48may be sent to one or more fluid sampling modules30and32. The fluid sampling modules30and32may include sample chambers that store the formation fluid. In the illustrated example, the electronics and processing system18and/or a downhole control system control the extendable probe assembly36and/or the drawing of a fluid sample from the formation42.

The modules26,24,28,30,32include sensors that may collect and transmit data50associated with the fluid properties and the composition of the reservoir fluid52to the electronics and processing system18at the surface16, where the data50may be stored and processed in a data processing system54of the electronics and processing system18. The data50may include a plurality of measurements representative of information associated with, for example, reservoir temperature, pressure, compositions, the gas-to-oil ratio (GOR), asphaltene content, density, viscosity, optical density, or a combination thereof of the reservoir fluid48.

The data processing system54may include a processor56, memory58, storage60, and/or display62. The memory58may include one or more tangible, non-transitory, machine readable media collectively storing one or more sets of instructions for operating the downhole tool10, calculating and estimating fluid properties of the reservoir fluid48, modeling the fluid behaviors using, for example, the equation of state models, and/or identifying certain reservoir realization scenarios associated with observed fluid behaviors. The fluid properties calculated and/or estimated by the equation of state models may include partial molar volume, molar mass, partial density, initial conditions, boundary conditions, solubility parameters, or a combination thereof. The memory58may simulate diffusion processes by utilizing a diffusion model to generate a diffusion composition path (e.g., output) at varying depths. The memory58may store comparison rules for comparing phase envelopes with the generated diffusion composition path and/or output. The memory58may store mixing rules associated with compositional characteristics of the reservoir fluid48, equation of state (EOS) models for equilibrium and dynamic fluid behaviors, reservoir realization scenarios (e.g., biodegradation, gas/condensate charge into oil, CO2charge into oil, fault block migration/subsidence, convective currents, among others), and any other information that may be used to determine the effects of the composition and fluid behaviors of the reservoir fluid48on reservoir productivity, and provide information about reservoir connectivity. The memory58may store instructions to generate a user visualization based at least in part on a comparison between the phase envelopes and the diffusion composition path. The data processing system54may use the fluid property and composition information of the data50to identify the plausible realization scenarios for the reservoir8, as discussed in further detail below. In certain embodiments, the data processing system54may apply filters to remove noise from the data50.

To process the data50, the processor56may execute instructions stored in the memory58and/or storage60. For example, the instructions may cause the processor to estimate fluid and compositional parameters of the reservoir fluid48, reservoir fluid gradients, diffusion composition paths, equilibrium and/or dynamic state of the reservoir fluid48, and realization scenarios associated with the fluid and compositional parameters and diffusion composition paths. As such, the memory58and/or storage60of the data processing system54may be any suitable article of manufacture that can store the instructions. By way of example, the memory58and/or the storage60may be read-only memory (ROM), random-access memory (RAM), flash memory, an optical storage medium, or a hard disk drive. The display62may be any suitable electronic display that can display information (e.g., logs, tables, cross-plots, reservoir maps, etc.) relating to properties of the well as measured by the downhole tool10and the plausible realization scenarios, including diffusion processes, associated with the reservoir8. It should be appreciated that, although the data processing system54is shown by way of example as being located at the surface16, the data processing system54may be located in the downhole tool10. In such embodiments, some of the data50may be processed and stored downhole (e.g., within the wellbore12), while some of the data50may be sent to the surface16(e.g., in real time). In certain embodiments, the data processing system54may use information obtained from petroleum system modeling operations, ad hoc assertions from the operator, empirical historical data (e.g., case study reservoir data) in combination with or lieu of the data50to determine plausible realization scenarios, including diffusion processes, occurring within the reservoir8.

As discussed above, the data50from the downhole tool10may be analyzed with the equation of state (EOS) models to determine how gradients and diffusion processes in reservoir fluid compositions are affected by various dynamic processes occurring within the reservoir8. The dynamic processes for the reservoir8may include gas/condensate charge, biodegradation, convective currents, fault block migration, and subsidence, among others. Gas charges into the reservoir8may result in asphaltene phase instability because addition of gas decreases the oil solubility parameters and lowers solvency capacity of the oil to dissolve asphaltenes.FIG. 2illustrates an embodiment of a realization scenario that may occur within the reservoir8. Moving from left to right, the diagram inFIG. 2illustrates the reservoir8saturated with immature oil82(e.g., black oil) and charged with gas84over time86. The immature oil82, also known as heavy/black oil, generally has a high concentration of high molecular weight hydrocarbons (e.g., asphaltenes, resins, C60+) compared to mature oil (e.g., light oil, gas), which has high concentrations of low molecular weight aliphatic hydrocarbons (e.g., methane (CH4), ethane (C2H6), propane (C3H8), C4, C5, C6+etc.). The longer the reservoir fluid (e.g., the reservoir fluid48) is within the formation42, certain high molecular weight hydrocarbons found in the immature oil82may breakdown into the low molecular weight aliphatic hydrocarbons that make up the mature/light oil. Additionally, over time, source rock (e.g., portion of formation42having hydrocarbon reserve) may be buried under several layers of sediment. As the sediment layers increase, a depth88of the source rock, reservoir temperature, and reservoir pressure also increase. The increased temperatures and pressures favor the generation of light hydrocarbons, which may enter the reservoir.

Over time, the low molecular weight aliphatic hydrocarbons (e.g., gas84) may be expelled from the source rock and travel through a high-permeability streak in the formation to the top of the reservoir unit. As shown in the middle diagram inFIG. 2, the gas84diffuses down into the reservoir8from top90to bottom92, thereby charging the immature oil82with the gas84. Late charge of gas84(e.g., diffusion of gas after the reservoir8has been saturated with immature oil) into the immature oil82destabilizes the reservoir8, resulting in a fluid gradient for several fluid properties of the immature oil82. For example, the late charge of gas84may cause fluid gradients in API gravity, gas-to-oil ratio (GOR), saturation pressure (Psat), and combinations thereof of the immature oil82. As shown in the middle diagram inFIG. 2, a GOR toward the top90is higher compared to a GOR toward the bottom92. In addition, asphaltenes94are generally insoluble in the gas84. Therefore, increased concentration of the gas84toward the top90of the reservoir8may cause the asphaltenes94to phase separate. Alternatively, the asphaltenes94may diffuse ahead of the gas front, and flow towards the bottom92(e.g., when the asphaltenes do not phase separate). Diffusion of the asphaltenes94ahead of the gas front may yield mass density inversions and gravity currents (convective currents), which may result in bitumen deposition upstructure and/or tar mats96at the bottom92of the reservoir8. For example, a flow of asphaltenes94to the bottom92may lead to a low concentration of asphaltenes94toward the top90compared to a concentration of asphaltenes94toward the bottom92, resulting in a concentration gradient for the asphaltenes94in the reservoir8.

FIG. 3illustrates an analysis method110for identifying possible areas of tar formation in accordance with an embodiment of the present techniques disclosed herein. A first step includes inputting (block112) data from the downhole fluid analysis tool and laboratory measurements including reservoir temperature, pressure, compositions, and asphaltene content. A second step includes estimating (block114) one or more basic fluid component properties such as partial molar volume, molar mass, partial density, initial conditions, boundary conditions, and solubility parameters. A third step includes simulating (block116) diffusion processes at varying depths, locations, and geologic times. A fourth step includes generating (block118) spinodal and binodal curves (including phase envelopes) in accordance with the techniques disclosed herein. A fifth step includes plotting (step120) phase envelope diagrams and diffusion composition paths for comparison. A sixth step includes identifying (step122) regions of asphaltene instability where possible tar mat formation occurs. The regions of instability may be represented by a visualization that is output for visual analysis.

A 1D diffusive model was derived for an N-component mixture, as described in U.S. Patent Application Ser. No. 62/234,433, entitled “Tar Mat Formation Prediction in Late-Charge Reservoirs,” filed on Sep. 29, 2015, which is hereby incorporated by reference in its entirety. The governing equation (N−1 independent variable) is expressed as:

where C and x are the vectors of the molarity (molar concentration) and mole fraction, t and Ctdenote the time and total molarity. The matrix [B] of the drag effects is given by the following expression of the Maxwell-Stefan diffusivities (Dijfor the i-j pair diffusivity):

The thermodynamic non-ideality of the mixture is computed by using the derivatives of the activity coefficients, which will be described in further detail below.

where γiand Δijare the activity coefficient of component i and the Kronecker delta function for the i-j component pair. As described herein, molarity can be converted to mole fraction by use of the following expression:

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

Initial conditions are set to be homogeneous compositions in the one dimension oil column, i.e.,
Ci(0,z)=Ci0, i=1, 2, . . . ,N(6)

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

A ternary mixture is taken as an example: gas (1)+asphaltene (2)+maltene (3), where components 1, 2 and 3 represent the gas, asphaltene, and maltene components, respectively. In the gas cap, pure methane is assumed (a pseudo-gas component can be presumed as well). It may be assumed that the gas component is in equilibrium with oil at the gas/oil contact (GOC) and no asphaltenes move up from the oil column to the gas cap because the gas component does not have solvency capacity to dissolve asphaltenes.

Changes to the GOC occur because of the addition of the gas to the oil column resulting in swelling of the column. The changes to the GOC can be estimated by the following equation:

where z0is the initial depth of the GOC. J1and C1are the gas flux and molar concentration at the GOC. G and O denote the gas and oil sides.

Initial and boundary conditions can be set differently according to actual situations. By combining the equations above with initial and boundary conditions, composition variations with time in the oil column can be calculated numerically. The calculation is conducted for solving the partial differential equations mentioned above.

As discussed above, gas charges into the reservoir8may results in asphaltene phase instability because addition of gas decreases the oil solubility parameters and lowers solvency capacity of the oil to dissolve asphaltenes. Estimating phase envelopes (e.g., phase boundaries) of the mixture of the reservoir fluid may help ensure the asphaltenes are stabilized in oil, as discussed in detail below. Prior to estimating phase envelopes of the reservoir fluid mixture, other measurements associated with the reservoir fluid are taken. The measurements may include reservoir temperature, pressure, compositions, the gas-to-oil ratio (GOR), asphaltene content, density, viscosity, or a combination thereof.

The phase envelopes (e.g., phase boundaries) may include phase instability boundaries including spinodal and binodal boundaries. The spinodal boundary condition may be construed as the limit of local stability with respect to small fluctuations, which is defined by the condition that the second derivative of Gibbs free energy is zero. The binodal boundary condition may be construed as the limit of a global minimum energy equilibrium state of the system. The binodal and spinodal boundaries meet at a critical (plait) point. As discussed in detail below, the methods disclosed herein describe generation of both spinodal and binodal loci.

Spinodal Curve Calculations

As may be appreciated, the spinodal criterion for an N-component system at a specified temperature and pressure is that the determinant of the Jacobian matrix of molar Gibbs free energy is equal to zero:

where g is the molar Gibbs free energy that is a function of temperature (T), pressure (P) and mole fraction (x1, x2, . . . , xN). The molar Gibbs free energy is related to the chemical potential by:

where μkis the chemical potential of component k. Therefore:

According to the Gibbs-Duhem equation in thermodynamics, the last term of Eq. 12 is zero, and

Equation (10) may then be rewritten as:

The chemical potential is further calculated by an activity coefficient model or an equation of state (EOS):
μi=μi0+RT ln fi=μiref+RT ln xiγi=μi0+RT ln xiφiP(16)

where superscripts ref and 0 are the reference and standard states, fi, γiand φiare the fugacity, activity coefficient, and fugacity coefficient of component i. Thus, Eq. 15 can be expressed as:
det|[Γ′]|=0   (17)

where the matrix [Γ′] is associated with [Γ] in diffusion Equation (1) by:

Therefore, the spinodal instability can be determined directly by evaluating the determinant value of matrix [Γ′]. If det|[Γ′]|<0, the mixture is unstable. Otherwise, the mixture is stable. If an EOS model is used, the activity coefficients are replaced by the fugacity coefficients.

For a ternary mixture, there are three unknowns (x1, x2and x3) at a specified temperature and pressure, but there are two equations which include Equation (17) 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. The Newton method is used to determine spinodal curves at a specified temperature and pressure. As 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 Equation (17) and the summation equation are satisfied.

Binodal (e.g., coexistence equilibrium) curves can be determined in terms of phase equilibrium criteria. As may be appreciated, vapor-liquid-liquid three phase equilibrium (VLLE) criteria are given by:
fiV=fiL1=fiL2, i=1, 2, . . . , N   (19)

where fiis the fugacity of component i and superscripts V, L1 and L2 denote the vapor, liquid 1 and liquid 2 phases. For simplicity, if the non-idealities of the gas phase and two liquid phases are described by the activity coefficient model (e.g., the model parameters are calculated by the cubic equation of state at reservoir conditions), the phase equilibrium criteria can be rewritten as:
yiγiV=xiL1γiL1=xiL2γiL2, i=1, 2, . . . , N   (20)

where yi, and xiare the mole fractions of component i.

For specified reservoir fluid compositions (z1, z2, . . . , zN), a tangent plane distance algorithm is used for phase stability check and a flash algorithm is employed for VLE, LLE and VLLE calculations. Because of convergence in the vicinity of the critical point, Michelsen's phase envelope algorithm is extended to calculate liquid-liquid coexistence curves and plait (critical) points. For LLE, there are N+1 equations which are given by:
gi(α, β)=ln xiL1+ln γiL1−ln xiL2−ln γiL2=ln Ki+ln γiL1−ln γiL2=0,i=1, 2, . . . ,N(21)

where Kiis the equilibrium constant which is dependent on temperature, pressure and compositions. The (N+1)thequation is the Rachford-Rice equation given by:

where F is the mole fraction of the L1 phase. Adding one more specification equation:
gN+2(α, β)=αk−S=0   (23)

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

Let zmbe a specified variable and zjone of the primary variables, the remaining compositions are calculated by:

which yields the following equation for a ternary system:
zi=1−zm−zj, i≠j,m(27)

where superscript 0 stands for the initial global compositions of components i, j and m.

Because the Newton-Raphson method is used to solve the set of nonlinear equations, good initial guesses are useful. Therefore, the same extrapolation method is utilized as Michelsen. A polynomial extrapolation is used to obtain a good initial estimation for kthiteration:
αk(S)=αk0+αk1S+αk2S2+αk3S3(28)

where αkj, j=0, 1, . . . , 3, are the coefficients of the polynomial regressed from the information of the last two points. This extrapolation method utilizes two points at the beginning. Therefore, 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. Equation (28) is then used afterward. At plait (critical) points, Ki=1 for i=1, 2, . . . , N.

Activity Coefficient Model

The activity coefficient is calculated by the Flory-Huggins regular solution model:

where δ is the solubility parameter and ljkis the binary interaction parameter between components j and k. For pure component j, ljj=0. If ljk=0, the Equation (29) is reduced to:

In one example, the reservoir temperature and pressure are set to 339 K and 30 MPa. The density and partial molar volume of gas and maltenes are calculated by the Peng-Robinson equation of state. The saturated gas solubility in maltenes at the reservoir condition is calculated by the Peng-Robinson equation of state, which is equal to 0.575. The solubility parameters of gas and maltene components are calculated by:
δi=17.347 ρi+2.904   (34)

where ρiis the density of component i in g/cm3and δiis the solubility parameter of component i in MPa0.5. The solubility parameter of asphaltenes in MPa0.5is calculated by:
δ=δ0(1−1.07×10−2(T−298.15))   (35)

where δ0is 21.85 MPa0.5for asphaltenes and T is the temperature in K.

Asphaltenes are present in the oil as nanoaggregates for black oil. According to the Yen-Mullins model, 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.

It is assumed that the basic parameters of each components are constant at the specified reservoir pressure and temperature as given in Table 1. The 100-m oil column initially has homogeneous mole fractions for substantially all the components as shown in column 1 of Table 1, and the gas cap 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.

The methane migration diffusion coefficients in source rocks is estimated ˜1.62×10−11m2/s, and is used as diffusion coefficient for methane-maltene pair: D13=D31=1.62×10−11m2/s. A measured methane diffusion coefficient of (˜1.3×10−10m2/s) in Athabasca bitumen at 50-75° C. and 8 MPa. Due to tortuosity of porous media, diffusion coefficients in typical sandstones is ˜10 times lower than those in free volume. Therefore, the diffusion coefficient is estimated for methane-asphaltene pair as D12=D21=0.9×10−11m2/s in porous media. Because there is no data available for maltene-asphaltene pair in the open literature, the diffusion coefficient for maltene-asphaltene pair is estimated at D23=D32=0.5×10−11m2/s in porous media.

The diffusion model is at first used to simulate the gas charge process into the oil column. Adjusting l13=l31=−0.01 in the activity coefficient model to match saturated gas solubility in maltenes to be equal to 0.575, which is the same as the value obtained by the Peng-Robinson equation of state. l12=l21=0.01, and l23=l32=−0.015. It is observed that density inversion is generated due to the gas charge into the oil column. The extent of density inversion depends on initial asphaltene concentration. The higher asphaltene content the larger density inversion. If the initial asphaltene mole fraction is set to 0.0024 (5.4 wt %), the maximum density inversion is <1 kg/m3whereas if it is set to 0.0028 (6.0 wt %), the maximum density inversion is 3.0 kg/m3. Increasing initial asphaltene content further results in asphaltene instability, which will be discussed in 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, the fluid density increases rapidly (e.g., almost linearly) and reaches a maximum point where fluid density is slightly (Δρ≈3.0 kg/m3) higher than the initial fluid density. Then the fluid density gradually decreases below the initial density value and subsequently increases to the initial density value. Thus, fluid density inversion is created. The density inversion generated by the gas charge can induce density driven convention (e.g., gravity currents), thus moving asphaltenes down to the base of the oil column. The fluid density inversion over time scale can be seen inFIG. 4. It may be appreciated that asphaltene phase instability may occur during the gas charge process, which depends on conditions. Therefore, it is helpful to analyze asphaltene phase instability thermodynamically.

As shown inFIG. 4, the density inversion of the simulated fluid density variation with depth at different times with initial asphaltene mole fraction of 0.0024 is shown by plot126. As shown inFIG. 5, the density inversion of the simulated fluid density variation with depth at different times with initial asphaltene mole fraction of 0.0028 is shown by plot128. The larger density inversion can be observed by the larger asphaltene mole fraction as generated by the diffusion model.

Phase Instability Analysis

The phase diagrams of the three constitutive binary mixtures for the ternary system described herein are calculated using the methods described above. The calculated binodal and spinodal curves are shown inFIGS. 6-8. 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 phases inside the spinodal curves, which may be referred to as spinodal decompositions, where a new phase is formed spontaneously. It may be appreciated that the metastable region is inbetween the binodal and spinodal curves.

To check the sensitivity of phase diagrams to the binary interaction parameters in the activity coefficient model for the three binary mixtures, the phase diagrams with and without binary interaction parameters are generated for comparison as shown inFIG. 4. It can be seen that positive binary interaction parameters enlarge metastable and unstable regions whereas negative values decrease metastable and unstable regions.FIG. 6depicts the binodal and spinodal curves for three binary mixtures. For example, lines136and138illustrate spinodal curves for a methane and maltene mixture130, while lines132and134illustrate binodal curves for the methane and maltene mixture130.FIG. 7depicts lines148and144illustrate spinodal curves for a methane and asphaltene mixture140, while lines146and148illustrate binodal curves for the methane and asphaltene mixture140.FIG. 8depicts lines152and154illustrate spinodal curves for a maltene and asphaltene mixture150, while lines156and158illustrate binodal curves for the maltene and asphaltene mixture150.

It may be appreciated that lijbinary interaction parameter 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 lijis similar to familiar kijin cubic equations of state. Unfortunately, in most cases, the lijparameter 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 mentioned previously, l13=l31=−0.01 was obtained by matching the saturated gas solubility in maltenes at the GOC (x1=0.575).

It may be appreciated that maltene and asphaltene components can be completely miscible at any given compositions as shown in the mixture150when temperature above 320 K by changing l23=l32from 0 to −0.015. It can be seen that the immiscible region has very high asphaltene contents (x2>0.4 when T >320 K with zero l23, corresponding to >90 wt % asphaltenes). Oil with such high asphaltene concentration (>90 wt %) is not very mobile owing to high viscosity. Asphaltene content from core extraction in the tar mat zone is around 60 wt %. The immiscible region is not anticipated for maltenes-asphaltenes without gas in reservoirs. Therefore, it is plausible to set a negative l23=l32. As such, it is estimated that l23=l32=−0.015 and l12=l21=0.01.

For a ternary (gas+maltene+asphaltene) mixture160, the calculated ternary diagram is given inFIG. 9with l12=l21=0.002, l13=l31=l23=l32=−0.002. InFIG. 9, line165is the spinodal curve. The remaining solid curves are binodal curves. Lines164and166illustrate binodal curves, including liquid-liquid equilibrium (LLE) curves on the top left and lines168and170illustrate binodal curves, including vapor-liquid equilibrium (VLE) curves close to the origin and on the bottom right. Tie-lines172,174, and176are also drawn for the LLE region, the VLE region, and the VLLE region. There is a three phase vapor-liquid-liquid equilibrium region on the bottom left inside the three open circles. In this region any compositions are separated into three vapor-liquid-liquid phases and the three phase compositions are located at the three open circles (e.g.,178,180,182). 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 (e.g., 3 components), P (e.g., 3 phases). If temperature and pressure are specified, then, F=0 (i.e., no degrees of freedom). Thus, 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 inbetween the binodal and spinodal curves. It should be noted that l23=l32=−0.002 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.

FIG. 10is an enlarged view ofFIG. 9. It can be seen that the simulated compositions are very close to the top-left side owing to very low asphaltene mole fractions. Line162represents the spinodal curve. The simulated composition paths of the diffusion are in the metastable and unstable region for initial asphaltene mole fractions of 0.0028 (e.g., line169) and 0.0024 (e.g., line192) in the 100-m thick oil column over 1 million years. Asphaltenes start to precipitate at the GOC because the mixture is in the metastable region168and then quickly falls into the unstable region (across the spinodal boundary194).

To analyze sensitivity of ternary phase diagrams to the lij's, the calculated ternary diagram is illustrated inFIG. 11with l12=l21=0.01, l13=l31=−0.01, l23=l32=−0.015. With this set of parameters, the liquid-liquid equilibrium region is smaller compared to the ternary diagram inFIGS. 9-10mainly because of l23=l32=−0.015 from l23=l32=−0.002. The vapor-liquid equilibrium region202increases and the VLLE three phase equilibrium region decreases.FIG. 12is the enlarged view of the interesting region from the initial reservoir compositions to the GOC. The simulated compositional paths of the diffusion over 1 million years in the 100-m thick oil column with initial asphaltene mole fractions of 0.0028 (e.g., line216) and 0.0024 (e.g., line214) are above the binodal curve (e.g., line218). That means there are no asphaltene instability issues.

FIGS. 13-16illustrate methane and asphaltene profile in 100 m and 35 m thick oil columns. When thickness of the oil column is lowered from 100-m thick to 35-m thick, the same parameters are used for both cases. There are no asphaltene instability issues at the top of the oil columns (e.g., profiles220and240) as shown inFIGS. 13 and 15. However, because the thickness of the oil column is much smaller, gas can diffuse down to the base of the oil column. Gas charges lower the oil solvency capacity to dissolve asphaltenes and expel asphaltenes out to the base. Therefore, asphaltenes are accumulated at the base of the 35-m thick oil column after 1 million years. As such, asphaltenes at the base are unstable (e.g., due to spinodal decomposition).

When the initial asphaltene mole fraction from 0.0028 (6.0 wt %) to 0.0018 (4.0 wt %) or the gas charge time for the 35-m thick oil column is shortened, the simulation can be completed. It can be seen that compositional path of the diffusion (e.g., profiles230and250) are below the binodal curve close to the base of the 35 m oil column as shown inFIGS. 14 and 15. The mixtures are in the metastable region. The simulated asphaltene profiles for both 100-m (e.g., profile240) and 35-m (e.g., profile250) thick oil columns are shown inFIGS. 15 and 16. The gas front does not reach the base in the 100-m oil column, whereas the gas front goes through the base and causes asphaltene instability, thus yielding a tar mat formation. As mentioned previously, density inversion induced convection 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.

Impacts of Parameters on Spinodal Curves

The impact of changing the values of lij's on the phase diagrams may be further understood with reference toFIGS. 17-22. To some extent, changing the values of the lij's on the phase diagrams may be equivalent to changing other parameters in the model. The sensitivity analysis is conducted by changing one parameter at a time, while the other parameters are kept unchanged.FIG. 17shows the impact of gas solubility parameters on spinodal curves. Gas with lower solubility parameters slightly increases the unstable region. The impact of asphaltene solubility parameters is illustrated inFIG. 18. It can be seen that higher asphaltene solubility parameters result in a bigger unstable region.FIG. 19depicts the effect of maltene solubility parameters on spinodal curves. It can be seen that lower maltene solubility parameters increase the unstable region. This may be further understood as the solubility theory in chemistry applies. In other words, substances with similar solubility parameters can be mutually dissolved. Otherwise they are partially or completely immiscible and separated into different phases.FIG. 20depicts the effect of temperature. An increase in temperature gives rise to a decrease unstable region.FIG. 21shows the impact of asphaltene sizes on spinodal curves. Asphaltene molecules with molar mass of ˜750 g/mol cause a large spinodal region and increase methane concentration in oil. Thus, generally no asphaltene instability issues exist because reservoir conditions are often far away from the unstable region. Increasing molar mass of asphaltenes (e.g., nanoaggregates) increases chances to have asphaltene instability because spinodal curves reach higher maltene and lower methane concentration regions.FIG. 22illustrates the impact of maltene molar mass on spinodal curves. Lighter oil slightly shifts the unstable region to the left-hand side.