Modelling the Diffusion Coefficients of Dilute Gaseous Solutes in Hydrocarbon Liquids

In this work, we present a model, based on rough hard-sphere theory, for the tracer diffusion coefficients of gaseous solutes in non-polar liquids. This work extends an earlier model developed specifically for carbon dioxide in hydrocarbon liquids and establishes a general correlation for gaseous solutes in non-polar liquids. The solutes considered were light hydrocarbons, carbon dioxide, nitrogen and argon, while the solvents were all hydrocarbon liquids. Application of the model requires knowledge of the temperature-dependent molar core volumes of the solute and solvent, which can be determined from pure-component viscosity data, and a temperature-independent roughness factor which can be determined from a single diffusion coefficient measurement in the system of interest. The new model was found to correlate the experimental data with an average absolute relative deviation of 2.7 %. The model also successfully represents computer-simulation data for tracer diffusion coefficients of hard-sphere mixtures and reduces to the expected form for self-diffusion when the solute and solvent become identical.


Introduction
The diffusion coefficient is one of the key transport properties controlling many natural and industrial processes such as transport though membranes, interfacial mass transfer and pore-scale processes in hydrocarbon reservoirs. In recent years, the diffusion coefficients of light gases in hydrocarbon liquids have received significant attention, especially in connection with enhanced oil recovery (EOR) and geological carbon dioxide storage (GCS) [1,2]. In EOR, natural gas, nitrogen or carbon dioxide is injected into an oil reservoir to maintain pressure, reduce oil viscosity and sweep the oil towards the production wells [3,4]. In GCS, carbon dioxide is injected into, for example, a depleted hydrocarbon reservoir where it may dissolve into the reservoir fluids (brine and residual hydrocarbons). The rate of gas dissolution is controlled in part by diffusion and the diffusion coefficients are relevant parameters for modelling the process [5]. The diffusion coefficients for binary gas-liquid systems have been studied by means of experimental and computational methods. Due to the scarcity of experimental measurements under high-pressure conditions, correlations based on hydrodynamic and kinetic theories have been adopted to estimate the diffusion coefficients [6]. Several attempts have been made to develop a universal correlation based on kinetic theory but, to date, the available methods are restricted in scope. Therefore, the objective of the present work was to extend a kinetic theory-based model to cover a wide range of solutes (inorganic and light hydrocarbon gases) and solvents (the broad class of hydrocarbon liquids) at infinite dilution.

Existing Modelling Approaches
Numerous theoretical and empirical models have been developed for the estimation of binary diffusion coefficients for both the general case and for gas-liquid systems in particular. The two most common approaches incorporating elements of theory are the hydrodynamic approach and the methods based on the theory of rough hardspheres. Both of these approaches are developed for the case in which one of the components is present at infinite dilution, corresponding to the so-called tracer diffusion coefficient.

Hydrodynamic Theory
The hydrodynamic theory is based on the Stokes-Einstein formula for the diffusivity of spherical particles in a continuum fluid and expresses the diffusion coefficient in terms of the thermal energy and the viscosity of the medium as follows: Here, k B is Boltzmann's constant, n SE is the Stokes-Einstein number determined by the boundary condition at the surface of the spherical particles, η is the viscosity of the fluid and a is the radius of the spherical particles. For a no-slip boundary condition, n SE = 6 [7] while, for a 'no-stick' boundary condition, n SE = 4 [8]. When applied to non-spherical particles, a becomes the effective hydrodynamic radius.
Although based on hydrodynamic theory for macroscopic particles, the Stokes-Einstein equation is readily applied to molecular solutes in molecular solvents, usually with the assumption that n SE = 4. The model then permits the diffusion coefficients of dilute solutes to be correlated with the viscosity of the solvent with a treated as an adjustable parameter. To obtain a good correlation, either a becomes temperature-dependent or the viscosity is raised to an empirical power [9].
This model is easily applied to correlate the diffusion coefficients of gaseous solutes in many different solvents as a function of temperature at constant pressure. When pressure is also a variable, matters are more complex and the observed behaviour is not captured well by making the hydrodynamic radius a function of temperature alone even when the pressure-dependence of the viscosity is included correctly. Cadogan et al. [10] and Taib et al. [11] found that a good correlation could be established in many cases by expressing the hydrodynamic radius as a function of the solvent density. However, such an approach was found to breakdown for highly viscous solvents such as squalane [10].

