Enhancing asphaltene precipitation modeling by cubic-PR solid model using thermodynamic correlations and averaging techniques

Cubic equation-of-state solid models are one of the most widely used models to predict asphaltene precipitation behavior. Thermodynamic parameters are needed to model precipitation under different pressures and temperatures and are usually obtained through tuning with multi asphaltene onset experiments. For the purpose of enhancing the cubic Peng–Robinson solid model and reducing its dependency on asphaltene experiments, this paper tests the use of aromatics and waxes correlations to obtain these thermodynamic parameters. In addition, weighted averages between both correlations are introduced. The averaging is based on reported saturates, aromatics, resins, asphaltene (SARA) fractions, and wax content. All the methods are tested on four oil samples, with previously published data, covering precipitation and onset experiments. The proposed wax-asphaltene average showed the best match with experimental data, followed by a SARA-weighted average. This new addition enhances the model predictability and agrees with the general molecular structure of asphaltene molecules.


Introduction
As Kokal and Sayegh described asphaltenes as "The cholesterol of petroleum" (Kokal and Sayegh 1995), asphaltene precipitation/deposition is recognized as one of the most significant problems in the oil industry. Asphaltene precipitation/deposition could occur in the reservoir, perforation tunnels, wellbore and tubulars, surface chokes, flowlines, processing facilities (separators), transmission pipelines, and stock tanks. Precipitation in any stage of the previous stages causes interruption, reduction, and possible stoppage of the production process, letting aside the costs of remedial/ removal actions. To avoid such problems, it is best to mitigate asphaltene precipitation either by avoiding the conditions of instability, or by adding inhibitors to prevent the precipitation process.
Asphaltenes are a fraction of four fractions more commonly known to the industry as SARA fractions. SARA analysis is usually performed on stock tank oils to report percentage, in wt% or mol%, of each fraction. Wax content may also be reported, usually in wt%.
Screening criteria have been proposed to check whether the risk of asphaltene precipitation exists or not (Wang et al. 2006;Shokrlu et al. 2011;Ahmed 2013). If there is potential of asphaltene precipitation, laboratory experiments are often conducted to quantitate this precipitation and/or to detect the conditions under which the precipitation initiates. These experiments are costly and require some relatively special setup (e.g., Solid Detection System or SDS equipment) that may not be found in every PVT laboratory.

Asphaltene stability
Asphaltene onset pressure (AOP) represents the pressure at which asphaltene precipitation occurs, keeping temperature and composition constant. Figure 1 shows the asphaltene precipitation envelope (APE) on a typical phase diagram of oil.
Experiments are usually performed to measure the AOP values at different temperatures, representing points lying on the upper asphaltene envelope. Other experiments aim to measure the quantity of precipitate (PPT), usually in wt%, versus (vs.) pressure at constant temperature and composition. These experiments simulate depletion or overall production and are sometimes repeated at different temperatures.

