Hydrogen Solubility in Heavy Undefined Petroleum Fractions Using Group Contributions Methods

— Hydrogen solubility in heavy unde ﬁ ned petroleum fractions is estimated by taking as starting point a method of characterization based on functional groups [Carreón-Calderón et al. (2012) Ind. Eng. Chem. Res. 51, 14188-14198 ]. Such method provides properties entering into equations of states and molecular pseudostructures formed by non-integer numbers of functional groups. Using Vapor-Liquid Equilibria (VLE) data from binary mixtures of known compounds, interaction parameters between hydrogen and the calculated functional groups were estimated. Besides, the incorporation of the hydrogen-carbon ratio of the unde ﬁ ned petroleum fractions into the method allows the corresponding hydrogen solubility to be properly estimated. This procedure was tested with seven unde ﬁ ned petroleum fractions from 27 to 6 API over wide ranges of pressure and temperature (323.15 to 623.15 K). The results seem to be in good agreement with experimental data (overall Relative Average Deviation, RAD < 15%).


INTRODUCTION
Hydrogen solubility predictions are essential part in several processes in the petroleum and chemical industries; for instance, hydrocracking and hydrogenation processes are used to transform heavy oil into more valuable products.Hydrogen solubility data are required for designing and operating such processes and it is also needed in the corresponding kinetic models (Saajanlehto et al., 2014).
Several authors have pointed out the importance of taking into account Vapor-Liquid Equilibria (VLE) in reaction modeling and its influence in conversion predictions (Chen et al., 2011;Hook, 1985;Pellegrini et al., 2008).
A detailed review about VLE of hydrogen and petroleum fraction systems was presented by Chávez et al. (2014).They presented experimental data and modeling methods, most of which were conducted for systems composed of hydrogen and defined hydrocarbons, with few systems formed by hydrogen and undefined petroleum fractions.This may be attributable to the lack of experimental data, making a hard task to develop accurate prediction methods.Besides, binary systems containing hydrogen all exhibit a type III phase behavior in the classification scheme of Van Konynenburg and Scott (Privat and Jaubert, 2013) which cannot be predicted with cubic Equations Of State (EOS) using classical mixing rules because of the quantum nature of hydrogen and the mixture size-asymmetry (Deiters, 2013).Despite these issues, most methods found in literature use cubic Equations Of State (EOS) because of their mathematical simplicity and their capability of modeling multicomponent mixtures over a wide range of temperatures and pressures.To overcome such limitations several authors have proposed to fit EOS parameters to VLE data (Qian et al., 2013); however, in general, their capability of prediction is reduced with such regressions.New approaches have been suggested in order to preserve the predictive features of cubic equations of state, such as that described by Jaubert et al. (2010).In such approach, binary interaction parameters are predicted by using a group contribution method, where classical cubic equations of state and excess free energy models À EOS/g E Á are combined, being applied successfully to binary mixtures of hydrogen with defined hydrocarbons (Qian et al., 2013).
A further challenge is the mathematical characterization of the undefined petroleum fraction; that is, the determination of its physical-chemical properties by mathematical means.In such characterization process, it is important to find a convenient number of pseudocomponents in which the undefined petroleum fraction will be divided.According to Lin et al. (1985), better results are obtained as the number of pseudocomponents increases when VLE of systems formed by hydrogen and undefined fractions are modeled.On the other hand, Chávez concluded that 30 pseudocomponents are enough to obtain good agreement between calculated and experimental data.The characterization procedure also involves estimation of the corresponding critical properties entering into cubic equations of state.These properties are frequently calculated from correlations (Firoozabadi, 1988;Lee and Kesler, 1975;Pedersen et al., 2004), but they do not always lead to the same results, because each one was developed for different sets of fluids at limited conditions.Thus, recently, group contribution methods have been recently extended to pseudocomponents, showing a good performance (Carreón-Calderón et al., 2012;Xu et al., 2015).Carreón-Calderón et al. (2012) presented a new characterization procedure based on group contribution methods.They used selected functional groups to assign a hypothetical chemical structure to an undefined petroleum fraction by a minimization process of its corresponding Gibbs free energy.This methodology was applied to simulate VLE of heavy (Carreón-Calderón et al., 2014) and light hydrocarbons (Uribe-Vargas et al., 2016).Good results were achieved in both cases without making use of correlations to calculate critical properties or adjusted binary interaction parameters.In particular, they modeled vapor-liquid equilibria of heavy petroleum fluids, where no significant differences were shown regardless of the number of pseudocomponents used in modeling.These results show the feasibility of applying this characterization procedure to VLE of hydrogen with heavy undefined petroleum fraction.In this work, we add hydrogen functional group to the set of functional groups original proposed and suggest modifications to the characterization procedure in order to determine the hydrogen solubility properly.The entire methodology was tested with experimental data found in literature: Venezuelan Heavy Coking Gas Oil (HGCO) (Ji et al., 2013), Karamay Atmospheric Residue (KRAR) (Ji et al., 2013), Liaohe Atmospheric Residue (LHAR) (Ji et al., 2013), Venezuelan Atmospheric Residue (VNAR) (Ji et al., 2013), Canadian Light Virgin Gas Oil (CLVGO) (Cai et al., 2001), Canadian Heavy Virgin Gas Oil (CHVGO) (Cai et al., 2001) and Athabasca Bitumen (AB) (Lal et al., 1999).

