Mass transport modelling of two partially miscible, multicomponent fluids in nanoporous media

Abstract

High-pressure fluid transport in nanoporous media such as shale formations requires further understanding because conventional continuum approaches become inadequate due to their ultralow permeability and complexity of transport mechanisms. We propose a species-based approach for modelling two partially miscible, multicomponent fluids in nanoporous media – one that does not rely on conventional bulk fluid transport frameworks but on species movement. We develop a numerical model for species transport of partially miscible, non-ideal fluid mixtures using the chemical potential gradient as the driving force. The model considers the binary friction concept to include the friction between fluid molecules as well as between fluid molecules and pore walls, and incorporates the key multicomponent transport mechanisms – Knudsen, viscous and molecular diffusion. Under single-phase conditions, the system under consideration is quantified by introducing multicomponent Sherwood number (Sh), Péclet number (Pe) and fluid–solid friction modulus (φ). Despite the complexity of fluid transport in nanopores, the steady-state single-phase transport results reveal the contribution of diffusion in nanopores, where all parameters collapse on a set of master curves for the multicomponent Sh with a dependence on multicomponent Pe and φ. Unsteady state, two-phase transport modelling of the codiffusion process shows that light and intermediate alkanes are produced much higher than heavy alkanes when the vapour phase appears. We demonstrate that the pressure gradient is also crucial in promoting CO₂ and alkane mixing during counterdiffusion processes. These results stress the need for a paradigm shift from classical bulk flow modelling to species-based transport modelling in nanoporous media.

1. Introduction

A thorough understanding of the mass transfer process in nanoporous media is crucial for better understanding and optimizing many natural and industrial processes, such as optimizing membrane and solid oxide fuel cell systems (Young 2007; Yuan & Sundén 2014), designing separation and catalytic processes (Krishna 2012), enhancing hydrocarbon production from shale reservoirs (Falk et al. 2015), geo-storage of carbon dioxide (CO₂) in shale formations (Lee, Cho & Lee 2020; Emami-Meybodi et al. 2025) and long-term integrity of CO₂-reservoir caprocks (Kampman et al. 2016; Adu-Gyamfi et al. 2022). Various conventional continuum approaches have been adopted to simulate multiphase, multicomponent fluid transport within nanoporous media at elevated pressures. However, due to their ultralow permeability, reliable mass transport models are missing to predict fluid transport behaviour in these porous media.

As summarized in table 1, depending on the dominant transport mechanism, a mass conservation equation is usually solved combined with a constitutive equation, such as Darcy's law, Fick's law, the Maxwell–Stefan (MS) approach or the dusty gas model (DGM). Generally, the total mass flux can be expressed as the sum of advective and diffusive contributions. The contribution of diffusive flux in the overall system behaviour is usually negligible for high-permeability porous media with high Péclet numbers (Fan et al. 2014). Hence, advection caused by the pressure gradient dominates mass transport in high-permeability porous media. According to Darcy's law, the advective flux (see table 1) can be related to gradient pressure and viscosity (Darcy, 1856). Various Darcy-based transport models with apparent permeability (Darabi et al. 2012; Wang, Wang & Chen 2018) have been developed to simulate multicomponent fluid transport in nanoporous media, particularly with the applications to hydrocarbon production from shale reservoirs and CO₂ injection into these formations (Chen, Balhoff & Mohanty 2014; Ozowe, Zheng & Sharma 2020). However, Darcy's law is limited to describing the bulk fluid motion and fails to predict transport in shale’s nanoporous matrix (Falk et al. 2015).

Table 1. The constitutive equations used to express molar flux, N_i, in fluid transport models. The subscript i represent the parameters of component i; N is molar flux; c is molar concentration; ${c_t}$ is total molar concentration; K is permeability; $\eta $ is viscosity; p is pressure; $\rho $ is fluid mass density; $\textrm{g}$ is gravitational constant; D is diffusion coefficient; x is mole fraction; R is gas constant; T is temperature; $\mu $ is chemical potential; ${n_c}$ is the number of components; $ \unicode{x0110}_{ij}$ is MS diffusion coefficient of pair component i and j; $ \unicode{x0110}_{iM}$ is Knudsen diffusion coefficients.

The advective transport mechanism is less effective than the diffusive transport mechanisms in nanoporous media (Cronin, Emami-Meybodi & Johns 2018; Liu & Emami-Meybodi 2021). Consequently, diffusive transport has become a crucial transport mechanism for species movement in these formations, which is widely represented by Fick's law. Fick's law was proposed by Fick (1855), inspired by the experiments of Graham (1850) on liquid and gas diffusion. The subsequent derivation and assumptions of Fick's law are elucidated through several approaches, including the continuum mechanics approach (Truesdell 1962; Valdes-Parada, Wood & Whitaker 2023), kinetic theory (Kerkhof & Geboers 2005a) and the thermodynamic approach (Curtiss & Bird 1999; Firoozabadi 2016). In nanoporous media, diffusion-based compositional models using Fickian-type modelling have been proposed (Cronin, Emami-Meybodi & Johns 2021) with a more relaxed set of assumptions. In these models, the total flux is expressed by diffusive flux (see table 1), which can be related to the effective/apparent Fickian diffusion coefficient (Liu & Emami-Meybodi 2022) and the concentration gradient. Moreover, several studies have investigated the impacts of both advective and diffusive fluxes by simulating production and solvent injection processes through Darcy's law and classical Fick's law (Moortgat et al. 2013; Zidane & Firoozabadi 2021). These investigations have shown that molecular diffusion plays a significant role in hydrocarbon recovery of solvent injection in nanoporous media. However, Fick's law, which was initially developed for binary component diffusion, may not fully capture the multicomponent mass transport in a liquid or vapour mixture, given the diffusional interactions between each pair of components (Krishna & Wesselingh 1997; Deen 2013). Therefore, using Fick's law may not suffice to describe multicomponent mass transport in nanoporous media (Krishna & Wesselingh 1997; Corral-Casas et al. 2022).

Based on the hypothesis that two equally counterbalanced driving forces exist, the MS approach is often used to express mass transport in multicomponent systems. The inception of this approach traces back to Maxwell (1860), who proposed an alternative expression of Fick's law for binary diffusion. Subsequently, Stefan extended this formulation into a generalized Maxwell equation, applicable to N-component diffusion, commonly referred to as the MS equation (Bird & Klingenberg 2013). The derivation of the generalized MS equation, encapsulating concentration, thermal, pressure and forced diffusion, is proposed in the principles of irreversible thermodynamic theory (Curtiss & Bird 1999) and continuum mechanics (Whitaker 2009). This study focused on concentration and pressure diffusion where the total flux related to the effective MS diffusion coefficient and the chemical potential gradient is directly calculated through the MS equation (see table 1). It can also be explicitly expressed in the form of Fick's law, known as the generalized Fick's law (Taylor & Krishna 1993; Moortgat & Firoozabadi 2013). The MS diffusion formulation with a matrix of composition-dependent diffusion coefficients has been used in simulating multicomponent fluid transport in porous media (Tian et al. 2021). However, the wall effect is not considered in MS formulation. Under the nanopore confinement conditions, the interaction between molecules and pore walls cannot be neglected. Therefore, modifications to the MS formulation are needed to incorporate the impact of molecular interactions with the nanopore walls.

