Self-diffusivity, M–S and Fick diffusivity of CO2 in Na-clay: The influences of concentration and temperature

Storing CO2 in underground saline aquifers is an important way to reduce CO2 emission in atmosphere, where gas/fluid diffusion in clay plays a key role in CO2 leakage and migration. Various diffusivities, self-diffusivity, Maxwell–Stefan (M–S) diffusivity and Fick diffusivity, in clay interlayer are investigated by molecular dynamics (MD). Self-diffusivity varies with CO2 concentration, and reaches the maximum value at 2 molecules/unit-cell. High fluid concentration leads to clay swelling, thereby increasing self-diffusivity. However, the fractional free volume of clay explains the trend of CO2 self-diffusivity, which does not decrease with CO2 concentration monotonously but reaches the maximum when CO2 concentration reaches 2. Displacement distribution of CO2 molecules is analysed to explore the microscopic diffusion mechanism, which is characterised by logarithmic normal distribution. The mean value of such distribution further explains the self-diffusivity dependence on CO2 concentration. M–S and Fick diffusivities of CO2 are calculated by MD for the first time, both of which increase with increasing CO2 and H2O concentration and temperature. Based on self-diffusivity and M–S diffusivity, a quantity representing the coupling strength between CO2 molecules is presented; it increases firstly with CO2 concentration but begins to decrease when CO2 concentration is beyond 2.

Carbon dioxide (CO 2 ) storage in underground saline aquifers is a promising way to reduce CO 2 emission in the atmosphere. This strategy provides long-term and large-scale storage of CO 2 . Clay minerals 1 , such as illite, chlorite, kaolinite and montmorillonite (MMT) play an important role in this process. On the one hand, owing to their porous (layered) structure, the clay minerals have remarkable capacity of adsorbing CO 2 2-4 . They provide a large amount of space for storing CO 2 5 . On the other hand, clay is almost impermeable, and therefore, the clay-enriched cap rocks show excellent sealing ability to retain injected CO 2 6-8 . Gas leakage and environmental impact are the main problems found after risk assessment of CO 2 storage. These problems are closely related to fluid (gas) transportation in clay. Gas diffusion in clay is the most important method for gas transportation because of the extremely low permeability of clay.
Many experimental studies have focused on the structure and interaction of CO 2 and clay mineral 9-12 . Loring et al. examined that CO 2 can migrate the interlayer region of MMT based on the in situ X-ray diffraction, magic-angle spinning in nuclear magnetic resonance spectroscopy and attenuated total reflection infrared spectroscopy 10,13 . Giesting et al. investigated the impact of CO 2 absorption on Ca-exchanged MMT expansion under different CO 2 pressures 14 . In addition, quasi-elastic neutron scattering experiments on hydrated clays have shown that hydrated cation diffusion mobility is probably a complex dynamic process, and the diffusion coefficients of the exchangeable cations were estimated 15 . Kozaki et al. determined the apparent diffusion coefficients of Cs + as functions of the temperature 16 . Sánchez et al. discussed that the self-diffusivity of water depends on the temperature and ionic strength of different kinds of clays 17 . However, no experiment on diffusivity of CO 2 in clay was reported.
With the current high-performance computational resources, molecular simulations have become powerful tools for understanding the molecular-scale structural 18,19 , thermodynamic 20 , mechanical 21 and dynamic [22][23][24][25] properties of clay. Grand-canonical Monte Carlo (GCMC) method was used to examine the adsorption of CO 2 with H 2 O 23 , CH 4 26-28 and organic molecules 25,29 in clay. Yang et al. used molecular dynamics (MD) to study the structure and self-diffusion coefficient (SDC) of CO 2 in uncharged clay-like slit pores 24 . Cygan et al. developed CLAYFF force field 30 for clay and three-site flexible potential modes of CO 2 31 . Myshakin et al. showed that the intercalation of CO 2 in clay caused the significant changes in the d-spacing and described that the distribution and diffusion of CO 2 , H 2 O and Na + were affected by the number of layers of clay [32][33][34] . Botan et al. found that CO 2 affects the diffusion of mobile species in clays 23 .
All the above-mentioned MD simulations mainly investigated the SDCs. Two important diffusion coefficients that are closely related to gas transportation, namely, Maxwell-Stefan diffusion coefficient (MDC) and Fick diffusion coefficient (FDC), have not been investigated. In fact, FDC, instead of SDC, is usually experimentally measurable and is used as parameter input in macroscopic gas transportation formulation, which predicts transportation of underground injected CO 2 . In this work, a series of MD simulations are performed to explore various diffusion coefficients of CO 2 in clay. The effects of gas concentration, water concentration and temperature on SDC, MDC and FDC are discussed in detail.