CHARACTERIZATION PROCEDURE
The approach presented by Carreón-Calderón et al. (2012) is used to determine the properties of the undefined petroleum fractions.In that methodology, a hypothetical structure is calculated with a set of selected functional groups, which depend on the group contribution method used.Carreón-Calderón et al. (2012, 2014) have suggested two group contribution methods; Joback-Reid (JR) (Joback and Reid, 1987) and Marrero-Gani (MG) (Marrero and Gani, 2001).JR method was suggested to simulate VLE of light petroleum fluids (Uribe-Vargas et al., 2016), while MG method showed better functionality in heavy petroleum fluids heavier than 14 API (Carreón-Calderón et al., 2014).In the latter, the C 7+ fraction, which has unknown composition, was split in 3, 6 and 9 pseudocomponents, founding no significant difference between calculations and experiments regardless of the fraction splitting.Accordingly, in this work, the undefined petroleum fractions are treated as a single pseudocomponent in order to have a simpler characterization, where their corresponding properties are calculated through MG method.
Table 1 shows the functional groups of the MG method and the equivalent groups of the GCVOL (Elbro et al., 1991;Ihmels and Gmehling, 2003) method for determining liquid densities, as originally proposed by Carreón-Calderón et al. (2012).A quick review of these functional groups was made and it was found that the last functional groups are not equivalent each other (aC 6 ¼ AC À C).Therefore, in this work the last functional group (aC) was changed to aC À C functional group from the same MG method.
An optimization process is used to estimate a hypothetical chemical structure of the undefined petroleum fraction, where its Gibbs free energy is the objective function to be minimized.Since the undefined fraction is treated like a pure component, minimize the liquid fugacity coefficient / L i À Á is equivalent to minimize the Gibbs energy.Thus, the objective function is: where T C , P C and x are the critical temperature, critical pressure and acentric factor, respectively, with FG being the number of functional groups types.The / L i expression depends on the cubic EOS selected; in this work, the Peng-Robinson cubic Equation Of State PR-EOS (Peng and Robinson, 1976) is used to determine the fugacity coefficient as follows: where Z L = PV L ⁄(RT) is the compressibility factor of the liquid phase, A and B are PR-EOS parameters.The properties of undefined petroleum fractions are usually given at standard conditions; hence, the pressure P and temperature T are set equal to 0.101325 MPa and 288.15 K, respectively; v L is the molar liquid volume and R is the gas constant.This minimization problem is constrained by the molecular weight and density of the corresponding undefined petroleum fraction: In the above expression, subscripts j and i indicate the functional group j of the undefined fraction i, MW is the molecular weight, t j , the number (coefficient) of a functional group, Dv j the volume increment, q L i the liquid density.The Dv j is expressed as a function of temperature by means of where A 1 , A 2 and A 3 are the parameters of the GCVOL approach (Ihmels and Gmehling, 2003).Finally, to solve the minimization problem, Equation (2) requires critical properties (T C , P C and x).These properties are calculated using Equations from (6) to (8), where, the temperature is in K, pressures in MPa and volume in m 3 /mol (Marrero and Gani, 2001).
In previous equations t Ci , p Ci and v Ci represent contributions to the critical temperature, the critical pressure, and the critical molar volume because of the presence of a defined functional group i in a molecule.
One of the major contributions of the characterization procedure described above is the assignment of a hypothetical chemical structure to an undefined petroleum fraction.This allows us to handle the undefined fraction as hypothetical pure component and, in principle, to make use of any thermodynamic model exclusively developed for known components and their mixtures.Thus, in addition to PR-EOS, Huron-Vidal (Michelsen, 1990) mixing rule together with the UNIversal quasi-chemical Functional group Activity Coefficients (UNIFAC) approach, as activity coefficient model, are used also in this work for calculating VLE of hydrogen and heavy petroleum fluids, where the group-interaction parameters W nm in the UNIFAC approach are given by The parameters a nm , b nm and v nm are specific parameters for a couple of functional groups n and m, whose values are obtained from literature (Horstmann et al., 2005).Notice that a nm 6 ¼ a mn , b nm 6 ¼ b mn , v nm 6 ¼ v mn .It is worth pointing out that hydrogen is handled as a single functional group.