The conventional no-slip boundary condition requires modification due to the distinctive transport mechanisms governing multicomponent fluid transport in nanopores. In nanopores, phenomena such as Knudsen diffusion, diffusion slip and thermal slip induce a velocity parallel to the pore wall. When the pore size decreases to the nanoscale, the mean free path may become comparable to or even greater than the pore diameter, so that collision between a molecule and the pore wall occurs more frequently with intermolecular collision (Do 1998; Gruener & Huber 2008; Kärger et al. 2012). The stochastic nature of this process gives rise to a Fickian expression for the flux in which the diffusivity depends only on the pore size and the mean molecular thermal velocity (Knudsen 1909). This is known as Knudsen diffusion. Additionally, in the context of multicomponent mixture transport, a concentration gradient parallel to the wall exists, resulting in a slip velocity along the wall known as diffusion slip. The diffusion slip velocity was first derived by Kramers & Kistemaker (1943) in connection with equimolar counterdiffusion. Subsequently, it has been applied to investigate multicomponent gas flow in capillaries and porous media (Young & Todd 2005; Mills & Chang 2013). Similarly, thermal diffusion is induced by a temperature gradient parallel to the wall. All the aforementioned mechanisms contribute to slippage of mass flux, a phenomenon that cannot be neglected for multicomponent transport in nanopores.

The DGM is a convenient approach for considering the effect of pore walls on fluid transport in porous media (Krishna & Wesselingh 1997; Mason 1983). In this approach, the rock grain is assumed to be a giant dummy molecule ‘dust’ in a mixture. The total flux (see table 1) can be related to the effective MS diffusivities, the effective Knudsen diffusivities and the chemical potential gradient (Krishna 1987; Krishna & Wesselingh 1997; Mason 1983). The DGM has succeeded in fields such as catalysis, membrane (Gao et al. 2011) and fuel cell modelling (Yuan & Sundén 2014), but typically only applies to low-pressure conditions. Moreover, DGM has faced questions concerning its inherent inconsistencies because of the different frames of reference chosen in different steps during its derivation (Kerkhof 1996; Bhatia, Bonilla & Nicholson 2011; Pant, Mitra & Secanell 2013). Recent studies have implemented DGM to study fluid transport in nanoporous media (Achour & Okuno 2022; Sena Santiago & Kantzas 2022) providing valuable insights. However, these studies have often considered idealized systems and boundary conditions, which may not accurately reflect real systems.

This study suggests a paradigm shift for mass transport modelling of a two-phase, partially miscible, multicomponent fluid in nanoporous media that considers key transport mechanisms for fluid transport in nanoporous media. We present a new perspective that relies on the transport of species based on the binary friction concept to reconcile Knudsen, viscous and molecular diffusion. A wider range of flow regimes, from continuum to Knudsen transport, can be accounted for under changing flow conditions. The gradient in chemical potential drives the component movement, irrespective of the number and miscibility of phases. For the first time, we use this species movement approach, instead of bulk fluid transport models, to simulate fluid transport in the nanoporous media under two-phase (vapour and liquid) fluid conditions where component mass transfer exists between phases, i.e. a partially miscible system. The proposed approach has the capability to simulate multiphase, multicomponent transport within nanopores under high-pressure conditions, taking into account non-ideal mixing and nanopore wall effects.

2. Methods

2.1. Species mass transport

Figure 1 compares multicomponent fluid transport in a pore based on species movement (figure 1b) instead of the conventional approach of adding the species transport to bulk fluid flow (figure 1a).

Figure 1. Multicomponent fluid transport in nanopores. (a) Multicomponent fluid transport based on the combination of advective transport and diffusive transport. Advective flux, determined by Darcy's law, arises from pressure gradients within porous media. Diffusive flux, influenced by composition gradients, utilizes the advective average velocity as a reference. The total flux is obtained by summing the advective and diffusive components. (b) Multicomponent transport based on species movement concept. The chemical potential gradient that inherently includes the pressure and composition gradients is regarded as the driving force. Three kinds of friction will balance it: friction between identical molecules; friction between different molecules; friction between molecules and walls.

A commonly used multicomponent transport modelling approach is based on the combination of advective transport and diffusive transport (see figure 1a). Accordingly, the species fluxes are separated into two parts: the advective transport, which automatically accompanies any bulk motion; and diffusive transport, which is related to molecular motions (Deen 2013). This study employs the molar form of the species transport and continuity equation. The species molar flux can be expressed as

(2.1)\begin{equation}{N_i} = {x_i}{c_t}\bar{v} + {J_i},\end{equation}

where N_i represents the molar flux of species i, c_t is the total molar concentration of the mixture, x_i is the mole fraction in the mixture, $\bar{v}$ is the molar average velocity and J_i is the molar diffusive flux. The first term on the right-hand side of the equation expresses the bulk fluid advective transport, while the second term accounts for the diffusive transport of species.

Darcy's law calculates the mean velocity of a mixture driven by pressure gradients, while Fick's law computes diffusive flux driven by composition gradients in porous media. Nevertheless, various phenomena, such as slippage, surface diffusion, Knudsen diffusion and nanopore confinement effects, affect species transports in nanopores. Determining whether these mechanisms should be attributed to advection or diffusion becomes challenging. Furthermore, some mechanisms, including slippage, surface diffusion and Knudsen diffusion, may be reconsidered in transport processes, given their similar association with molecular and wall interactions. In addition, conventional continuum descriptions, such as Darcy's, may fail to predict transport in nanoporous media (Falk et al. 2015; Kavokine, Netz & Bocquet 2021).

We develop a model based on species transport (see figure 1b) without relying on bulk fluid flow by considering the chemical potential gradient as the driving force and three kinds of friction: friction between different molecules; friction between identical molecules; and friction between molecules and walls. Friction between different molecules is associated with molecular diffusion, expressed by MS equations. Friction between identical molecules arises from varying velocities in different layers related to viscous flow. Friction between molecules and walls is linked to Knudsen diffusion and slippage. The Knudsen number remains below 0.02 within a nanoporous medium owing to the extremely small mean free path for liquid phase. As a result, Knudsen diffusion is not considered in the liquid phase. However, experiments and molecular dynamic simulation have observed the fast mass transport of liquid (water and hydrocarbon) in nanopores (Holt et al. 2006; Wang, Javadpour & Feng 2016). Therefore, a slip modulus is considered to account for the fast mass transport of liquid in the nanopores.

2.2. Conceptual model

First, we establish a conceptual model (figure 2a) to simulate and analyse the steady-state multicomponent alkane mixture transport in an individual nanopore. Then, we extend our investigation to a more complex fluid transport system in nanoporous media by considering a multicomponent mixture transport in a one-dimensional (1-D) nanoporous medium consisting of a bundle of capillaries connected to a high-permeability porous medium (figure 2b). The presence of a connected high-permeability porous medium plays a pivotal role in defining the inner boundary condition within various multicomponent transport systems. This interconnected medium is a key aspect of diverse multicomponent transport systems such as the anode and cathode in solid oxide fuel cell systems (Pant et al. 2013), bulk tanks utilized in separation and catalytic devices (Veldsink, Versteeg & van Swaaij 1994) and fractures encountered within shale reservoirs (Hoteit & Firoozabadi 2009; Tian et al. 2021).