Asphaltene precipitation models
Modeling is an important step to track precipitation conditions and the amount of precipitated asphaltenes as the oil moves from the reservoir through the production system. These values are used in the mitigation/inhibition of asphaltene precipitation. Modeling is coupled with the performed experiments to ensure reliability. Several models have been proposed in the literature to predict the behavior of the asphaltene precipitation process (Subramanian et al. 2016). These models could be classified based on how the precipitation process is approached.
Thermodynamic calculations approach is more common in modeling asphaltenes. In this approach, thermodynamic equilibrium calculations for all the components and/or pseudo-components of the system are performed to predict the phase behavior of oil at specified conditions. This approach consists of two main procedures that differ in how the asphaltene component is 1 3 treated in the thermodynamic equilibrium calculations. The two procedures are defined by the solubility or colloidal approaches. Solubility approach is based on the lyophilic theory, which states that the hydrocarbon fluid holds asphaltenes dissolved. Precipitation occurs if this dissolving power drops below a critical level. On the other hand, colloidal approach represents the lyophobic theory, which assumes asphaltenes as insoluble in oil. However, resins can stabilize asphaltene particles by adsorbing onto their surface, preventing them from aggregation.
The solubility approach deals with quantifying asphaltene solubility, by calculating either the solubility parameters of asphaltenes and crude oil, or the binary interaction coefficients (BICs) between asphaltene and the other oil constituents (Subramanian et al. 2016). Consequently, there are two types of modes underneath the solubility approach: solubility parameter mode and equation-of-state (EOS) mode. Different models are included under the solubility parameter mode (Hirschberg et al. 1984;Thomas et al. 1992;Cimino et al. 1995;Chung 1992;Buckley et al. 1998;Yarranton and Masliyah 1996). These models deal with asphaltene as a single component. A fewer number of models (Mansoori et al. 1988;Kawanaka et al. 1991) have been proposed to consider the heterogeneous nature of asphaltene mixtures.
The second mode is the EOS mode. Only EOSs and, sometimes, some additional distinct calculations are needed to model asphaltene precipitation behavior. EOS mode consists of three schemes or three types of EOSs: Statistical Associating Fluid Theory (SAFT) EOSs (Ting et al. 2003;Gonzalez et al. 2005;Panuganti et al. 2012), Cubic Plus Association (CPA) EOSs (Du and Zhang 2004;Sabbagh et al. 2006;Li and Firoozabadi 2010;Shirani et al. 2012), and cubic EOSs (Thomas et al. 1992;Nghiem et al. 1993;Nghiem and Coombe 1997;Gupta 1986;Pedersen and Christensen 2007;Kohse et al. 2000). SAFT and cubic EOS models are often used to model asphaltene precipitation behavior (Abouie et al. 2016a, b), while EOS models (in general) are considered the most widely used (Subramanian et al. 2016). These models assume the asphaltene is treated as one homogeneous component.
The objective of this study is to enhance the PR-cubic solid model of asphaltenes, decreasing its dependence on asphaltene experiments and improving its predictability.

Model equations
The model used in this work is the cubic solid model developed by Ngheim et al. (1993), Nghiem and Coombe (1997) and Kohse et al. (2000). The model has several assumptions that can be summarized as follows (Subramanian et al. 2016): • Homogenous asphaltene behavior • Molecular interactions are not accounted for • Reversible asphaltene precipitation, although this assumption can be modified according to the procedure proposed by Nghiem et al. (2001) and Kohse and Nghiem (2004) • Asphaltene molecular structure or geometry is not accounted for The relation used to estimate the solid fugacity is given by: where the '*' superscript indicates reference conditions, v a is the molar volume of asphaltene at reference conditions, R is the universal gas constant, subscript "tp" indicates triplepoint conditions, ΔH is the enthalpy, and ΔC P the isobaric heat capacity difference between liquid and solid.
For isothermal processes (T = T*), the above equation reduces to:

Thermodynamic parameters
Triple-point temperature and pressure, enthalpy, and heat capacity are thermodynamic parameters of the modeled solid (asphaltenes). For asphaltenes, the triple-point pressure is considered atmospheric. The triple-point temperature is taken as the fusion temperature. In the application of the model, the waxes correlation of fusion temperature is the only one that is really used (Eq. 3). Enthalpy of fusion (ΔH f ) and heat capacity difference (ΔC P ) are used as matching parameters to adjust the APE using the measured AOPs Kohse et al. 2000;Tavakkoli et al. 2014).
where MW is the molecular weight and the subscript 's' stands for solid (asphaltene).

Fluid characterization
Splitting the reported plus fraction of the sample into single carbon numbers (SCNs), assigning critical properties to the obtained SCNs, and lumping the pseudo-components are done following the methodology mentioned in Nghiem et al.
Binary interaction coefficients (BICs) of hydrocarbon-hydrocarbon, excluding Asphaltene-hydrocarbon, are given by the following formula (Li et al. 1985): where v c and e are critical volume and tuning exponent parameter, respectively. Finally, the Asphaltene-hydrocarbon BICs are considered a single value for the light components symbolized as δ asp and zero for the heavy components up to the non-precipitating component (C 31A+ ). It should be noted that δ asp is greater than the corresponding BIC values of the non-precipitating component, because asphaltene is less homogenous with the light components than the nonprecipitating component. For the non-hydrocarbon BICs, the values reported in the SPE Monograph (Whitson and Brule 2000), correlations mentioned by Ahmed (2013), or any appropriate values could be used.
Tuning is performed to match the saturation pressure and available asphaltene experiments. The tuning parameters are the exponent parameter e, Asphaltene-light hydrocarbon BIC δ asp , and asphaltene molar volume v a . Both e and δ asp affect the saturation pressure; therefore, they should be tuned simultaneously.
Extra information about fluid characterization and flash calculations is found in supplementary material. Shoukry et al. (2019) proposed to use waxes (Won 1996) and aromatics (Pan et al. 1997) correlations for T f and ΔH f (Eqs. 6-8), along with the ΔC P correlation developed by Pedersen et al. (1991) (Eq. 9). The use of thermodynamic