Methods
Theory of diffusion. As stated above, diffusion coefficients mainly have three types, SDC, FDC and MDC, which are useful in characterising molecular transport and mobility within porous materials. Self-diffusion means random motions or mixing of particles in the thermodynamic equilibrium. In MD simulation, SDC, D iself , of component i is computed by Einstein equation 35 : where n i is the number of the molecules of component i, d is the dimension of the system and r l , i (t) is the position vector of molecule l at time t. FDC is the particle motions driven by chemical or concentration gradient, resulting in the net mass transport. Fick's law of diffusion defines the transport diffusivity, D iFick , as the proportionality factor between the flux N and the concentration gradient: Chemical potential gradient is considered as the fundamental driving forces for diffusion [36][37][38][39][40] . For single-component diffusion in porous material, the transport equation can be expressed as follows: where μ i is the chemical potential of gas i. Moreover, L is obtained from the MD simulations using 41 the following: where V is the clay volume, k B is the Boltzmann constant and T is the temperature. The content in <…> represents an average on the cross-displacement correlation function (CDCF). The notation < > represents an average on a number of independent runs of same ensembles starting from different initial conditioins. The Maxwell-Stefan (M-S) theory of diffusion was already available from Maxwell 42 and Stefan 43 , in which gas diffusion is driven by chemical gradient. Using the M-S theory, the following expression can be derived for diffusion of the single component in a nanoporous material: where ρ is the clay density expressed as the number of unit cells per cubic meter, Θ i,sat is the saturation loading of the component i, μ i is the chemical potential expressed in joules per molecule and N i is the molar flux of the component i expressed in molecules per square meter per second. The fractional occupancy θ i are defined as follows: where Θ i is the loading of the component i. Combining equation (3) and equation (5), the MDC, D iM-S , is calculated by the following: Under isothermal condition, the molecules are considered to move with average velocity subject to a driving force ∇μ j ,μ = μ 0 + k B T ln f, where f is the fugacity and ∇μ = k B T ∇f/f. Thus, equation (2) is equivalent to equation (3). The FDC and MDC are related by  4 44 , which comprises 2:1 or tetrahedral−octahedral−tetrahedral layers. The substitution of octahedral Al 3+ by Mg 2+ and tetrahedral Si 4+ by Al 3+ leads to net negative layer charge. The interlayer Na + cations are balanced by these negative changes. Three-dimensional periodic boundary condition is applied to a simulation box comprising 32 (8 × 4 × 1) united cells. Clay model with different number of CO 2 and H 2 O are allowed to freely swell to their equilibrium d-spacing. All simulated systems are listed in Table 1. Cn (n = 1, 2,…) denotes the system in which the number of H 2 O is fixed, and the number of CO 2 is represented by n. Hn (n = 1, 3,…) represent systems with different numbers of H 2 O. Tn (n = 1, 3,…) represents systems with varying temperature.
In this article, all simulations are performed with the Accelrys Materials Studio software.The CLAFF force field 30,31 (Cygan, Liang et al. 2004, Cygan, Romanov et al. 2012) is used [30][31][32] . The Ewald summation method is applied for the electrostatic interactions. The atom-based summation method is used in the van der Waals interactions with a cut-off of 12.5 Å. Gas diffusions are simulated as follows: Gases are firstly inserted into the above equilibrium configuration by GCMC method to get a new clay-gas cell. The cell is minimised and equilibrated by NPT and NVT. Then, the cell is run in a NVT ensemble (0.5 ns) followed by a long-time (5 ns) NVE run. The trajectories in the NVE run are saved for diffusion coefficient computation. Integration time step is set to 1 fs. Since diffusion coefficients are based on statistics of numerous molecule trajectories, all the diffusivities in this article are computed by averaging on more than 10 indenpent MD runs to get reliable results.Gas concentration, water component and temperature are used to investigate their impacts on gas diffusions.