Figure 2. Conceptual model of multicomponent fluid transport in nanopores. (a) The conceptual model for multicomponent fluid transport in an individual nanopore with a radius of R. A constant chemical potential drop was imposed across the nanopore by settling different pressure and mole fractions at the inlet and the outlet. The species velocity profile can be calculated within the nanopore by solving a microscopic species momentum balance equation. (b) The conceptual model for multicomponent mixture transport within a 1-D nanoporous medium consists of a bundle of capillaries connected to a high-permeability porous medium.

In a nanoporous medium, the influence of pressure-driven viscous flow decreases significantly as the pore size decreases. Conversely, molecular diffusion takes on a more important role in multicomponent systems. Simultaneously, the fluid–wall interaction within nanopores cannot be neglected, resulting in the no-slip boundary assumption being inadequate. For gas transport in nanopores, Knudsen diffusion should be considered, while slippage should be considered for liquid transport in nanopores. The intricate interplay of these key transport mechanisms within nanoporous media is comprehensively addressed through a species transport approach.

The proposed conceptual model relies on several assumptions to simplify the investigation of fluid transport in nanoporous media:

1. (i) the fluid is assumed to be a homogeneous distribution in an individual nanopore;

2. (ii) the high permeability porous medium is considered to have infinite conductivity when compared with the ultralow permeability of the nanoporous medium;

3. (iii) the nanoporous medium is viewed as a collection of parallel capillaries;

4. (iv) chemical reactions or adsorption phenomena are ignored;

5. (v) the system is assumed to be isothermal during codiffusion and counterdiffusion processes;

6. (vi) the fluid is not subjected to external body forces such as gravity or electric forces.

2.3. Steady-state, single-phase multicomponent fluid transport in a nanopore

This study accounts for key transport mechanisms, encompassing viscous flow, molecular diffusion, Knudsen diffusion in the gas phase and slippage in the liquid phase. Notably, the coupling of viscous flow and molecular diffusion is achieved through a species momentum balance equation. Additionally, Knudsen diffusion and slippage are incorporated via a slip boundary condition at the pore wall.

We focused on the species momentum balance instead of the bulk mixture momentum balance. Multicomponent fluid transport through an individual nanopore is given by a microscopic species momentum balance equation (Kerkhof & Geboers 2005b), equivalently, the generalized MS equation (Bhatia et al. 2011),

(2.2)\begin{equation}{\rho_i}\frac{{\partial {\boldsymbol{v}_i}}}{{\partial t}} ={-} {\rho _i}({\boldsymbol{v}_i}\boldsymbol{\cdot}\boldsymbol{\nabla }{\boldsymbol{v}_i}) - {c_i}{\nabla_T}{\mu _i} + {\rho _i}{\tilde{\boldsymbol{F}}_i} + {c_t}RT\sum_{\scriptstyle j = 1 \atop \scriptstyle j {\ne} i }^n {\frac{{{x_i}{x_j}}}{{{D_{ij}}}}({\boldsymbol{v}_j} - {\boldsymbol{v}_i})} + \boldsymbol{\nabla}\boldsymbol{\cdot }(2{\eta_i}{\boldsymbol{S}_i}),\end{equation}

where ${\rho _i}$ is the mass density of component i, ${\boldsymbol{v}_i}$ is the velocity of component i, ${c_i}$ is the molar concentration, ${c_t}$ is the total molar concentration, ${\mu _i}$ is the chemical potential of component i, ${\boldsymbol{\nabla }_T}{\mu _i}$ is the gradient of the chemical potential of component i at constant temperature (T), $\tilde{\boldsymbol{F}}$ is the external body force, R is the universal gas constant, ${x_i}$ is the mole fraction of component i, ${D_{ij}}$ is the MS diffusion coefficient, ${\eta _i}$ is the partial viscosity of component i and ${\boldsymbol{S}_i}$ is the rate of deformation tensor.

In the absence of a radial chemical potential gradient and external body force under steady-state conditions, the liquid transport equation for a capillary can be simplified from (2.2) and expressed as the following equation with four boundary conditions (Kerkhof, Geboers & Ptasinski 2001; Bhatia et al. 2011):

(2.3a)\begin{gather}-{c_i}{\left. {\frac{{\textrm{d}{\mu_i}}}{{\textrm{d}\zeta}}} \right|_T} + {\eta _i}\frac{1}{r}\frac{\partial}{{\partial r}}\left( {r\frac{{\partial {v_{i,\zeta}}}}{{\partial r}}} \right) + {c_t}RT\sum_{\scriptstyle j = 1 \atop \scriptstyle j {\ne} i }^{{n_c}}{\frac{{{x_i}{x_j}}}{{{D_{ij}}}}({v_{j,\zeta }} - {v_{i,\zeta }})} = 0,\end{gather}

(2.3b)\begin{gather}{v_{i,\zeta }}(r = R) ={-} G\frac{{\textrm{d}{v_{i,\zeta }}}}{{\textrm{d}r}},\end{gather}

(2.3c)\begin{gather}\frac{{\partial {v_{i,\zeta }}}}{{\partial r}}(r = 0) = 0,\end{gather}

(2.3d)\begin{gather}{\mu _i}(\zeta = 0) = {\mu _i}(\kern 1.5pt {p_{in}},{x_{i,in,i = 1,2, \ldots ,{n_c}}}),\end{gather}

(2.3e)\begin{gather}{\mu _i}(\zeta = L) = {\mu _i}(\kern 1.5pt {p_{out}},{x_{i,out,i = 1,2, \ldots ,{n_c}}}),\end{gather}

where G is the slip modulus, $\zeta $ represents the distance along the capillary direction, r is the distance along the radius direction, L is the length of the capillary and R is the radius of the capillary. Under isothermal conditions, the chemical potential is calculated based on the Peng–Robinson equation of state (PR-EOS) with a given pressure and fluid composition. A constant slip modulus is considered to account for the fast mass transport of liquid in the nanopores. For the gas phase, Knudsen diffusion is considered with a Maxwell slip boundary condition, and the slip modulus can be expressed as (Kerkhof & Geboers 2005b)

(2.4)\begin{equation}G = 1.19\frac{\eta }{p}{\left( {\frac{{8RT}}{{{\rm \pi} M}}} \right)^{1/2}}.\end{equation}

As the last term on the left-hand side of (2.3a) depends on the velocities of other components, it is necessary to solve all the equations simultaneously. These equations can be solved using a numerical method described in Appendix A.

2.4. Unsteady state, two-phase multicomponent mixture transport in a nanoporous medium

2.4.1. Governing equations

We extend our investigation to solve the pressure and composition field within a 1-D nanoporous medium by integrating transport equations with mass conservation equations. Here, we present the governing equations for a nanoporous medium, treating the entire high permeability porous media as a boundary condition. The mass conservation equation is applied to the liquid and vapour phases. As a result, there are n_c mass conservation equations with ${n_c}$ components in both mediums. The governing equations for the nanoporous medium can be expressed as (Deen 2013)