Proposed modifications
correlations quickens the modeling process. The waxes correlation of enthalpy is listed as: where the aromatics correlations are given by: and the heat capacity correlation is presented as: The heat capacity difference correlation is the same for waxes and aromatics (Schlumberger 2009).
This study computes weighted averaging methods of the parameters, calculated using both waxes and aromatics correlations. These averages depend on the SARA analysis of the samples and the reported wax content (if available). The following simple equation along with Table 1 demonstrates the proposed averaging techniques.
where θ is the fusion temperature and/or the enthalpy of fusion.
The above modification seeks the optimum combination of already available correlations, so that the best match to the measured laboratory values can be achieved.
The different proposed averages do not represent/change the original reported SARA fractions, nor they interfere with the sample's composition. They are just used in computing the thermodynamic parameters (fusion temperature and enthalpy of fusion) required for Eq. 1.

Fluid samples
Four oil samples are tested in this study, covering two types of experiments (AOP and PPT wt% vs. pressure). All the samples are taken from published papers/researches. The samples were chosen due to: (1) having enough data for modeling/calculations and (2) covering a good span of API and asphaltene content. Table 2 summarizes the reservoir
The reference conditions for samples with AOP experiment are taken at the AOP with the lowest temperature. In case of PPT wt% versus pressure set of experiments, the reference pressure is taken as the highest pressure point of the set.

PPT wt% versus Pressure (Sample 1)
A total final of 12-component characterization is performed on Sample 1 using the PPT wt% versus pressure at one temperature. Then, this characterization is used to predict the PPT wt% versus pressure at another temperature. The tuning of this sample depends on the PPT wt% versus pressure data measured at 369.2 K. This is an isothermal process, so thermodynamic parameters have no effect. Table 4 shows the lumped components' compositions along with their molecular weights. The obtained characterization of this sample is then used to fully predict the PPT wt% versus pressure curve at 386.2 K, so that the experimental values reported at 386.2 K are only used for error calculation from the predicted ones.

AOP (Samples 2, 3, and 4)
The three samples (2, 3, and 4) have been fully characterized with the result output characterization reported by Abouie et al. (2016a, b). The characterized samples are present in supplementary material. This characterization is used directly in the proposed model for predicting the upper asphaltene envelope and comparing it with the measured AOP values. It is worth mentioning that the reported characterization in the two papers was done using other asphaltene data (lower onset pressures). In addition, the characterization was done to ensure a consistent shape of precipitate versus pressure. Consequently, measured AOPs have not been utilized in the characterization procedure, and their prediction using the reported characterization is valid.

PPT wt% versus pressure (Sample 1)
The tuned values for e, δ asp (taken till the C 6 -C 7 component), and asphaltene molar volume v a are 0.56, 0.21, and 0.908 (m 3 /kg-mole), respectively. Figure 2 shows the tuning results for this sample. Only one curve is present as the tuning process is isothermal (thermodynamic parameters are irrelevant as shown in Eq. 2). The MW of asphaltenes for this sample is 656. Table 5 represents the different combinations of thermodynamic parameters (to be tested in Eq. 1) for prediction of PPT wt% versus pressure at the other temperature (386.2 K). Figure 3 shows the prediction results for Sample 1. The different curves represent different methods of obtaining the thermodynamic parameters (waxes correlations, aromatics correlations, and different averaging techniques). The curves are almost parallel, but starting from a different AOP according to the method of obtaining thermodynamic parameters. The predicted versus measured values are shown in Table 6.