Kinetic Theory
The kinetic theory of dilute gases was extended to dense pure fluids by Enskog under the assumption of hard spherical molecules [8]. Thermophysical properties of real liquids can be estimated using hard-sphere models because the repulsive part of the intermolecular potential dominates the properties at high densities [9]. The smooth-hard-sphere (SHS) theory [12] can be applied to correlate the dimensionless viscosity, thermal conductivity and self-diffusion coefficients of pure fluids obtained from molecular simulations as functions of reduced volume V*. Here, V* is defined as V/V 0 , where V is molar volume and V 0 is the molar core volume which, for hard- As an example, the dimensionless self-diffusion coefficient D * 11 is defined by where D 11 is the self-diffusion coefficient, n is number density, and [nD 11 ] 0 = (3∕8 2 )(RT ∕ M) 1∕2 is the kinetic-theory expression for the dilute-gas limit of nD 11 . This dimensionless self-diffusion coefficient is correlated as a function of V* only: where F 11 (V*) is a universal function [13,14]. The same theory can be applied to real spherical molecules when V 0 is treated as a weakly temperature-dependent parameter. It is found that both molecular simulation data for smooth hard-spheres and experimental data for real monatomic gases conform to this model. Non-spherical molecules have internal energy that may be exchanged in collisions and such systems are treated as rough hard-spheres (RHS) [15]. In the RHS theory, an additional temperature-and density-independent 'roughness factor' appears such that, for example, the dimensionless self-diffusion coefficient is defined by where A 11 is the roughness factor and 0 ≤ A 11 ≤ 1. The universal function F 11 remains the same and experimental data for a given substance are correlated by determining the optimal values of A 11 and of a small set of parameters that describe V 0 as a function of T. The tracer diffusion coefficients D 12 for RHS systems can be treated theoretically in a similar way [12]. The method has been further developed and applied to real systems by Cadogan et al. [10] in terms of a dimensionless tracer diffusion coefficient D * 12 for solute 1 in solvent 2 defined as follows: here, A 12 is the corresponding roughness factor, is the reduced molar mass, V 0,2 is the molar core volume of the solvent, , and V 0,1 is the molar core volume of the solute. For hard spherical molecules, D * 12 is rigorously a function of the reduced volume V* and of the two ratios V 0,2 ∕ V 0,1 and 2 and σ i is the diameter of molecule of type i. However, analysis of the tracer diffusion coefficient data for smooth hard-sphere mixtures, obtained by molecular dynamics simulations [16], showed that D * 12 can be represented by a function of just two variables: the reduced volume V * = V ∕ V 0,2 and the 'asymmetry' ratio χ defined by [10] Therefore, in general, we have where F 12 (V * , ) is a universal function to be determined. An attractive feature of Eqs. 5 and 7 is that they treat tracer diffusion in binary systems in a similar way to the RHS theory for the transport properties of pure substances. In particular, the approach reduces to the established RHS correlation for the self-diffusion coefficient when the solute and solvent become identical.
The molar core volumes V 0,i (T) in Eq. 5 are expected to be identical with those determined from the pure-component properties. As discussed in the literature, these are weak functions of temperature and may be obtained by analysis of the transport properties of each pure substance in the high-density region where the rough hard-sphere theory is valid. Proceeding on this basis, Cadogan et al. [10] analysed their experimental mutual diffusion coefficient data for dilute solutions of CO 2 in various hydrocarbon liquids. These systems were characterised by a nearly constant asymmetry ratio χ = 2.1 and, for that class of mixtures, the data were represented by means of a simple polynomial: � .
The model therefore required just one system-dependent parameter, A 12 , and four universal constants a i , in addition to the molar core volumes of the pure components. The purpose of the present work is to extend this methodology to encompass a wider range of non-polar solutions and to obtain the dependence of F 12 on both V* and χ.
The RHS theory for tracer diffusion in binary systems has also been studied by Akgerman and co-workers in a series of papers published between 1987 and 1997 [17][18][19][20]. The models derived in that work corresponds to an assumption that F 12 is a linear function of V* and that the dependence upon χ can be absorbed into a prefactor. The final form of their model is very simple: where V D is a close-packed molar volume (generally different to V 0 ) at which diffusion effectively ceases and β is a scaling parameter. In an attempt to make the model predictive, β and V D have been correlated in various forms. For example, Eaton et al. [20] correlated β and V D with the critical constants of the solute and solvent in terms of universal constants which were adjusted by comparison with experimental data for the diffusion coefficient of 1-octene in supercritical alkane solvents. When compared with more than 1500 data points from the literature for other binary systems, the average absolute relative deviation (∆ AARD ) was found to be approximately 15 %. However, it turns out that experimental data (e.g. for CO 2 in squalane) do not conform to the linear dependence upon V embodied in Eq. 9. Therefore, we prefer the more general representation offered by Eq. 7.