(2.5)\begin{equation} \frac{\partial \left[ {\phi {\displaystyle\sum_\beta} {({c_\beta }{S_\beta }{x_{\beta ,i}})} } \right]}{{\partial t}} ={-} \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{N}_i}\quad \textrm{for } i = 1,2, \ldots ,{n_c},\end{equation}

where the subscript $\beta $ is the phase index, which represents the liquid phase ($\beta = l$) or vapour phase ($\beta = v$), ${\boldsymbol{N}_i}$ is the molar flux of component i (mole per unit area per unit time), $\phi $ is the porosity, ${x_{\beta ,i}}$ is the mole fractions of component i in the $\beta $ phase, ${c_\beta }$ is the molar concentration of the $\beta $ phase (mole per unit volume), ${S_\beta }$ is the saturation of the $\beta $ phase and ${n_c}$ is the total number of components.

The molar flux within the nanoporous medium is determined by utilizing the binary friction model (BFM), which is explained in detail in the subsequent section.

2.4.2. The BFM

Assuming the nanoporous medium is comprised of a bundle of capillaries with uniform radii, the area-averaged molar velocity across an individual capillary is the same as that of the entire nanoporous medium. The molar flux across an interface is computed by multiplying the area-averaged molar velocity by the interface area. The velocity profile within a single capillary can be calculated by solving a microscopic species momentum balance equation, as shown in § 2.3. Accordingly, the area-averaged molar velocity across an individual capillary is computed by integrating the velocity across interactions along the radial direction and subsequently dividing by the area,

(2.6)\begin{equation}\langle {v_i}\rangle = \frac{2}{{{R^2}}}\int_0^R {{v_i}(r)r\,\textrm{d}r} ,\end{equation}

where $\langle {v_i}\rangle $ is the area-averaged molar velocity of component i.

Applying the aforementioned methodology, the pressure and composition fields are solved simultaneously with the velocity field across multiple interactions. However, the velocity field is not necessarily solved if the species momentum balance equation can be simplified to directly relate flux with pressure and composition. Fortunately, Kerkhof (1996) proposed the BFM, grounded on the assumption of ‘homogeneous viscous mixture flow’, offering such a method to compute the average velocity and flux from the species momentum balance equation. Kerkhof & Geboers (2005b) conducted a comparison between the BFM and (2.6), alongside multicomponent experiments, revealing generally insignificant discrepancies. As the BFM is more general for multicomponent and simplify the calculation process, it is used to calculate the molar flux within the matrix and through the matrix fracture interface in this paper.

The constitutive equation that relates flux to nanoporous medium properties for ${n_c}$ components can be expressed by (Kerkhof (1996); Pant et al. (2013))

(2.7)\begin{equation}{-} \frac{{{x_i}}}{{RT}}{\nabla _T}{\mu _i} = \sum_{\scriptstyle j = 1 \atop \scriptstyle j {\ne} i }^{{n_c}} {\frac{{{x_j}{\boldsymbol{N}_i} - {x_i}{\boldsymbol{N}_j}}}{{{c_t}D_{ij}^e}}} + {\boldsymbol{f}_{im}}{\boldsymbol{N}_i}\quad \textrm{for}\;i = 1,2, \ldots ,{n_c},\end{equation}

where ${\boldsymbol{f}_{im}}$ is the modified wall friction factor, which is a function of the partial viscosities and the slip modulus. The chemical potential gradients at constant temperature are expanded to explicitly include the contribution of pressure and composition in previous works (Krishna & Wesselingh 1997; Ma & Emami-Meybodi 2024b),

(2.8)\begin{equation}{\nabla _T}{\mu _i} = \sum_{j = 1}^{{n_c} - 1} {{\varGamma _{ij}}\boldsymbol{\nabla }{x_i}} + {\bar{V}_i}\boldsymbol{\nabla }p,\end{equation}

where ${\varGamma _{ij}}$ is thermodynamic factors and ${\bar{V}_i}$ is partial molar volume. However, we do not expand the chemical potential gradient in this work because it is more convenient to calculate the flow rate under two-phase conditions.

The modified wall-friction factor is used to describe the interaction between species and pore walls, taking into account both viscous and wall–molecule collision effects. When dealing with non-ideal gas, weighted coefficients are used to account for both the friction from intermolecular collisions and the friction resulting from collisions between molecules and pore walls. The ratio of the frequency of molecule-wall collisions to the total collision frequency and the ratio of intermolecular collision frequency to total collision frequency are used to establish the weighting coefficients (Wu et al. 2015; Pang, Fan & Chen 2021; Liu & Emami-Meybodi 2023). As a result, the modified wall-friction factor can be expressed as follows:

(2.9)\begin{equation}{\boldsymbol{f}_{im}} = \frac{1}{{{c_t}}}{\left( {{\omega_K}D_{iK}^e + {\omega_v}\frac{\boldsymbol{K}}{{{\kappa_i}}}} \right)^{ - 1}}\end{equation}

where ${\omega _K}$ and ${\omega _v}$ are Knudsen number dependent weight factors for Knudsen diffusion and viscous flow, respectively. Due to the small Knudsen number, there is no Knudsen diffusion in the liquid phase; hence ${\omega _K} = 0$, ${\omega _v} = 1$. Here $\boldsymbol{K}$ represents permeability and ${\kappa _i}$ represents the fractional viscosity coefficient.

The effective diffusion coefficient is adjusted by saturation, and the phase permeability is adjusted by relative permeability, ${k_{r\beta }}$, to account for multiphase fluid transport in porous media. We use square brackets to represent the variable as a matrix or vector, wherein each element corresponds to an individual component or pair of components. Equation (2.7) can be cast into ${n_c}$-dimensional matrix notation to obtain the following explicit expression for the flux in each phase $\beta $:

(2.10a)\begin{gather} {[}{\boldsymbol{N}_\beta }] ={-} [{\boldsymbol{T}_\beta }][{\nabla _T}\mu ],\end{gather}

(2.10b)\begin{gather} {[}{\boldsymbol{T}_\beta }] = \frac{{{c_\beta }}}{{RT}}{[{\boldsymbol{B}_\beta }]^{ - 1}}[{\boldsymbol{X}_\beta }],\end{gather}

where $[{\boldsymbol{T}_\beta }]$ represents transmissibility matrix of components, $[{\boldsymbol{X}_\beta }]$ is a diagonal matrix with mole fraction as non-zeros value at diagonal position, $[{\nabla _T}\mu ]$, are the gradient of the chemical potential vector, $[{\boldsymbol{B}_\beta }]$ is a ${n_c} \times {n_c}$ matrix,

(2.11a)\begin{gather}{\boldsymbol{B}_{\beta ,ij}} = \sum_{\scriptstyle k = 1 \atop \scriptstyle k {\ne} i }^{{n_c}} {\frac{{{x_{\beta ,k}}}}{{D_{\beta ,ik}^e}}} + \frac{1}{{{\omega _K}D_{\beta ,iK}^e + {\omega _v}\dfrac{{\boldsymbol{K}{k_{r\beta }}}}{{{\kappa _{\beta ,i}}}}}}\quad \textrm{for}\;i = j,\end{gather}