SOLUBILITY OF HYDROGEN IN DEFINED HYDROCARBON COMPONENTS
Before proceeding with the calculation of VLE of systems composed of hydrogen and heavy petroleum fluids, the hydrogen solubility in defined hydrocarbons was estimated in order to evaluate the thermodynamic models and their corresponding binary parameters described above.
Table 2 shows the 13 binary systems involved in this study, which are formed by nine paraffinic and four aromatic compounds.The range of experimental temperatures and pressures are also included and the comparison of calculated with experimental values are given in Table 3.
As observed, hydrogen solubility estimation is smaller than the experimental values in all cases.The difference is higher as the carbon number of the component increases; the Relative Average Deviation (RAD) increases from 13% for pentane to 75% for hexatetracontane, whereas, in the case of aromatic compounds increases from 32% to 52% for benzene and pyrene, respectively.Because of these results, group-interaction parameters W nm were calculated again by fitting calculations to the available VLE experimental data in Table 2.This fitting process requires to identify the functional groups forming the defined components.In case of paraffinic compounds, they can be represented by -CH 3 and -CH 2 functional groups, and aromatic compounds by ACH and AC functional groups.Because -CH 3 and AC belongs to the main group 1 (-CH 2 ) and 3 (ACH) respectively, only interaction parameters between -CH 2 and H 2 , ACH and H 2 , and ACH and CH 2 are required.Table 4 shows the original a nm parameters between main functional groups.
In the fitting process, for simplicity, parameter a nm is adjusted while b nm and v nm parameters keep their original values.The fitted group-interaction parameters are obtained by a minimization process, where the objective function (F Obj ) is given by Here P Sat Exp  fitting process is a group of new parameters (a nm , a mn ) for each binary mixtures (H 2 -hydrocarbon).Figures 1a and 1b show the trend of a nm and a mn parameters for paraffinic and aromatic compounds as function of the carbon number.
The correlated parameters take the form of: where, C k is the carbon number of the defined compounds.(2003); they used a Statistical Associating Fluid Theory (SAFT) model to predict hydrogen solubility in heavy alkanes, where the corresponding cross-binary interaction parameters were obtained as a function of the carbon number as well.
Figures 2a and 2b show experimental (r) and calculated hydrogen solubility data for pyrene and hexatetracontane, respectively.The original group-interaction parameters underestimate hydrogen solubility (----) whereas the predicted values are in better agreement with experimental data (-).
Table 3 summarizes the RAD of all defined components with the original group-interaction parameters and with those given by Equations from ( 11) to ( 14).