Basis of the New Correlation
To extend the RHS model, we considered experimental data for the tracer diffusion coefficients of a range of gaseous solutes in hydrocarbon solvents, together with molecular dynamics simulations for systems of smooth hard-spheres. The simulation data anchor the correlation because, for these SHS systems, A 12 = 1 by definition. Table 1 summarises the data retrieved from the literature for the solute-solvent systems that were used in this study, including the temperature and pressure ranges, the number of data points, the measurement techniques and the estimated uncertainty of the measurements. Most of the systems were studied using the Taylor dispersion method which is generally a reliable technique with relative uncertainties of between 1 and 3 %. Higher uncertainties were found for systems studied using capillary cells and Mach-Zehnder interferometers. The most well studied systems are solutions of CH 4 and CO 2 in hydrocarbon liquids, with data extending over substantial ranges of temperature and pressure. The other systems present data mainly at atmospheric pressure.

Experimental Data
In order to analyse the experimental data, we required correlations for the molar core volumes as functions of temperature. These were mostly obtained from the literature and, for ease of reference, the relevant equations are reproduced in Table 2. New correlations for V 0 (T) were developed for nitrogen and argon in a form that presents reasonable extrapolation behaviour at high temperatures. For these components, it was found that the available viscosity data could be represented by the extended hard-sphere mode of Ciotta et al. using a simple-three-parameter model for V 0 (T) as follows: The development of these V 0 correlations is detailed in the Supplementary Information. Table 3 summarises the available molecular dynamics (MD) data for self-and mutualdiffusion coefficients in hard-sphere systems. Only three studies report data for mutual diffusion coefficients of binary SHS systems and, among these, the results of Easteal and Woolf [36] are the most comprehensive and accurate. The MD data are reported in the form of the ratio where n is number density, [nD 12 ] E is given by Enskog's expression,

Molecular Simulation Data
is the value of nD 12 in the dilute-gas limit, 12 = 1 2 ( 1 + 2 ) , and g 12 (σ 12 ) is the value of the unlike radial distribution function at contact. Hence, in terms of C 12 the dimensionless tracer diffusion coefficient is given by In order to recover D * 12 we obtained g 12 (σ 12 ) from the relation where = n 3 2 ∕6 = √ 2 ∕(6V * ) is the solvent packing density [37,38].