(2.11b)\begin{gather}{\boldsymbol{B}_{\beta ,ij}} ={-} \frac{{{x_{\beta ,i}}}}{{D_{\beta ,ij}^e}}\quad \textrm{for}\;i \ne j.\end{gather}

Several parameters are needed to calculate transmissibilities, including MS diffusivities, permeance, partial viscosity, Knudsen diffusivity and weight factors, which can be obtained using the method delineated in Appendix B.

The final governing equation is obtained by expanding (2.10) for each individual component and subsequently substituting into (2.5):

(2.12)\begin{equation}\frac{\partial \left[ {\phi {\displaystyle\sum_\beta} {({c_\beta }{S_\beta }{x_{\beta ,i}})} } \right]}{\partial t} ={-} \boldsymbol{\nabla }\boldsymbol{\cdot }\left\{ {\sum_\beta {\sum_{j = 1}^{{n_c}} {{\boldsymbol{T}_{\beta ,ij}}{\nabla_T}{\mu_j}} } } \right\}.\end{equation}

The nanoporous medium is discretized into multiple grids to solve the pressure and composition fields, with the assumption of local thermodynamic equilibrium within each grid. In each grid, ${n_c}$ independent variables can be solved from ${n_c}$ mass conservation equations (2.12) for all components. We note that the pressure ($\,p$) and overall composition (${z_1}, \ldots ,{z_{{n_c}}}$) are the primary variables that need to be solved. Meanwhile, chemical potential serves as an intermediate (secondary) variable that facilitates the computation of flux and the construction of governing equations. In (2.12), the chemical potential gradient is calculated through

(2.13)\begin{equation}{\nabla _T}{\mu _i} = RT{\nabla _T}\,\textrm{ln}\,{f_i},\end{equation}

where f is the fugacity. Under local thermodynamic equilibrium conditions, the fugacity is calculated based on the PR-EOS with given pressure and fluid composition (Peng & Robinson 1976; Michelsen & Mollerup 2004),

(2.14)\begin{equation}{f_i} = {x_i}p\,\textrm{exp}\left[ {(Z - 1) - \textrm{ln}(Z - B) + \frac{A}{{2B\sqrt 2 }}\,\ln \frac{{Z + (1 - \sqrt 2 )B}}{{Z + (1 + \sqrt 2 )B}}} \right],\end{equation}

where Z is the compressibility factor solved from PR-EOS, A is the PR-EOS energy parameter and B is the PR-EOS covolume parameter.

2.4.3. Initial and boundary conditions

The boundary conditions of our system exist at the interface of the nanoporous medium- and high-permeability porous medium (i.e. $\zeta = 0$), where we specify the inner boundary condition, and at the outer end of nanoporous medium, where a closed outer boundary condition is specified. For the 1-D system under consideration, the following initial and boundary conditions are required to solve (2.12):

(2.15a)\begin{gather}p(\zeta ,0) = {p_0},\end{gather}

(2.15b)\begin{gather}{z_i}(\zeta ,0) = {z_{i0}}\quad \textrm{for}\;i = 1,2, \ldots ,{n_c} - 1,\end{gather}

(2.15c)\begin{gather}p(0,t) = {p_{HP}}(t),\end{gather}

(2.15d)\begin{gather}{z_i}(0,t) = {z_{HP,i}}(t)\quad \textrm{for}\;i = 1,2, \ldots ,{n_c} - 1.\end{gather}

(2.15e)\begin{gather}\frac{{\partial p}}{{\partial \zeta }}(L,t) = 0,\end{gather}

(2.15f)\begin{gather}\frac{{\partial {z_i}}}{{\partial \zeta }}(L,t) = 0\quad \textrm{for}\;i = 1,2, \ldots ,{n_c} - 1,\end{gather}

where ${p_0}$ is the initial pressure, ${z_{i0}}$ is initial overall mole fraction. The outer boundary conditions are based on the assumption that the matrix block has a no-flow boundary at one end.

The entire high-permeability porous medium is treated as the inner boundary, where ${p_{HP}}(t)$ and ${z_{HP,i}}(t)$ represent the time-varying pressure and overall mole fraction within the high-permeability porous medium. In the case of codiffusion, we assume a constant pressure of high permeability porous media, ${p_{HP}}(t) = {p_{HP,co}} < {p_0}$, and a portion of the fluid is removed from the high-permeability porous media with the same fluid composition, ${z_{q,i}}(t) = {z_{HP,i}}(t)$. Similarly, in the case of counterdiffusion, we assume a constant pressure of high permeability porous media, ${p_{HP}}(t) = {p_{HP,counter}} > {p_0}$, and a certain amount of fluid will be added to the high-permeability porous media. The overall composition of the added fluid is held constant, specified as ${z_{q,i}}(t) = {z_{I,i}}$. The fluid overall composition within high-permeability porous media varies with time during codiffusion and counterdiffusion processes. A material balance equation for the high-permeability porous medium is solved together with the governing equations, (2.12), to obtain the fluid composition within the high-permeability porous medium at a given time. Further information regarding the inner boundary conditions can be found in Appendix  C.

3. Results

3.1. Steady-state, single-phase multicomponent fluid transport in a nanopore

We investigate the species velocity profiles of a ternary alkane mixture – methane (CH₄), n-butane (C₄H₁₀) and n-dodecane (C₁₂H₂₆) – within a nanopore under steady-state conditions by solving (2.3). For simplicity, we use C1, C4 and C12 to refer to CH₄, C₄H₁₀ and C₁₂H₂₆, respectively. The chosen alkanes represent light (C1), intermediate (C4) and heavy (C12) alkane species; hence, the impacts of molecule types and their diffusivities on transport can be studied. Figure 3 depicts the effects of pressure gradients (0, 10 and 100 MPa m⁻¹) at different pore radii of 1, 5 and 10 on the velocity profiles. The velocity values are normalized by the velocity of C1 at the centre of nanopores, i.e. v ₀ = v _C1 (r = 0). The mass-averaged velocity is calculated based on all the species’ velocities to represent the advection (bulk fluid transport). Considering the mass-averaged velocity as a reference velocity, the difference between species velocity and mass-averaged velocity represents diffusion.

Figure 3. Species velocity profiles in a nanopore. Species and mass-averaged velocity profiles are calculated at pressure gradients of 0, 10 and 100 MPa m⁻¹ and pore radii of 1, 5 and 10 nm. The sign of the velocity represents the transport direction, with ‘+’ indicating the same direction as the pressure gradient and ‘−’ indicating the opposite direction. The velocity values are normalized by the velocity of CH₄ at the centre of nanopores, i.e. ${v_0} = {v_{C{H_4}}}$ (r = 0). The distance along the pore diameter is normalized with the pore radius.