SOLUBILITY OF HYDROGEN IN UNDEFINED PETROLEUM FRACTIONS
As first step, critical properties and molecular pseudostructures of the seven undefined petroleum fractions studied in a) b) Figure 1 Group-interaction parameters as function of carbon number a) paraffinic compounds, b) aromatic compounds.
TABLE 4 Original group-interaction parameters (a nm , a mn ) between functional groups used in this study (Horstmann et al., 2005).this work were determined using the characterization procedure described in Section 1. Table 6 shows density, molecular weight and hydrogen-carbon ratio of the seven undefined fractions used in this work.As will be discussed later, the last two properties play an important role in modeling such systems with hydrogen.In all cases, the experimental error reported in the literature for hydrogen solubility was ±5% (Cai et al., 2001;Ji et al., 2013), except for AB system where such experimental error was not reported (Lal et al., 1999).
Once the molecular pseudostructure is formed with the functional groups from Table 1, it is possible to make a distinction between paraffinic, naphthenic and aromatic contributions through to the no-integer coefficients (t j ) resulting from the minimization process.Notice that although nine functional groups are used in the characterization procedure (see Tab. 1), there are only two main functional groups according to the UNIFAC method: first six functional belong to the main -CH 2 functional group, whereas the rest of them belong to the main ACH functional group.Accordingly, the calculated coefficients t j can be added in order to obtain the carbon number according to these two main functional groups.Therefore, Equations from ( 11) to ( 14) can be used to estimate group-interaction parameters between the H 2 and the functional groups calculated by the characterization procedure.
As illustration, Figure 3 shows the experimental and predicted hydrogen solubility in CHVGO (at 603.15 K) which has a molecular weight of 350 g/mol.The black diamonds (r) represent the experimental values and the crossed bars represent the experimental reported error of ±5%.Calculations using the original group-interaction parameters (----) underestimates the hydrogen solubility, as it is also observed for the defined hydrocarbons from Table 2; in this case the RAD is 51%.When our proposed correlations, Equations from ( 11) to ( 14), are used to estimate group-interaction parameters, the predicted values (Á Á Á) remains below experimental data, but the RAD is reduced up to 25%.To improve the accuracy of hydrogen solubility prediction, a modification of the original minimization process of characterization is suggested.This modification consists of including another constraint, in addition to Equations ( 3) and (4).We believed that this constraint should somehow account for the chemical nature of the undefined fractional.Hence, the Hydrogen-Carbon ratio (H/C) is proposed as a third constraint.According to experimental observations made by Ji et al. (2013), this parameter may be related to the aromaticity of undefined fraction.
Figure 4 shows the coefficients (t j ) obtained from the original characterization procedure and the modified one with the hydrogen-carbon ratio as constraint.The latter gives smaller values for the aromatic functional groups and larger values for all other functional groups.According to Park et al. (1996) hydrogen solubility decreases as the number of aromatic rings increases.Hence, we believe that the inclusion of the H/C ratio may help to construct a  pseudostructure in better agreement with its aromatic nature.
Figure 3 shows the predicted hydrogen solubility (-) using the modified characterization procedure.As seen, predicted hydrogen solubilities are in better agreement with experimental data (RAD = 5%).
The modified characterization procedure is summarized in Figure 5, where the stages are described in a flowchart.Here (H⁄C) j is the hydrogen-carbon ratio of a functional group j and (H⁄C) Exp is the experimental hydrogen-carbon ratio of a i undefined petroleum fraction.As observed in Table 6, three undefined fractions have molecular weights below 360 g/mol, whereas the rest of them have values above 500 g/mol.Although good results were achieved in the former, hydrogen solubility was overestimated in the latter.When the calculated pseudostructures are reviewed, it is found that they are formed mainly by paraffinic and cyclic functional groups and by a small portion of aromatic ones.This seems to be the reason why the hydrogen solubility is overestimated.Hence, a sensitivity analysis was made to know the effect of the H/C ratio in the characterization process from Figure 5.
In Figure 6a, three different pseudostructures of AB depending on the H/C ratio are shown.When H/C parameter is set equal to the experimental value, the pseudostructure is mainly formed by paraffinic and naphthenic functional groups and hydrogen solubility is overestimated (RAD > 50%).In contrast, hydrogen solubility is underestimated (RAD > 20%) when H/C parameter is set equal to  30% below the experimental value.In latter case, aromatic functional groups predominate over paraffinic and naphthenic ones.Figure 6b shows the trend of the RAD as function of the used H/C ratio.It can be seen that there is minimum value when H/C parameter is set around 80% of the experimental value; its corresponding pseudostructure is shown in Figure 6a.At this point, it is important to emphasize that the characterization procedure do not intend to make a characterization at molecular level, but look for an average structure that allow to reproduce bulk properties (Carreón-Calderón et al., 2012, 2014;Uribe-Vargas et al., 2016).
As with the A/B undefined petroleum fraction, the same sensitivity analysis was made for the VNAR, KRAR, and LAHR fractions.Figure 7 illustrates that minimum values of the RAD are achieved when H/C parameter lies between 80 and 85% of the experimental value.Therefore, we suggested to set the H/C ratio in this interval in order to obtained better results for heavy undefined petroleum fraction with molecular weights above 500 g/mol.
The results of hydrogen solubility prediction in Athabasca Bitumen (AB) at different temperatures and pressures are illustrated in Figures 8a and 8b.Here, the experimental value of the H/C ratio is 1.51, but, according to that previously discussed, it was set equal to 1.208, obtaining an overall RAD of 6.8%.In Table 7 all results are summarized, the overall RAD being equal to 14.12%.