Results and Discussion
MD simulation can be used to obtain not only gas diffusion coefficients but also the microscopic details of the diffusion process, which is helpful to understand underlying diffusion mechanism in a confined space. Figure 1 depicts a simulation snapshot of the diffusion of CO 2 and H 2 O in MMT clay, from which it is seen that fluid molecules are uniformly distributed in the interlayer of MMT. The direction perpendicular to the walls of interlay is defined as z-axis and the mutually orthogonal x-and y-axes form a plane parallel to the walls. The short length along the z-axis greatly hinders the diffusion of gases along this direction. Thus, only the diffusion behaviour of fluids along the x-y plane is investigated in this study. Unless stated otherwise, the diffusion coefficients in this work refer to the diffusion in the x-y plane.
Self-diffusion coefficients and its Variation with CO 2 concentration. Diffusion coefficients are calculated based on the statistical average of the motion of large amounts of particles, and their accuracy relies on the sample size. For the sake of accuracy, each system is run 10 times independently. The resulting MSDs are averaged to calculate SDC. Diffusion coefficients are related to many factors including type and concentration of the gas, size of the pores, temperature and so on. Injection pressure of CO 2 , moisture content and underground temperature are key factors in the geological storage of CO 2 . Therefore, this study focuses on the impacts of CO 2 concentration, water concentration and temperature on CO 2 diffusion coefficient. Firstly, the influence of CO 2 concentration on its SDC is studied. The number of water molecules is set to 5 molecules per unit cell and the temperature to 313.15 K (approximately the temperature of 1000 meters underground). The Monte Carlo method is used to insert varying amounts of CO 2 molecules into the system. The system is then relaxed to equilibrium through NPT and NVT run, with resulting interlay distances (d-spacing) shown in Table 1. The equilibrated system is then subjected to dynamic run in NVE ensemble to compute SDC, as shown in Fig. 2a. With the increase of CO 2 concentration, the SDCs of CO 2 and H 2 O firstly increase and then decrease. At the CO 2 loading value of 64 (C6 system), the SDCs of both reach maximum values. The SDC of H 2 O remains higher than that of CO 2 , which originates from the smaller thermodynamic diameter of H 2 O (2.65 Å) compared with that of CO 2 (3.3 Å). However, many studies had shown that the SDCs of gases in porous media decrease with the increase of their concentrations 36,38,45 , which is different from the aforementioned result in Fig. 2a. To explore the cause of the variation of SDC, we further analyse the effect of CO 2 on the clay structure. As shown in Table 1, the d-spacing of MMT continues to increase with CO 2 concentration, which could provide fluids with larger space for diffusion, resulting in larger diffusion coefficient. SDC is expected to increase with CO 2 concentration  monotonously. However, our simulations reveal a higher SDC in high-concentration system C6 compared to the low-concentration system C5. To further explore the cause of this unusual SDC variation, fractional free volume (FFV) of the system is defined as the ratio of interlayer space unoccupied by fluids to the total volume of MMT: free total where V free is the free volume and V total is the total volume of the MMT. The FFV can be directly determined by molecular probes 46,47 . As shown in Fig. 3a, the FFV increases with the CO 2 concentration firstly and peaks at 64 CO 2 molecules, which is in accordance with the position of maximum of CO 2 SDC. This finding indicates two counter-interactions, as follows. Firstly, the increase of CO 2 concentration could cause the expansion of d-spacing in the system and provide fluids with more space for diffusion. Secondly, increasing CO 2 molecules occupy more space, thereby resulting in the decrease of free volumes. The fluid molecules become more crowded and diffuse more slowly. Therefore, fluid SDC depends on the free space in the interlayer. At the initial stage of CO 2 injection, the rapid expansion of d-spacing increases the free volume in the interlayer, which is beneficial for fluid diffusion and increases CO 2 SDCs. The free volume reaches a maximum at 64 CO 2 molecules. The free volume decreases despite the continued expansion of d-spacings, which hinders gas diffusion and decreases the SDCs of CO 2 and H 2 O (Fig. 3a).
Variation of SDCs with moisture and temperature. Water content varies with the depth of storage aquifer. A MMT system with 2 CO 2 molecules per unit cell is used to investigate quantitatively the impact of water concentration on the diffusion. Moreover, the temperature is fixed at 313 K. Figure 2b shows Fig. 3b. The FFV increases with moisture content, showing the same trend as the SDCs. The increase of water molecules increases the free volume in the interlayer and supports the diffusion of fluids. Temperature is an important factor in the diffusion of gases. As shown in Fig. 2c, increasing temperature can accelerate the thermal motion of molecules and increase their kinetic energy, thereby gradually increasing fluid SDC. Calculations show that the relationship between SDC and temperature follows the classical Arrhenius formula: where D 0 is a constant, and E a is the activation energy of diffusion. According to the equation, the activation energy for CO 2 diffusion in MMT is 20.26 kJ/mol (see inset of Fig. 2c), which is greater than its activation energy of diffusion in pure water 17.8904 kJ/mol 48 . The MMT layers impose restriction on CO 2 diffusion, causing high activation energy in MMT.
Details such as the paths and displacements of molecules provided by MD simulation can be used to study the microscopic mechanism of gas diffusion. Previous work indicated that the diffusion of gases in porous media (such as polymeric compounds and metal organic frameworks) involves two processes: limited movements within one micropore (small displacements) and jumps among micropores (large displacements) 37 . Layered MMT provides anisotropic confined space without interconnected micropores. Diffusion of gas in such structure is different from random movements of gas in unlimited environment. We further analyse the displacement distribution of CO 2 molecules, as shown in Fig. 4. All the distributions are seen to have a single peak, and no peak is shown at large displacements. The displacements generally fall within 0.6-1.0 Å. This finding indicates that no inter-pore jumps occur during gas diffusion. Figure 4 also demonstrates that, as CO 2 concentration increases, the displacement peak (i.e. the prevailing displacement) moves towards the left and becomes wider, reaching the maximum width in the C5 system. This finding implies increases of both the prevailing displacement and the number of fast-speed molecules. Therefore, with the increase of CO 2 concentration, SDC of CO 2 increases and reaches a maximum at C5 system. Similarly, the increase of moisture and temperature increases prevailing , an important characteristic parameter for log-normal distributions, is calculated. The results are shown in Fig. 5b-d. Figure 5b shows the relationship between mean displacement and CO 2 concentration. In this study, R av increases with CO 2 concentration firstly, reaches a maximum at C5 and then decreases,  showing the same trend as the SDC of CO 2 . The average displacement of molecules in unit time (1 ps) reflects the relative speed of diffusion, and hence, larger R av results in faster diffusion speed and a higher SDC. In a word, R av is positively related to the SDCs of CO 2 and H 2 O. This finding is further validated by systems with varying moisture contents and temperatures. As shown in Fig. 5c,d, R av increases with moisture and increases linearly with temperature, causing the increase of SDCs of CO 2 and H 2 O (Fig. 2).