Figure 3 shows that in the absence of a pressure gradient ($\boldsymbol{\nabla }p = 0\;\textrm{MPa}\;{\textrm{m}^{ - 1}}$) where the mass-averaged velocity is zero, species transport occurs due to composition gradients. Based on the Fickian diffusion theory, the species velocity profiles should remain fixed along the pore diameter independent of pore size. However, the results in figure 3 demonstrate that wall friction affects diffusion and is responsible for the parabolic velocity profile behaviour, especially in the nanopore with R = 1 nm. Furthermore, diffusion is significantly hindered due to fluid molecules and pore wall interactions as pore size decreases from R = 10 to 1 nm. The wall friction effect is limited to the near-wall region at larger pores. In the case involving a nanopore with a radius of 10 nm under a pressure gradient of 0 MPa m⁻¹, a virtually flat velocity profile is observed, with only a small region near the wall exhibiting a velocity gradient. Another notable observation for $\boldsymbol{\nabla }p = 0\;\textrm{MPa}\;{\textrm{m}^{ - 1}}$ cases is that the velocity at the centre of nanopores does not increase with pore size due to intermolecular friction, i.e. fluid molecule interactions. However, in the presence of a pressure gradient, $\boldsymbol{\nabla }p > 0$, the velocity at the centre of the pore is found to increase with an increase in pore size under the same pressure gradient, consistent with Poiseuille's law (Sutera & Skalak 1993).

As mentioned earlier, the difference between the species velocity and the mass-averaged velocity reflects the diffusion contribution. The results presented in figure 3 indicate that the difference between each component's velocity and the mass-averaged velocity is more pronounced in the 1 nm pore compared with the 10 nm pore under the same pressure gradient, demonstrating the significant contribution of diffusion in smaller pores. Furthermore, the magnitude between the advective velocity and each component's velocity increases as the pressure gradient decreases, signifying an increase in diffusion contribution with pressure gradient reduction.

The contribution of diffusion is quantitatively assessed through the multicomponent Sherwood number (Sh), which is the ratio of total to diffusive mass transfer rate (see Appendix D). Unlike the Sherwood number traditionally employed in dilute solutions or ideal gas mass transport systems, the Sherwood and Péclet numbers in the present work are for each individual component within multicomponent mass transport systems. A higher Sh means a higher advection contribution to the total flux. Figure 4(a) shows the Sh as a function of the pressure gradient for nanopores with R = 1, 5, 10, 25 and 50 nm. As expected, Sh increases as the pressure difference increases because the pressure difference predominantly governs advection. Therefore, there is a greater advective flow and advection contribution. Accordingly, the diffusion contribution of macropores (pore diameter > 50 nm) is negligible, even under small pressure gradient conditions. On the contrary, the contribution of diffusion remains significantly high at small pore sizes, even at high-pressure gradients.

Figure 4. Effect of pressure gradient and pore size on diffusive transport contribution. Effect of (a) pressure gradient and (b) mole fraction gradient on Sherwood number at different pore sizes.

Figure 4(b) shows the Sh as a function of the mole fraction gradient for nanopores with R = 1, 5 and 10 nm under a pressure gradient of 10 MPa m⁻¹. The sign of Sh conveys the relationship between total and diffusive flux direction, with ‘+’ indicating they flow in the same direction and ‘−’ indicating the opposite direction flow. While Sh always has a positive value during the species cocurrent transport that occurs when concentration (mole fraction) and pressure gradients are in the same direction, Sh can be either positive or negative during the species countercurrent transport when a pressure gradient opposes the concentration gradient. For instance, in a 5 nm pore, when the composition gradient of C1 reaches $\boldsymbol{\nabla }{x_{C1}} ={-} 10\;{\textrm{m}^{ - 1}}$, the concentration diffusion supersedes the pressure diffusion, aligning the overall flux with the diffusive direction, thus imparting a positive value to Sh. As the composition gradient decreases to the point where diffusive and advective fluxes balance, the overall velocity of C1 approaches zero, i.e. Sh → 0. Further reduction in the composition gradient causes the advective flux to surpass the diffusive flux, prompting the overall flux in the opposite direction as diffusive flux, leading to a negative Sh value. We note that the diffusive transport rate, which resides in the Sh denominator, decreases as the composition gradient decreases, resulting in the absolute value of Sh on an upward trajectory. Additionally, at the same composition gradient, the total flux significantly increases with increasing pore size due to increased advection, causing Sh to increase with pore size. Notably, the value of Sh tends to infinity when the mole fraction gradient equals zero, i.e. zero diffusive transport rate. We note both species cocurrent and countercurrent transport occur simultaneously in the multicomponent fluid transport system, with some species experiencing cocurrent transport while others experience countercurrent transport. We consider the absolute value of Sh in the following analysis and investigation.

Besides pressure gradient and mole fraction gradient, pore size, diffusion coefficient, partial viscosity and species type (alkane length) influence species transport within nanopores. Péclet number (Pe), the ratio of advective to diffusive rates, is typically used to consider all influencing parameters comprehensively (Patankar 1980). We extend the definition of Pe to multicomponent fluid systems by calculating the diffusive transport rate with the MS approach (see the Methods section). Dimensional analysis then shows that there is an implicit function between Sh and Pe for any component in multicomponent systems,

(3.1)\begin{equation}S{h_i} = f(P{e_i},{\varphi _i})\quad \textrm{for}\;i = 1,2, \ldots ,{n_c},\end{equation}

where ${n_c}$ is the total number of components and φ_i is the fluid–solid friction modulus calculated from the ratio between diffusive and shear stress friction coefficient (see Appendix D).

Figure 5(a) shows that the Sherwood number for all parameters (pore size, pressure gradient, diffusion coefficients, species type) collapses on a master curve with dependence on the Péclet number under a fluid–solid friction modulus of 100. The Sh–Pe plot shows two distinctive curves corresponding to species cocurrent and countercurrent transport processes. During the species cocurrent transport process, where diffusive and advective transports align in the same direction, Sh increases in parallel with an upsurge in Pe (see blue cycles in figure 5a). When Pe reaches a significant magnitude, a linear relationship between Sh and Pe emerges. During the species countercurrent transport process, where diffusive and advective transports are in opposite directions, Sh decreases as Pe increases until Sh → 0, where advective transport counterbalances the diffusive transport. However, Sh starts increasing as Pe increases, and ultimately, Sh–Pe curves for both cocurrent and countercurrent transport processes collapse onto one straight line.

Figure 5. Relationship between Sherwood number and Péclet number. (a) Relationship between Sherwood number and Péclet number with fluid–solid friction modulus of 100. The relationship between Sh and Pe illustrates three distinct fluid transport systems: diffusion-dominated; transition; advection-dominated. The dashed lines represent two limits of diffusion-dominated and advection-dominated systems given by (3.2a) and (3.2c), respectively. (b) Effect of fluid–solid friction module on the relationship between Sh and Pe.

We adopt the scaling relation of Sh ~ Pe^m to scale the diffusion contribution for the systems with various parameters (Marle 1981). Figure 5(a) illustrates the correlation between the Sherwood and Péclet numbers, revealing three distinct fluid transport systems: diffusion-dominated; transition; advection-dominated. The asymptotic scaling relations can be expressed as follows:

(3.2a)\begin{gather}S{h_i} \approx \frac{1}{{{\omega _i}}},\quad P{e_i} < 0.1,\end{gather}