Results and Discussion
For the chosen combinations of solutes and solvents, the asymmetry ratio χ covers a range from approximately 0.8 to 3.1. These ratios are only weak functions of temperature and this is illustrated in Fig. 1 for five solutes (CH 4 , C 3 H 8 , N 2 , CO 2 and Ar) in octane. For a given solute, the values of χ in different alkane solvents are very similar, although the corresponding values in the aromatic solvent methylbenzene are about 15 % smaller. For the hard-sphere systems, the data set of Easteal and Woolf was limited to combinations of σ 1 /σ 2 ≤ 1 and M 1 /M 2 ≤ 1, representative of 'light' solutes in 'heavy' solvents. Figure 2 shows the simulation data for SHS as a function of χ for three different values of reduced volume. From this we deduced that the dependence upon χ is approximately linear. Therefore, the trial function for F 12 (V * , ) was a simple extension of Eq. 8: For hard-spheres, A 12 = 1 by definition while, for real fluids, A 12 is a constant for each solute-solvent pair. The parameters a i (i = 0, 3) and b 0 together with the roughness factor for each real-fluid solute-solvent pair were then adjusted in a non-linear optimisation. The objective function to be minimised was the un-weighted sum of the squares of the relative deviations of the data from the model. Because of the limited range of validity of the correlations for V 0 (T), the data in n-alkane solvents up to dodecane were limited to a maximum temperature of 411 K.
This simple model proved to provide a good account of the data as illustrated in Figs. 3 and 4 with deviations mostly within ± 5 %. Figure 3 compares the data for different solutes with the model as functions of the reduced volume V* of the solvent, where the model curve is plotted for the mean value of χ applicable to each solute gas. Figure 4 shows the relative deviations of the data from the model as functions of T, V*, M 1 /M 2 and χ. The comparison shows mostly good agreement but some evidence of systematic deviations for larger values of V*. However, generally, the deviation plots do not suggest that additional parameters in the universal function would improve the overall representation. In Fig. 3, the SHS data are for selfdiffusion (χ = 1), while all of the SHS data considered are included in Fig. 4. Table 4 reports the values of the universal coefficients determined in the fit, while Table 5 lists the roughness factors obtained and the values of the average absolute  * or self-diffusion coefficients D 11 * as function of the reduced volume V* of the solvent for different solutes. Experimental data for different solvents: hexane; heptane; octane; nonane; decane; dodecane; hexadecane; squalane; methylbenzene. Molecular simulation data for the self-diffusion coefficients of smooth hard-spheres (SHS).
---RHS theory plotted for a representative value of χ for each solvent system; ------correlation for the dimensionless self-diffusion coefficient of SHS [13] Application of the model to a new compound requires at least one experimental value of the tracer diffusion coefficient and sufficient pure-component data to determine V 0 (T). The model is only correlative because we presently have no means of predicting the roughness factors. Figure 5 shows some of the roughness factors determined for different solutes in alkane solvents. For CO 2 as the solute, there is a clear  trend with carbon number but the same cannot be said for CH 4 or C 2 H 6 . It is possible that new experimental data might result in smoother trends but, for now, the evidence is that A 12 does not depend systematically upon the carbon number of the solvent for alkane-alkane systems. On the other hand, the smooth trend for CO 2 -alkane systems suggest that the diffusion coefficients for CO 2 in other alkanes could be predicted. Although the model has been developed for hydrocarbon solvents with inorganic and light hydrocarbon solutes, the method might be expected to apply to other nonpolar solute-solvent systems. As presented, the correlation is restricted to light solutes. Preliminary investigation suggests that a more complex form of the universal function would be required to encompass also 'heavy' solutes in 'light' solvent and this might be developed using available experimental data for large molecular solutes in light supercritical solutes and also the molecular simulation data that were excluded from the present analysis (those with σ 1 /σ 2 > 1 and/or M 1 /M 2 > 1).

Conclusions
The model presented in this paper allows tracer diffusion coefficients to be rationalised within the well-established hard-sphere theory used to correlate and predict the transport properties of pure substances. To model a particular system one requires the molar core volume of both the solute and the solvent, as determined from pure-component transport properties, and the value of the roughness factor A 12 specific to the solute-solvent pair. Once A 12 is known, the model may be used to predict tracer diffusion coefficient over wide ranges of temperature and density. The method has been tested for non-polar systems comprising a light solute in a heavier solvent at reduced densities 1/V* ≥ 0.5. For CO 2 in alkane solvents, the roughness factor shows a clear systematic dependence upon the carbon number of the solvent but generally no clear relationship between A 12 and the nature of the solute and solvent emerges. We note that the available experimental data against which to test the model is quite restricted and that new measurements over extended ranges of temperature and pressure would be helpful. Future work might also focus on extending the model to wider ranges of size and mass ratios.

Acknowledgements
The authors gratefully acknowledge the financial support of Saudi Aramco for a scholarship.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.

Fig. 5
Roughness factors A 12 for various solutes in alkane solvents of carbon number C n : CO 2 ; CH 4 ; C 2 H 8