Variation of MDC and FDC with concentration, moisture content and temperature.
To calculate MDC, the Onsager coefficient L should be determined using equation (4) in advance. The CDCF can thus be computed from the slope of the curve, which is used to calculate MDC from equation (7) and FDC using equation (9). The thermodynamic factor in equation (9) requires the adsorption isotherm of CO 2 , which can usually be simulated using the GCMC method. The conventional GCMC method assumes a rigid absorbent (μVT ensemble). However, MMT clay swells after adsorbing CO 2 . Swollen MMT can continue to adsorb more gas. Therefore, the conventional GCMC is not suitable for calculating adsorption isotherms of CO 2 in MMT. In this study, we adopt a new elastic-adsorption approach to simulate CO 2 isotherm by considering the swelling of MMT. The MMT structure with the initial 224 water molecules is taken as an example (see Fig. 6). Using the conventional GCMC method, we can obtain significantly different adsorption isotherms at different fixed d-spacing values, i.e. the initial and swollen d-spacing. In these isotherms, the curves rapidly approach to saturation at low pressure, and the amount of CO 2 adsorption is almost kept constant when the pressure rises above 1 MPa (shown by the black line in Fig. 6), which is inconsistent with experimental results 49 . However, the CO 2 adsorption isotherm obtained by the above elastic-adsorption method increases gradually with the CO 2 pressure, which is more qualitatively coincident with the Langmuir formula and the above-mentioned experiment results. The maximum adsorption θ sat can be obtained by fitting with the Langmuir equation, and then, the thermodynamic factor can be calculated by equation (10). The results show that thermodynamic factor always increases with increasing CO 2 concentration. For low CO 2 concentrations, thermodynamic factor is not influenced by H 2 O apparently, whereas at higher CO 2 concentrations, the thermodynamic factor gradually decreases with H 2 O concentration.
The effect of CO 2 concentration, H 2 O content and temperature on MDC and FDC is also explored. The red dots in Fig. 7 represent the variation of CO 2 MDC, which increase with CO 2 concentration, H 2 O content and temperature. The black dots in Fig. 7 represent the FDC. The thermodynamic factor is generally high. Thus, the FDC is usually larger than MDC. Similar to MDC, FDC increases with CO 2 concentration, H 2 O content and temperature.
The significance differences of CO 2 diffusivities among CO 2 concentration moisture content and temperature are calculated as 0.033 (SDC), 0.027 (MDC) and 0.012 (FDC) respectively, indicating the reliablity of the above results.
Relationship between SDC and MDC: The coupling effect of CO 2 molecules. Krishna and Paschek et al. proved that the difference between MDC and SDC is caused by the self-exchanging coupling between molecules of a diffusing gas 40 . The strength of the coupling between gas molecules can be represented by the self-exchange coefficient D ii , which is calculated by the following formula: Using the aforementioned SDC and MDC results, the impact of CO 2 concentration on D ii is shown in Fig. 8. With increasing CO 2 , D ii firstly increased and then decreased, reaching the maximum at 2 molecules/unit cell. D ii increases with increasing CO 2 molecules, thereby leading to weak correlation among CO 2 molecules. When Figure 6. Adsorption isotherms of CO 2 by conventional GCMC method (black) and flexible adsorption method (red) in C1 system. the CO 2 concentration exceeds 2, D ii begins to decrease, resulting in stronger correlation. This result is different from the case described in ref. 40, where D ii monotonically increased with gas concentration in MOF. The coupling between molecules may have been influenced by two factors. On one hand, the correlation is strengthened by increasing the CO 2 concentration. On the other hand, the rising concentration causes the MMT to expand, thereby increasing the volume of the layer and weakening the correlation. The final strength of the coupling may depend on the FFV within the MMT.