(3.2b)\begin{gather}S{h_i} = 0,\quad P{e_i} \approx \frac{1}{{{\omega _i}}},\end{gather}

(3.2c)\begin{gather}S{h_i} \approx P{e_i},\quad P{e_i} > 100,\end{gather}

where ω_i is the average mass fraction of component i.

Equations (3.2a) and (3.2c) represent the diffusion-dominated and advection-dominated systems, respectively. Notably, the relationship between Sh and Pe exhibits two prominent limits at Pe_i = 0.1 and Pe_i = 100. In transition systems with 0.1 < Pe_i < 100, diffusion and advection play significant roles. Within this Pe range, (3.2b) represents a diffusion–advection neutral point (Sh_i = 0) for the species countercurrent transport process that occurs when diffusive and advective fluxes are equal but in opposing directions.

The fluid–solid friction modulus, φ, serves as a conduit to the realm of shear stress impact on diffusive transport. Figure 5(b) examines the influence of φ on Sh as a function of Pe for a three-component (ternary) fluid system. High values of φ (φ_i > 100) represent a system in which shear stress has negligible impact on diffusive transport, and if the system is also diffusion dominant, then Sh equals to 1/ω_i. However, as φ decreases, Sh gradually diminishes since the effect of shear stress on diffusive transport increases. In transition systems, the pore wall effect hinders the diffusive flux, leading to the diffusion–advection neutral point decreasing as φ decreases. On the other hand, in advection-dominated systems, Sh and Pe consistently maintain a linear relationship. Broadly, as the pore wall effect increases, meaning that φ decreases, the Sh–Pe curve moves down and towards the left. This result is consistent with the phenomenon reflected in species velocity profiles (figure 3). The fluid–solid friction modulus (φ) reflects the influence of shear stress on diffusive transport, which is intricately proportional to the square of the pore radius. Diffusion is substantially constrained as the pore size decreases, a consequence of the shear stress encountered due to the presence of pore walls. Conversely, in larger pores, the influence of fluid–solid interaction is confined to a narrow near-wall region, thereby exerting limited impact.

3.2. Unsteady-state, partially miscible two-phase multicomponent fluid transport in the nanoporous medium

In this section, our analysis focuses on codiffusion and counterdiffusion processes for ternary alkane mixtures and CO₂ under unsteady-state high pressure and elevated temperature conditions. An equal molar composition ${z_0} = [0.34,0.33,0.33]$ of C1, C4 and C12 is selected to represent the initial ternary alkane mixture in the nanoporous medium. For codiffusion, the nanoporous medium is initially saturated with the ternary alkane mixture at ${p_0} = 34.47\;\textrm{MPa}$ (5000 p.s.i.) connected to a high-permeability porous medium with a constant pressure of ${p_{HP}}(t) = 6.89\;\textrm{MPa}$ (1000 p.s.i.). The ternary alkane mixture is produced from nanopores, where all species move in the same direction. For counterdiffusion, the nanoporous medium is initially saturated with the ternary mixture at an initial pressure of ${p_0} = 6.89\;\textrm{MPa}$, while CO₂ is introduced to the high-permeability porous medium at ${p_{HP}}(t) = 34.47\;\textrm{MPa}$. Different species may move in the opposite direction due to the composition gradient. Note that a vapour phase appears when the pressure is below the saturation pressure (10.68 MPa), and CO₂ exhibits partial miscibility with the alkane mixture. The PR-EOS (Peng & Robinson 1976) is employed to calculate the phase behaviour of the alkane and CO₂ mixture.

Figure 6(a,b) shows species’ molar rates and recovery factors with time, respectively, for the codiffusion process. The recovery factor for each species is calculated from the ratio of the cumulative produced molar flux to the initial moles of the species in the nanoporous medium. The composition of produced fluid from the nanoporous medium differs from the initial equimolar composition, especially when the vapour phase appears. The species flux rates remain the same during the early codiffusion when liquid remains a single phase. As the pressure drops in the nanoporous medium increases, a vapour phase emerges within the porous medium, resulting in changes in the species fluxes (see figure 6a). The species flux is inversely proportional to its molecular weight: the light alkane species (C1) has the highest molar rate, then the intermediate species (C4) and the heavy species (C12) has the lowest. These findings are consistent with the molecular simulations conducted by Falk et al. (2015) and Perez & Devegowda (2020) and the field data reported by Freeman et al. (2012). Similar trends are observed for the recovery factors, and the C1 recovery factor is 15 % after 10 years of simulation, which is significantly higher than C12 with a recovery factor of 9 %.

Figure 6. Partially miscible two-phase multicomponent transport during the codiffusion process. (a) Molar flow rates of C1, C4 and C12 with time during the codiffusion process. Species fluxes are the same at the early time, but they separate from each other as a vapour phase appears. (b) Species molar recovery with time during the codiffusion process. The significant separation between light component (C1) and intermediate (C4) and heavy (C12) species is revealed, with higher recovery for C1 than for C4 and C12. (c) Schematic diagram of the interplay of partially miscible vapour and liquid phase transport phenomena.

The simulation results, including the species fluxes, reveal that the interplay of the appearance of the vapour phase due to pressure drop, mass transfer across phases, partial miscibility of species and multicomponent transport mechanisms within the liquid and vapour phase dictates the species mass transport within the nanoporous medium. Figure 6(c) depicts how the interplay of these phenomena results in the higher production of lighter species. Under the initial equal molar composition in the liquid phase, the transmissibilities of all species are in the same order of magnitude. Hence, species fluxes are the same initially. The pressure within the nanoporous medium decreases as fluid is produced, causing light and intermediate alkane components to form a vapour phase (see figure 6c). Species transmissibilities in the vapour phase are higher than in their transmissibilities in the liquid phase due to the significant difference in their diffusion coefficients and viscosities. Moreover, the vapour phase primarily consists of light and intermediate components (C1 and C4), with a small amount of C12 because of low volatility. Accordingly, due to the higher flux in the vapour phase than in the liquid phase, light components are produced faster due to the higher species flux in the vapour phase. In contrast, heavy species remain less substantive as the flux of species in the liquid phase is slower. The vapour–liquid interface moves deeper into the nanoporous medium as the light species continuously vaporizes from the liquid phase into the vapour phase. This process forms a continuous vapour phase within the nanoporous medium, enhancing the recovery of light alkane species.

Figure 7(a,b) shows the pressure profile and CO₂ profile within the nanoporous medium after different time counterdiffusion, respectively. The mixing extent of the CO₂ and alkane mixture is an important factor in helping us to investigate the multicomponent transport mechanism during the counterdiffusion process. Here, we adopt the Waggoner, Castillo & Lake (1992) approach to define the ‘mixing region’, which is the distance over which the cross-sectional average of the introduced fluid composition decreases from 0.9 to 0.1 (see figure 7b). The length of the mixing region quantitatively describes the mixing extent between the CO₂ and alkane mixture. It is important to note that the fluid in the mixing region can be single-phase or two-phase, which differs from the miscible zone concept (Lake et al. 2014). As depicted in figure 7(b), in the case of counterdiffusion, CO₂ diffuses into the nanoporous medium and mixes with alkanes, creating a ‘mixing region’, which grows as CO₂ penetrates deeper into the nanoporous medium. Over time, the CO₂ mole fraction becomes very high in the near-fracture region, forming a CO₂ bank that pushes the mixing region deeper into the nanoporous medium. The ‘mixing region’ is mainly formed at the beginning of the counterdiffusion process with a significant pressure gradient. Even though all the species move in the same direction under the pressure gradient, the pressure gradient is still crucial in promoting the mixing of CO₂ and alkane. We calculated the distance of investigation using (3.3) to quantify the rate of pressure diffusion through the nanoporous medium (Behmanesh et al. 2015),

