1 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.

2 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 hard-spheres. 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.

2.1 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:

$$D = \frac{{k_{{\text{B}}} T}}{{n_{{{\text{SE}}}} \pi \eta a}}.$$
(1)

Here, kB is Boltzmann’s constant, nSE 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, nSE = 6 [7] while, for a ‘no-stick’ boundary condition, nSE = 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 nSE = 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].

2.2 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/V0, where V is molar volume and V0 is the molar core volume which, for hard-sphere of diameter σ, is \(N_{{\text{A}}} \sigma ^{3} {\kern 1pt} /{\kern 1pt} \sqrt 2\), where NA is Avagadro’s constant. As an example, the dimensionless self-diffusion coefficient \(D_{{11}}^{ * }\) is defined by

$$D_{{11}}^{ * } = \frac{{[nD_{{11}} ]}}{{[nD_{{11}} ]_{0} }}\left( {\frac{V}{{V_{0} }}} \right)^{{2/3}} = \left( {2N_{{\text{A}}} } \right)^{{1/3}} \frac{{8\sqrt \pi }}{3}\left( {\frac{M}{{RT}}} \right)^{{1/2}} \left( {\frac{{D_{{11}} }}{{V^{{1/3}} }}} \right),$$
(2)

where D11 is the self-diffusion coefficient, n is number density, and \([nD_{{11}} ]_{0} = (3{\kern 1pt} /{\kern 1pt} 8\sigma ^{2} )(RT{\kern 1pt} /{\kern 1pt} \pi M)^{{1/2}}\) is the kinetic-theory expression for the dilute-gas limit of nD11. This dimensionless self-diffusion coefficient is correlated as a function of V* only:

$$D_{{11}}^{ * } = F_{{11}} (V^{ * } ),$$
(3)

where F11(V*) is a universal function [13, 14]. The same theory can be applied to real spherical molecules when V0 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

$$D_{{11}}^{ * } = \left( {2N_{{\text{A}}} } \right)^{{1/3}} \frac{{8\sqrt \pi }}{{3A_{{11}} }}\left( {\frac{M}{{RT}}} \right)^{{1/2}} \left( {\frac{{D_{{11}} }}{{V^{{1/3}} }}} \right),$$
(4)

where A11 is the roughness factor and 0 ≤ A11 ≤ 1. The universal function F11 remains the same and experimental data for a given substance are correlated by determining the optimal values of A11 and of a small set of parameters that describe V0 as a function of T.

The tracer diffusion coefficients D12 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:

$$D_{{12}}^{ * } = \left( {2N_{{\text{A}}} } \right)^{{1/3}} \frac{{8\sqrt \pi }}{{3A_{{12}} }}\left( {\frac{{2M_{{12}} }}{{RT}}} \right)^{{1/2}} \left( {\frac{{V_{{0,12}} }}{{V_{{0,2}} }}} \right)^{{2/3}} \left( {\frac{{D_{{12}}^{{}} }}{{V^{{1/3}} }}} \right).$$
(5)

here, A12 is the corresponding roughness factor, \(M_{{12}} = M_{1} M_{2} {\kern 1pt} /{\kern 1pt} (M_{1} + M_{2} )\) is the reduced molar mass, V0,2 is the molar core volume of the solvent, \(V_{{0,12}} = (V_{{0,1}}^{{1/3}} + V_{{0,2}}^{{1/3}} )^{3} /8\), and V0,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}} {\kern 1pt} /{\kern 1pt} V_{{0,1}}\) and \(M_{2} {\kern 1pt} /{\kern 1pt} M_{1}\) where \(V_{{0,i}} = N_{{\text{a}}} \sigma _{i}^{3} {\kern 1pt} /{\kern 1pt} \sqrt 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^{ * } {\kern 1pt} = V{\kern 1pt} /{\kern 1pt} V_{{0,2}}\) and the ‘asymmetry’ ratio χ defined by [10]

$$\chi = M_{1} V_{{0,2}} {\kern 1pt} /{\kern 1pt} (M_{2} V_{{0,1}} ).$$
(6)

Therefore, in general, we have

$$D_{{12}}^{ * } = {\text{F}}_{{12}} (V^{ * } ,\chi ),$$
(7)

where \({\text{F}}_{{12}} (V^{ * } ,\chi )\) 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 CO2 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:

$$F_{{12}} (V^{ * } ,\chi = 2.1) = \sum\limits_{{i = 0}}^{3} {a_{i} \left( {V*} \right)^{i} } .$$
(8)

The model therefore required just one system-dependent parameter, A12, and four universal constants ai, 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 F12 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 F12 is a linear function of V* and that the dependence upon χ can be absorbed into a pre-factor. The final form of their model is very simple:

$${{D_{{12}} } \mathord{\left/ {\vphantom {{D_{{12}} } {\sqrt T }}} \right. \kern-\nulldelimiterspace} {\sqrt T }} = \beta (V - V_{{\text{D}}} ),$$
(9)

where VD is a close-packed molar volume (generally different to V0) at which diffusion effectively ceases and β is a scaling parameter. In an attempt to make the model predictive, β and VD have been correlated in various forms. For example, Eaton et al. [20] correlated β and VD 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 CO2 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.

3 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, A12 = 1 by definition.

3.1 Experimental Data

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 CH4 and CO2 in hydrocarbon liquids, with data extending over substantial ranges of temperature and pressure. The other systems present data mainly at atmospheric pressure.