Conclusions
In this study, MD is used to study various diffusion coefficients of CO 2 in MMT clay under varying conditions. With increasing CO 2 concentration, the SDCs of CO 2 and H 2 O display a peak, respectively. Water content and temperature are positively correlated with the SDC. Many previous studies show that SDCs of gases in porous media usually decrease monotonically with the gas concentration. To further explain the unusual results, the FFV within MMT and the displacement distribution of CO 2 are analysed. FFV has an important effect on the diffusion of gas molecules in MMT. Increasing CO 2 concentration causes the expansion of MMT, thereby increasing the  internal FFV of MMT and providing gas molecules with more space for diffusion. This finding explains the initial increase of the fluid SDC in MMT. FFV begins to decrease when the CO 2 concentration rises above 2, which hinders the gas diffusion and thereby decreasing SDC. Similarly, the increase of water content also causes the MMT to expand and create more free space for CO 2 diffusion, leading to increasing CO 2 SDC. The fluid SDC in MMT depends on the FFV of MMT. Furthermore, the displacement distribution of CO 2 follows log-normal distribution, and the mean of the distribution shows the same trend as the SDC, which provides a good explanation of the effects of CO 2 concentration, moisture and temperature on SDC from another perspective.
Most previous works focused on the self-diffusion coefficients of gases in MMT. FDCs are of practical importance for the geological storage of CO 2 . In this work, the MDC and FDC of CO 2 in MMT are investigated by MD for the first time, both of which increase with CO 2 concentration, moisture and temperature. Using the SDC and MDC, we calculate the self-exchange coefficient D ii , which represents the coupling effect among CO 2 gas molecules. This coefficient firstly increases with CO 2 concentration, reaching a maximum with concentration 2, and then decreases, showing the same trend as FFV. This implies that, with increasing CO 2 concentration, the FFV increases, thereby weakening the correlation between CO 2 molecules until the CO 2 concentration reaches 2. Then, the FFV decreases and the CO 2 correlation coupling strength begins to increase.