CONCLUSIONS
The characterization procedure presented by Carreón-Calderón et al. (2012, 2014) was extended to predict hydrogen solubility in heavy undefined petroleum fractions including H 2 as another functional group, and the hydrogen-carbon ratio as third constraint in the characterization process.Group-interaction parameters between hydrogen and the main functional groups of the UNIFAC approach were calculated through a set of correlations based on experimental vapor-liquid equilibria data of binary mixtures of hydrogen and defined hydrocarbons.These modifications were enough for prediction purposes with molecular weights below 500 g/mol.For undefined petroleum fraction with molecular weights above 500 g/mol, the hydrogen-carbon ratio was modified to 80% of its experimental value.The inclusion of the hydrogen-carbon ratio was revealed as important property in order to construct the hypothetical structures in better accordance with its molecular nature.
) B Parameter of the equation of state b 1 -b 4

i
and P Sat Cal i are the experimental and calculated saturation pressure on the bubble point respectively and N represent the number of data points.The result of this Figure 2 Solubility of hydrogen a) pyrene at 433.2 K, b) hexatetracontane at 402.4 K.

Figure 5 Flow
Figure 5Flow chart of the modified characterization process.

Figure 4
Figure 4Coefficients calculated by characterization process.
Figure 6 Sensitivity analysis of H/C parameter used in the characterization process for AB.a) Pseudostructure arrangements of AB depending on percent of the H/C ratio, b) RAD vs. percent of H/C ratio used in the minimization process.

Figure 7
Figure 7Sensitivity analysis of RAD as a function of H/C parameter in the characterization process for different cuts.

TABLE 1
Original functional groups.

TABLE 2
List of 13 binary mixtures used in this study.
a H 2 ÀCH 2 and a CH 2 ÀH 2 are the group-interaction parameters between -CH 2 and H 2 functional groups, a H 2 ÀACH and a ACHÀH 2 are group-interaction parameters between ACH and H 2 functional groups.Finally, coefficients a n , b n , c n , d n are shown in Table5.Similar results have been obtained by Florusse et al.

TABLE 5
Group-interaction parameters (a nm ) between functional groups used in this study.

TABLE 6
Relevant physical properties of the undefined petroleum fractions used in this work.
Solubility of hydrogen in CHVGO petroleum fraction at 603.15 K as a function of group-interaction parameters.

TABLE 7
RAD on prediction of H 2 solubility in undefined hydrocarbon fractions.