AOP (Samples 2, 3, and 4)
Figures 4, 5, and 6 show the upper APEs modeled using the different correlations and averages for Samples 2 through 4, respectively. Similar to Fig. 3, different methods of computing the thermodynamic parameters result in different curves. Table 7 shows the calculated α and β values of the samples, in addition to the comparison between the obtained AOP values through modeling (using the waxes correlations, aromatics correlations, and the proposed averages) and laboratory measured values.

Discussion
The typical procedures of the cubic-PR solid model need two sets of experiments: one to obtain the standard model and characterization parameters (δ asp , e, and v a ) and the other (usually AOPs) to obtain thermodynamic parameters (ΔH f and ΔC P ) through tuning. Using the thermodynamic correlations eliminates the need of extra experimentations. Nevertheless, the proposed averages between waxes and aromatics correlations are shown to be superior in terms of matching the measured lab data.

3
The waxes and aromatics correlations were tested against an actual measured value of heat of fusion of asphaltene (Zhang et al. 2003;Gray et al. 2004). The measured value fell between the values obtained from the two correlations. These values are within the range of the values calculated for the samples in this study.
Using the aromatics correlations had variable results depending on the sample. Samples 3 and 4 exhibited better matching with the aromatics' correlations than waxes'. Almost all the experimental data show that using aromatics correlations appears to underestimate, in overall, values compared to measured ones. These underestimated values reach a limit that model's extrapolation into the two-phase envelope is made to check the value of AOP, which is unlikely to be the true case. On the other hand, the waxes correlations seem to be overestimating, with noticeable upward bending of the predicted upper APE, making the envelope less probably to close. Similar observations were reported in the work of Shoukry et al. (2019).
Several studies have been made regarding the status of carbon atoms in asphaltene molecules (Andrews et al. 2011;Molina et al. 2017), where the measurements include either direct carbon aliphaticity and aromaticity (as percentage) or a factor (referred to as aromaticity factor). These measurements quantify, accurately, the amount of aromatic and   saturated carbons in an asphaltene molecule. The measurements point out that asphaltene molecules in crude oils contain both aromatic and saturate carbons. Therefore, it is only logical that averaging between waxes and aromatics properties will give the best match. Which average to use depends on the oil sample. According to the samples in this study, Wax-Asp average appears to be the best average (least mean absolute percentage error or MAPE). This is probably because it shows true averaging of using waxes correlations, and modeling asphaltenes. In case the wax content is not reported, S-ARA averaging seems to be the better average. Table 8 shows the MAPEs of all samples using the different proposed methods. For samples, with computed AOPs, the reference AOP is not considered in the calculations of the MAPE.
Using an already performed characterization into the model (like for Samples 2, 3, and 4) exposes the applicability of the proposed modifications. The observations become, hence, less dependent on the characterization methods.
It is worth mentioning that the cubic EOS solid model is efficient to track asphaltene phase behavior; however, it fails to model asphaltene gradients. For gradient analysis, it is important to incorporate the Yen-Mullins model which specifies the nano-colloidal structures of asphaltenes (Freed et al. 2010).

Conclusions
The following conclusions can be drawn from this study: • An enhanced PR-cubic solid model has been developed to model asphaltene precipitation. The developed model includes two sets of thermodynamic parameters correlations for waxes and aromatics and three methods for empirical averaging of these parameters among the two sets. The averaging is based on the SARA analysis and wax content of the sample. • Using the aromatics correlations provides more rational upper APE trends than waxes correlations. • The averaging methods provided the least errors for all samples in the study. • Wax-Asp average (Average 3) is recommended, followed by S-ARA average (Average 1), for calculating the thermodynamic properties (fusion temperature and enthalpy). • The concept of averaging between waxes and aromatics properties agrees with the molecular structure of asphaltene molecules. • The proposed approach in cubic-PR solid model provides better predictability of AOP and precipitation amount outside the experimental range.
More testing on a larger dataset of oil samples is needed to consolidate the concept of averaging thermodynamic parameters within asphaltene precipitation modeling. The concept can, as well, be used/tested in other types of models [e.g., CPA models (Sattari et al. 2016)].

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.