(3.3)\begin{equation}\textrm{DOI} = 0.194\sqrt {(K/\phi {\eta _I}{c_{tI}})t} ,\end{equation}

where ${\eta _I}$ and ${c_{tI}}$ are the mixture's viscosity and the total compressibility at the initial pressure, respectively.

Figure 7. Partially miscible two-phase multicomponent transport during the counterdiffusion process. (a) Pressure profiles along the porous medium during codiffusion after 1, 6, 24, 60 and 120 months. (b) The distance of investigation movement during the counterdiffusion process. The pressure diffusion approaches the boundary after four years of injection. (c) The CO₂ mole fraction profiles along the nanoporous medium. Three regions (CO₂ bank; mixing region; alkane mixture bank) are separated after 120 months of counterdiffusion based on the boundary defined by CO₂ mole fractions of 0.1 and 0.9. (d) Mixing region and CO₂ bank fronts movement during the counterdiffusion process. The mixing region length keeps increasing during the counterdiffusion process.

Figure 7(c,d) illustrate the extent of CO₂ and alkane mixing changing over time during counterdiffusion. The pressure diffusion process occurs much faster than changes in the compositions. As the distance of investigation approaches the boundary, the mixing region front advances slowly, indicating that the pressure difference is a crucial factor promoting CO₂ and hydrocarbon mixture mixing. The pressure gradient drives all components in the same direction, but the differing transmissibilities of each component cause them to have different velocities. This velocity difference drives the mixing between CO₂ and hydrocarbon mixture.

As the pressure diffusion approaches the boundary, the pressure gradient gradually decreases, decreasing the driving force. However, the mixing region still slowly expands due to the diffusion caused by the composition gradient. The diffusion caused by only the composition gradient is a slow process – see the CO₂ mole fraction profile between 60 and 120 months in figure 7(b). That is the reason why most of the CO₂ is still confined to the near high-permeability medium even after a decade.

During the counterdiffusion of CO₂ into the alkane mixture, the length of the mixing region increases. However, the increase rate in early time is much faster than that in late time. The early time curve represents the formation of the mixing region with both pressure and composition gradient as the driving force, while the late time curve represents the expansion of the mixing region after the solvent bank is formed with only the composition gradient as the driving force.

4. Discussion

The transmissibilities calculated from the method in this paper satisfy Onsager's reciprocal relations (Onsager 1931). In addition to experimental data, the species' properties can be accurately determined by matching molecular dynamics simulations. The method based on species movement in this paper proved a great bridge connecting macroscope multicomponent transport with molecular dynamic simulation since both are based on potential energies and molecular interactions. Different transport mechanisms can be lumped into the interaction (friction) between different molecules, the interaction between identical molecules and the interaction between molecules and walls. For example, Krishna & van Baten (2005) successfully lumped surface diffusion into the wall-friction coefficient by matching molecular dynamic simulation results to simulate multicomponent transport across zeolite membranes.

The species movement approach presented in this paper follows phenomenological law (Deen 2013; Firoozabadi 2016) and is a set of more general constitutive equations than Darcy's and Fick's laws. The proposed paradigm shifts in predicting multicomponent transport by species movement can reconcile different transport mechanisms when the pore size decreases to the nanoscale and provide a possibility to couple nanoconfined fluid phase behaviour (Ma & Emami-Meybodi 2024a,c) to better understand multicomponent transport problems in many natural and industrial processes.

The contribution of diffusion and advection is quantified by the $P{e_i} = \tilde{v}L/\sum\nolimits_{j = 1}^{n - 1} {D_{ij}^e\Delta {x_j}} $. A small $Pe$ indicates diffusion dominance, while a large $Pe$ signifies advection dominance. We assume the average velocity in nanopores can be determined using the Haggen–Poiseuille equation (Holt et al. 2006), where the average velocity is proportional to the square of pore diameter, $\tilde{v}\sim {d^2}$. In contrast, the effective diffusion coefficient is related to the porous medium by $D_{ij}^e = {Q_m}{D_{ij}} = (\phi /\tau ){D_{ij}}$, which is not directly dependent on pore diameter. As the pore diameter decreases, there might be slight reductions in porosity ($\phi $) and tortuosity ($\tau $); however, the magnitude is limited. In contrast, the average velocity decreases significantly (up to several magnitudes) with pore diameter decreasing. Consequently, viscous flow (advection) undergoes substantial hindrance, resulting in an increased contribution from diffusion. Knudsen diffusion or slippage is incorporated as a slip boundary condition with a slip length denoted as ${L_s} = 2\sqrt {RG} $. The value of slip length is typically in nanoscale (Holt et al. 2006; Whitby et al. 2008; Javadpour 2009) and dependent on the fluid and pore wall solid characteristics. This parameter is lumped into

(4.1)\begin{equation}{\varphi _i} = \frac{{{x_i}{C_t}RT(R_p^2 + 4RG)}}{{{\eta _i}}}\sum_{\scriptstyle j = 1 \atop \scriptstyle j {\ne} i }^n {\frac{{{x_j}}}{{{D_{ij}}}}} .\end{equation}

When the pore size is large, the influence of Knudsen diffusion and slippage can be disregarded. It is only when the pore size diminishes to a scale comparable to the slip length that Knudsen diffusion and slippage exert significant contributions.

5. Conclusions

We proposed a species-based modelling approach for fluid transport in nanoporous media and developed a numerical multiphase, multicomponent model. The model provides a compositional approach to accurately capture the composition and pressure changes within a nanoporous medium. Our investigation of multicomponent transport in individual nanopores showed that wall friction affects the diffusive behaviour of species, particularly at small pores. The scaling relationship between Sh and Pe illustrates three distinct regimes: diffusion-dominated regime (Pe < 0.1); transition regimes (0.1 < Pe < 100); and advection-dominated regime (Pe > 100). In the case of codiffusion, light alkane species are produced more than the heavy alkanes when the gas phase appears because the mobility of light components in the vapour phase is higher than that of heavy components in the liquid phase. In the case of counterdiffusion, the pressure gradient remains a crucial factor in promoting the mixing of CO₂ and hydrocarbon mixture. The pressure gradient drives all components in the same direction, but the differing transmissibilities of each component cause them to have different velocities, driving the mixing between CO₂ and the alkanes. The composition gradient is another driving force for ‘mixing region’ expansion when the pressure gradient decreases. However, even after a decade, most injected CO₂ remains confined to the near high-permeability medium due to the slow species counterdiffusion rate driven by the composition gradient.