Table 1 Summary of experimental diffusion coefficients data from the literature, where T is temperature, p is pressure, N is number of data points and Ur is expanded relative uncertainty (k = 2)
Full size table

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 V0(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 V0(T) as follows:

$$V_{0} = a - b\exp ( - T/c),$$
(10)
Table 2 Correlations for the molar core volume for different compounds as functions of temperature, where θ = T/K and C = carbon number
Full size table

The development of these V0 correlations is detailed in the Supplementary Information.

3.2 Molecular Simulation Data

Table 3 summarises the available molecular dynamics (MD) data for self- and mutual-diffusion 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

$$C_{{12}} = \frac{{[nD_{{12}} ]}}{{[nD_{{12}} ]_{{\text{E}}} }},$$
(11)
Table 3 Molecular dynamics simulation data for the self- and mutual-diffusion coefficients of smooth hard-sphere systems, where σi is the diameter and mi the mass for molecules of type i, V is the molar volume, \(V_{{0,i}} = N_{{\text{a}}} \sigma _{i}^{3} {\kern 1pt} /{\kern 1pt} \sqrt 2\), N is the number of data points and Ur is the expanded relative uncertainty (k = 2)
Full size table

where n is number density, [nD12]E is given by Enskog’s expression,

$$[nD_{{12}} ]_{{\text{E}}} = \frac{{[nD_{{12}} ]_{0} }}{{g_{{12}} (\sigma _{{12}} )}},$$
(12)

\([nD_{{12}} ]_{0} = (3{\kern 1pt} /{\kern 1pt} 8\sigma _{{12}}^{2} )(RT{\kern 1pt} /{\kern 1pt} 2\pi M_{{12}} )^{{1/2}}\) is the value of nD12 in the dilute-gas limit, \(\sigma _{{12}} = \tfrac{1}{2}(\sigma _{1} + \sigma _{2} )\), and g12(σ12) is the value of the unlike radial distribution function at contact. Hence, in terms of C12 the dimensionless tracer diffusion coefficient is given by

$$D_{{12}}^{ * } = \frac{{C_{{12}} }}{{g_{{12}} (\sigma _{{12}} )}}\left( {\frac{V}{{V_{0} }}} \right)^{{2/3}} .$$
(13)

In order to recover \(D_{{12}}^{ * }\) we obtained g12(σ12) from the relation

$$g_{{12}} (\sigma _{{12}} ) = \frac{1}{{1 - \xi }} + \frac{{3\xi }}{{(1 + \sigma _{2} /\sigma _{1} )(1 - \xi )^{2} }} + \frac{{\xi ^{2} }}{{(2\sigma _{2} /\sigma _{1} )(1 - \xi )^{3} }},$$
(14)

where \(\xi = \pi n\sigma _{2}^{3} {\kern 1pt} /{\kern 1pt} 6 = \pi \sqrt 2 {\kern 1pt} /{\kern 1pt} (6V^{ * } )\) is the solvent packing density [37, 38].

4 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 (CH4, C3H8, N2, CO2 and Ar) in octane.

Fig. 1
figure 1

Asymmetry ratios χ for five solutes in octane as functions of temperature T

Full size image

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 M1/M2 ≤ 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 \({\text{F}}_{{12}} (V^{ * } ,\chi )\) was a simple extension of Eq. 8:

$$F_{{12}} (V^{ * } ,\chi ) = \sum\nolimits_{{i = 0}}^{3} {a_{i} (V^{ * } )^{i} } + b_{0} \chi .$$
(15)
Fig. 2
figure 2

Dimensionless tracer diffusion coefficients of SHS systems as a function of the asymmetry ratio χ for various reduced volumes: V* = 1.5; V* = 1.7; V* = 1.9 [16]

Full size image

For hard-spheres, A12 = 1 by definition while, for real fluids, A12 is a constant for each solute–solvent pair. The parameters ai (i = 0, 3) and b0 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 V0(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*, M1/M2 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 self-diffusion (χ = 1), while all of the SHS data considered are included in Fig. 4.

Fig. 3
figure 3

Dimensionless tracer diffusion coefficients D12* or self-diffusion coefficients D11* 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]

Full size image
Fig. 4
figure 4

Relative deviations (Dexp − DRHS)/Dexp of experimental or simulated tracer diffusion coefficients from the RHS model as function of temperature T, reduced volume V/V0,2, molar mass ratio M1/M2 and asymmetry ratio χ, where 1 denotes solute and 2 denotes solvent. Symbols indicate solutes: × CO2; CH4; C2H6; C3H8; N2; Ar; smooth hard-sphere

Full size image

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 relative deviation (ΔAARD) and the maximum absolute relative deviation (ΔMARD). Entries without values of ΔAARD and ΔMARD are for systems with only a single data point available from which A12 was determined.

Table 4 Values of the parameters ai and b0 in the universal curve for tracer diffusion coefficients, together with the average absolute relative deviation (ΔAARD) and the maximum absolute relative deviation (ΔMARD) for the fit
Full size table
Table 5 Roughness factors A12 for different systems with absolute average relative deviations ∆AAD and maximum absolute relative deviations ∆MAD for the diffusion coefficient
Full size table

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 V0(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 CO2 as the solute, there is a clear trend with carbon number but the same cannot be said for CH4 or C2H6. It is possible that new experimental data might result in smoother trends but, for now, the evidence is that A12 does not depend systematically upon the carbon number of the solvent for alkane-alkane systems. On the other hand, the smooth trend for CO2-alkane systems suggest that the diffusion coefficients for CO2 in other alkanes could be predicted.

Fig. 5
figure 5

Roughness factors A12 for various solutes in alkane solvents of carbon number Cn: CO2; CH4; C2H8

Full size image

Although the model has been developed for hydrocarbon solvents with inorganic and light hydrocarbon solutes, the method might be expected to apply to other non-polar 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 M1/M2 > 1).

5 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 A12 specific to the solute–solvent pair. Once A12 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 CO2 in alkane solvents, the roughness factor shows a clear systematic dependence upon the carbon number of the solvent but generally no clear relationship between A12 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.