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 Reference Young2007; Yuan & Sundén Reference Yuan and Sundén2014), designing separation and catalytic processes (Krishna Reference Krishna2012), enhancing hydrocarbon production from shale reservoirs (Falk et al. Reference Falk, Coasne, Pellenq, Ulm and Bocquet2015), geo-storage of carbon dioxide (CO2) in shale formations (Lee, Cho & Lee Reference Lee, Cho and Lee2020; Emami-Meybodi et al. Reference Emami-Meybodi, Ma, Zhang, Rui, Rezaeyan, Ghanizadeh, Hamdi and Clarkson2025) and long-term integrity of CO2-reservoir caprocks (Kampman et al. Reference Kampman, Busch, Bertier, Snippe, Hangx, Pipich, Di, Rother, Harrington and Evans2016; Adu-Gyamfi et al. Reference Adu-Gyamfi, Ampomah, Tu, Sun, Erzuah and Acheampong2022). 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. Reference Fan, Schlick, Umbanhowar, Ottino and Lueptow2014). 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, Reference Darcy1856). Various Darcy-based transport models with apparent permeability (Darabi et al. Reference Darabi, Ettehad, Javadpour and Sepehrnoori2012; Wang, Wang & Chen Reference Wang, Wang and Chen2018) have been developed to simulate multicomponent fluid transport in nanoporous media, particularly with the applications to hydrocarbon production from shale reservoirs and CO2 injection into these formations (Chen, Balhoff & Mohanty Reference Chen, Balhoff and Mohanty2014; Ozowe, Zheng & Sharma Reference Ozowe, Zheng and Sharma2020). 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. Reference Falk, Coasne, Pellenq, Ulm and Bocquet2015).
Table 1. The constitutive equations used to express molar flux, Ni, 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 Reference Cronin, Emami-Meybodi and Johns2018; Liu & Emami-Meybodi Reference Liu and Emami-Meybodi2021). 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 (Reference Fick1855), inspired by the experiments of Graham (Reference Graham1850) 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 Reference Truesdell1962; Valdes-Parada, Wood & Whitaker Reference Valdes-Parada, Wood and Whitaker2023), kinetic theory (Kerkhof & Geboers Reference Kerkhof and Geboers2005a) and the thermodynamic approach (Curtiss & Bird Reference Curtiss and Bird1999; Firoozabadi Reference Firoozabadi2016). In nanoporous media, diffusion-based compositional models using Fickian-type modelling have been proposed (Cronin, Emami-Meybodi & Johns Reference Cronin, Emami-Meybodi and Johns2021) 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 Reference Liu and Emami-Meybodi2022) 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. Reference Moortgat, Firoozabadi, Li and Espósito2013; Zidane & Firoozabadi Reference Zidane and Firoozabadi2021). 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 Reference Krishna and Wesselingh1997; Deen Reference Deen2013). Therefore, using Fick's law may not suffice to describe multicomponent mass transport in nanoporous media (Krishna & Wesselingh Reference Krishna and Wesselingh1997; Corral-Casas et al. Reference Corral-Casas, Li, Borg and Gibelli2022).
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 (Reference Maxwell1860), 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 Reference Bird and Klingenberg2013). 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 Reference Curtiss and Bird1999) and continuum mechanics (Whitaker Reference Whitaker2009). 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 Reference Taylor and Krishna1993; Moortgat & Firoozabadi Reference Moortgat and Firoozabadi2013). 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. Reference Tian, Zhang, Lei, Yin, Kazemi and Wu2021). 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 Reference Do1998; Gruener & Huber Reference Gruener and Huber2008; Kärger et al. Reference Kärger, Ruthven and Theodorou2012). 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 Reference Knudsen1909). 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 (Reference Kramers and Kistemaker1943) in connection with equimolar counterdiffusion. Subsequently, it has been applied to investigate multicomponent gas flow in capillaries and porous media (Young & Todd Reference Young and Todd2005; Mills & Chang Reference Mills and Chang2013). 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 Reference Krishna and Wesselingh1997; Mason Reference Mason1983). 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 Reference Krishna1987; Krishna & Wesselingh Reference Krishna and Wesselingh1997; Mason Reference Mason1983). The DGM has succeeded in fields such as catalysis, membrane (Gao et al. Reference Gao, Chen, Yu and Asumana2011) and fuel cell modelling (Yuan & Sundén Reference Yuan and Sundén2014), 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 Reference Kerkhof1996; Bhatia, Bonilla & Nicholson Reference Bhatia, Bonilla and Nicholson2011; Pant, Mitra & Secanell Reference Pant, Mitra and Secanell2013). Recent studies have implemented DGM to study fluid transport in nanoporous media (Achour & Okuno Reference Achour and Okuno2022; Sena Santiago & Kantzas Reference Sena Santiago and Kantzas2022) 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 Reference Deen2013). This study employs the molar form of the species transport and continuity equation. The species molar flux can be expressed as
\begin{equation}{N_i} = {x_i}{c_t}\bar{v} + {J_i},\end{equation}
where Ni represents the molar flux of species i, ct is the total molar concentration of the mixture, xi is the mole fraction in the mixture, 
$\bar{v}$ is the molar average velocity and Ji 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. Reference Falk, Coasne, Pellenq, Ulm and Bocquet2015; Kavokine, Netz & Bocquet Reference Kavokine, Netz and Bocquet2021).
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. Reference Holt, Park, Wang, Stadermann, Artyukhin, Grigoropoulos, Noy and Bakajin2006; Wang, Javadpour & Feng Reference Wang, Javadpour and Feng2016). 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. Reference Pant, Mitra and Secanell2013), bulk tanks utilized in separation and catalytic devices (Veldsink, Versteeg & van Swaaij Reference Veldsink, Versteeg and van Swaaij1994) and fractures encountered within shale reservoirs (Hoteit & Firoozabadi Reference Hoteit and Firoozabadi2009; Tian et al. Reference Tian, Zhang, Lei, Yin, Kazemi and Wu2021).
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:
(i) the fluid is assumed to be a homogeneous distribution in an individual nanopore;
(ii) the high permeability porous medium is considered to have infinite conductivity when compared with the ultralow permeability of the nanoporous medium;
(iii) the nanoporous medium is viewed as a collection of parallel capillaries;
(iv) chemical reactions or adsorption phenomena are ignored;
(v) the system is assumed to be isothermal during codiffusion and counterdiffusion processes;
(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 Reference Kerkhof and Geboers2005b), equivalently, the generalized MS equation (Bhatia et al. Reference Bhatia, Bonilla and Nicholson2011),
\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 Reference Kerkhof, Geboers and Ptasinski2001; Bhatia et al. Reference Bhatia, Bonilla and Nicholson2011):
\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}
\begin{gather}{v_{i,\zeta }}(r = R) ={-} G\frac{{\textrm{d}{v_{i,\zeta }}}}{{\textrm{d}r}},\end{gather}
\begin{gather}\frac{{\partial {v_{i,\zeta }}}}{{\partial r}}(r = 0) = 0,\end{gather}
\begin{gather}{\mu _i}(\zeta = 0) = {\mu _i}(\kern 1.5pt {p_{in}},{x_{i,in,i = 1,2, \ldots ,{n_c}}}),\end{gather}
\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 Reference Kerkhof and Geboers2005b)
\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 nc mass conservation equations with 
${n_c}$ components in both mediums. The governing equations for the nanoporous medium can be expressed as (Deen Reference Deen2013)
\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,
\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 (Reference Kerkhof1996) 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 (Reference Kerkhof and Geboers2005b) 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 (Reference Kerkhof1996); Pant et al. (Reference Pant, Mitra and Secanell2013))
\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 Reference Krishna and Wesselingh1997; Ma & Emami-Meybodi Reference Ma and Emami-Meybodi2024b),
\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. Reference Wu, Li, Wang, Chen and Yu2015; Pang, Fan & Chen Reference Pang, Fan and Chen2021; Liu & Emami-Meybodi Reference Liu and Emami-Meybodi2023). As a result, the modified wall-friction factor can be expressed as follows:
\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 $:
\begin{gather}
{[}{\boldsymbol{N}_\beta
}] ={-} [{\boldsymbol{T}_\beta }][{\nabla _T}\mu
],\end{gather}
\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,
\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}
\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):
\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
\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 Reference Peng and Robinson1976; Michelsen & Mollerup Reference Michelsen and Mollerup2004),
\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):
\begin{gather}p(\zeta ,0) = {p_0},\end{gather}
\begin{gather}{z_i}(\zeta ,0) = {z_{i0}}\quad \textrm{for}\;i = 1,2, \ldots ,{n_c} - 1,\end{gather}
\begin{gather}p(0,t) = {p_{HP}}(t),\end{gather}
\begin{gather}{z_i}(0,t) = {z_{HP,i}}(t)\quad \textrm{for}\;i = 1,2, \ldots ,{n_c} - 1.\end{gather}
\begin{gather}\frac{{\partial p}}{{\partial \zeta }}(L,t) = 0,\end{gather}
\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 (CH4), n-butane (C4H10) and n-dodecane (C12H26) – within a nanopore under steady-state conditions by solving (2.3). For simplicity, we use C1, C4 and C12 to refer to CH4, C4H10 and C12H26, 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−1) 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 0 = 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−1 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 CH4 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−1, 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 Reference Sutera and Skalak1993).
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−1. 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 Reference Patankar1980). 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,
\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 ~ Pem to scale the diffusion contribution for the systems with various parameters (Marle Reference Marle1981). 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:
\begin{gather}S{h_i} \approx \frac{1}{{{\omega _i}}},\quad P{e_i} < 0.1,\end{gather}
\begin{gather}S{h_i} = 0,\quad P{e_i} \approx \frac{1}{{{\omega _i}}},\end{gather}
\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 Pei = 0.1 and Pei = 100. In transition systems with 0.1 < Pei < 100, diffusion and advection play significant roles. Within this Pe range, (3.2b) represents a diffusion–advection neutral point (Shi = 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 CO2 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 CO2 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 CO2 exhibits partial miscibility with the alkane mixture. The PR-EOS (Peng & Robinson Reference Peng and Robinson1976) is employed to calculate the phase behaviour of the alkane and CO2 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. (Reference Falk, Coasne, Pellenq, Ulm and Bocquet2015) and Perez & Devegowda (Reference Perez and Devegowda2020) and the field data reported by Freeman et al. (Reference Freeman, Moridis, Michael and Blasingame2012). 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 CO2 profile within the nanoporous medium after different time counterdiffusion, respectively. The mixing extent of the CO2 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 (Reference Waggoner, Castillo and Lake1992) 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 CO2 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. Reference Lake, Johns, Rossen and Pope2014). As depicted in figure 7(b), in the case of counterdiffusion, CO2 diffuses into the nanoporous medium and mixes with alkanes, creating a ‘mixing region’, which grows as CO2 penetrates deeper into the nanoporous medium. Over time, the CO2 mole fraction becomes very high in the near-fracture region, forming a CO2 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 CO2 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. Reference Behmanesh, Clarkson, Tabatabaie and Heidari Sureshjani2015),
\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 CO2 mole fraction profiles along the nanoporous medium. Three regions (CO2 bank; mixing region; alkane mixture bank) are separated after 120 months of counterdiffusion based on the boundary defined by CO2 mole fractions of 0.1 and 0.9. (d) Mixing region and CO2 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 CO2 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 CO2 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 CO2 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 CO2 mole fraction profile between 60 and 120 months in figure 7(b). That is the reason why most of the CO2 is still confined to the near high-permeability medium even after a decade.
During the counterdiffusion of CO2 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 Reference Onsager1931). 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 (Reference Krishna and van Baten2005) 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 Reference Deen2013; Firoozabadi Reference Firoozabadi2016) 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 Reference Ma and Emami-Meybodi2024a,Reference Ma and Emami-Meybodic) 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. Reference Holt, Park, Wang, Stadermann, Artyukhin, Grigoropoulos, Noy and Bakajin2006), 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. Reference Holt, Park, Wang, Stadermann, Artyukhin, Grigoropoulos, Noy and Bakajin2006; Whitby et al. Reference Whitby, Cagnon, Thanou and Quirke2008; Javadpour Reference Javadpour2009) and dependent on the fluid and pore wall solid characteristics. This parameter is lumped into
\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 CO2 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 CO2 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 CO2 remains confined to the near high-permeability medium due to the slow species counterdiffusion rate driven by the composition gradient.
Acknowledgements
The authors also thank the member companies of the Subsurface Energy Recovery and Storage Joint Industry Partnership (SERS JIP) at The Pennsylvania State University for their financial support. Dr H. Emami-Meybodi holds the Dr C.H. Bowman and L.A. Holleran Early Career Professorship in Petroleum and Natural Gas Engineering and the E. Willard and R.S. Miller Faculty Fellowship at Pennsylvania State University.
Funding
This work is supported by the National Key Research and Development Program of China under grant (2022YFE0206700). Acknowledgment is also made to the National Natural Science Foundation of China (52204057), Science Foundation of China University of Petroleum, Beijing (2462021YJRC012 and 2462021BJRC003), and the donors of the American Chemical Society Petroleum Research Fund (ACS PRF 65892-ND9).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Numerical scheme
A.1. Numerical solution of the velocity profile in nanopores
Assuming that mole fraction xi does not depend on the radius (Young & Todd Reference Young and Todd2005), the momentum balance equation, (2.3), is rearranged in a concise form by combining parameters that are independent of the radius,
\begin{gather}{B_i} + {\eta _i}\frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial {v_i}}}{{\partial r}}} \right) + \sum_{\scriptstyle j = 1 \atop \scriptstyle j {\ne} i }^n {{A_{ij}}({v_j} - {v_i})} = 0,\end{gather}
\begin{gather}{B_i} ={-} {c_i}{\left. {\frac{{\textrm{d}{\mu_i}}}{{\textrm{d}\kern 0.06em x}}} \right|_T},\end{gather}
\begin{gather}{A_{ij}} = {c_t}RT\frac{{{x_i}{x_j}}}{{{D_{ij}}}}.\end{gather}Expanding the rate of momentum change of the viscous term using the chain rule, the momentum balance equation can be expressed as follows:
\begin{equation}{\eta _i}\frac{{{\partial ^2}{v_i}}}{{\partial {r^2}}} + \frac{{{\eta _i}}}{r}\frac{{\partial {v_i}}}{{\partial r}} + \sum_{\scriptstyle j = 1 \atop \scriptstyle j {\ne} i }^n {{A_{ij}}({v_j} - {v_i})} + {B_i} = 0.\end{equation} The above equation is discretized using the finite-difference method. The radius is divided into 
$m - 1$ segments with m grids, and the residual equation at grid m can be expressed as follows:
\begin{align}
{R_{i,m}} &= \left[ {\dfrac{{{\eta_i}}}{{{{(\Delta
r)}^2}}} + \dfrac{{{\eta_i}}}{{2{r_m}\Delta r}}}
\right]{v_{i,m + 1}} - \dfrac{{2{\eta _i}}}{{{{(\Delta
r)}^2}}}{v_{i,m}} + \left[ {\dfrac{{{\eta_i}}}{{{{(\Delta
r)}^2}}} - \dfrac{{{\eta_i}}}{{2{r_m}\Delta r}}} \right]\nonumber\\
&\quad + {v_{i,m - 1}} + \sum_{\scriptstyle j =
1 \atop \scriptstyle j {\ne} i }^n
{{A_{ij}}({v_{j,m}} - {v_{i,m}})} + {B_i}.
\end{align}We can also get the discrete boundary conditions with second-order precision,
\begin{gather}{R_{i,N}}
= {v_{i,N}} - {L_s}\frac{{4{v_{i,N - 1}} - {v_{i,N - 2}} -
3{v_{i,N}}}}{{2\Delta r}} = 0,\end{gather}
\begin{gather}{R_{i,1}}
= \frac{{4{v_{i,2}} - {v_{i,3}} - 3{v_{i,1}}}}{{2\Delta r}}
= 0.\end{gather} Total 
$m \times {n_c}$ velocities need to be calculated for 
${n_c}$ components, which can be solved by combining the (A3) and (A4) with the Newton–Raphson method.
A.2. Numerical solution of pressure and composition in the nanoporous medium
The nanoporous medium is divided into several grids, and the pressure and composition are solved numerically. Central differencing is applied to discretize the governing equations in space, and backward differencing is used for their discretization in time (Abou-Kassem, Islam & Farouq Ali Reference Abou-Kassem, Islam and Farouq Ali2020). The discretized equations are then written in a residual form,
\begin{align}
R_{n,i}^{t + 1} &= {\left( {\dfrac{{{V_b}}}{{\Delta t}}}
\right)_n}[({x_i}{c_l}\phi {S_l} + {y_i}{c_v}\phi
{S_v})_n^{t + 1} - ({x_i}{c_l}\phi {S_l} + {y_i}{c_v}\phi
{S_v})_n^t]\nonumber\\
&\quad - \sum_{m \in \psi } {\left\{
{\sum_\beta {\sum_{j = 1}^{{n_c}} {\left(
{\dfrac{A}{{\Delta \zeta }}{S_\beta }{T_{\beta ,ij}}}
\right)_{nm}^{t + 1}[({\mu_{j,T}})_m^{t + 1} -
({\mu_{j,T}})_n^{t + 1}]} } } \right\}} \nonumber\\
& \quad \quad
\textrm{for}\;i = 1,2, \ldots ,{n_c},\;n = 1,2, \ldots
,{n_G},
\end{align}
where 
$\psi $ is a set of the grid that is connected with grid n, A is the area perpendicular to the 
$\zeta $ direction, 
$\Delta \zeta $ is the distance along the 
$\zeta $ direction, 
${V_b}$ is the bulk volume, superscript t is the time index, t + 1 represents the current time step, t represents the last time step, subscript n and m are the grid index, β is the phase index.
Each grid has nc independent primary variables and nc independent conservation equations. The primary variables of each grid are the pressure (p) and total mole fraction of each component zi (nc − 1 independent value). We use the sum of all nc conservation equations corresponding to the pressure
\begin{equation}R_{n,overall}^{t + 1} = \sum_{i = 1}^{{n_c}} {R_{n,i}^{t + 1}} .\end{equation} The system of governing equation for grid n becomes 
$\boldsymbol{R} = (R_{n,overall}^{t + 1},R_{n,1}^{t + 1}, \ldots ,R_{n,{n_c} - 1}^{t + 1})$ with the exact ordering of primary variables as 
$\boldsymbol{v} = (\,{p_n},{z_{n,1}}, \ldots ,{z_{n,{n_c} - 1}})$.
Figure 8 depicts a comprehensive flow chart outlining the various modules involved in solving the governing equation system. These modules include the input module, control module, secondary variables' calculation module, solver module and output module. The input module reads the grid information and initial porous media properties. Then, the pressures and total mole fractions of fracture and each grid is chosen as the initial primary variables (
${\boldsymbol{v}^{\boldsymbol{I}}}$). The control module controls the main time loop, starting with the user-specified time steps. In the secondary variables’ calculation module, liquid phase and vapour phase mole fraction (
${x_i},{y_i}$), molar density (
${c_l},{c_v}$) and fugacity (
$\,{f_i}$) are calculated based on a volume-liquid-equilibrium calculation using PR-EOS. Additionally, the saturation (
${S_l},{S_v}$), partial viscosity fraction (
${\kappa _{l,i}},{\kappa _{v,i}}$), effective MS diffusivities (
$D_{ij}^e$) and effective Knudsen diffusivities (
$D_{iK}^e$) are updated in the secondary variables calculation module. Subsequently, the residual equations system combining the conservation equation of fracture and total 
${n_G}$ grids is constructed. Then the Newton–Raphson method is used to calculate the primary variables for each time step in the solver module. When the governing equations at the current time step are solved, the next time step begins with the primary variables of the current time step as the initial primary variables. Finally, the data is stored in the output module, and the simulation is ended when the user-specified final time is reached.
Figure 8. The numerical scheme for solving governing equation system (Ma & Emami-Meybodi Reference Ma and Emami-Meybodi2024b).
Appendix B. Transport coefficients
B.1. Bulk fluid properties
In order to ensure that the total shear stress on all components of a uniform fluid is equal to that on the mixture as a whole, the partial viscosity is defined as (Bhatia & Nicholson Reference Bhatia and Nicholson2008)
\begin{equation}{\eta _i} = {\omega _i}\eta ,\end{equation}where ωi is the weight fraction of component i. Additionally, the fractional viscosity coefficient, which is employed for the calculation of the modified wall friction factor, is given by the following equation (Kerkhof et al. Reference Kerkhof, Geboers and Ptasinski2001):
\begin{equation}{\kappa _i} = \frac{{{\eta _i}}}{{{c_i}RT}}.\end{equation}For the bulk binary diffusion coefficient, the Fuller, Schettler & Giddings (Reference Fuller, Schettler and Giddings1966) approach has been widely used to calculate the ideal gas binary diffusion coefficient. However, considering the reservoir pressure and temperature, we use a unified model proposed by Leahy-Dios & Firoozabadi (Reference Leahy-Dios and Firoozabadi2007) to calculate the non-ideal multicomponent diffusion coefficient under bulk conditions.
B.2. Porous-specific parameters
The aforementioned equation only refers to bulk fluid properties (viscosity and diffusion coefficient). Several porous-specific parameters must be known or estimated to predict the fluxes of each phase within the porous medium. These parameters include intrinsic permeability, K, obstruction factor, 
${Q_m}$, and the Knudsen radius, 
${Q_k}$. Experiments and molecular simulation are the best way to get these porous-specific parameters. However, the primary focus of this study is not the experimental determination of these parameters. Instead, we assume a porous medium as a bundle of sinuous yet parallel capillaries, and the following correlations are provided for the estimation of these parameters.
The intrinsic permeability is calculated by the Carman–Kozeny equation (Peters Reference Peters2012),
\begin{equation}K = \frac{{\phi {d^2}}}{{32\tau }},\end{equation}
where d is the average pore diameter, and 
$\tau $ is the tortuosity.
A correction factor is required to account for the obstruction by the soil grains. The bulk binary diffusion coefficient, 
${D_{ij}}$, is replaced by the effective binary diffusion coefficient, 
$D_{ij}^e$, which is given by (Abu-El-Sha'r & Abriola Reference Abu-El-Sha'r and Abriola1997)
\begin{equation}D_{ij}^e = {Q_m}{D_{ij}},\end{equation}
where 
${Q_m}$ is the obstruction factor, which is assumed to be solely a function of the porous medium. For a bundle of capillaries, 
${Q_m}$ can be expressed as (Epstein Reference Epstein1989)
\begin{equation}{Q_m} = \frac{\phi }{\tau }.\end{equation}Similarly, the Knudsen diffusivity is determined by the mean thermal velocity and a factor related to porous medium geometry (Jackson Reference Jackson1977; Gruener & Huber Reference Gruener and Huber2008),
\begin{equation}D_{iK}^e = {Q_k}\sqrt {\frac{{8RT}}{{{\rm \pi} {M_i}}}} ,\end{equation}
where 
${M_i}$ is the molecular weight of component i, and 
${Q_k}$ is the Knudsen radius, which is proportional to the average pore radius and is assumed constant for a given medium as a first degree of approximation. With an assumption of the Maxwell slip boundary condition, the Knudsen radius for a bundle of capillaries can be expressed as (Pant et al. Reference Pant, Mitra and Secanell2013)
\begin{equation}{Q_k} = 0.89\frac{{\phi {d_a}}}{{3\tau }}.\end{equation}B.3. Weighting coefficient
The weighting coefficient for Knudsen diffusion is the ratio of the wall–molecule collision frequency to the total collision frequency (Wu et al. Reference Wu, Li, Wang, Chen and Yu2015),
\begin{gather}{\omega _K} = \frac{1}{{1 + \dfrac{2}{{Kn\left( {1 + \dfrac{1}{\varepsilon }} \right)}}}},\end{gather}
\begin{gather}{\omega _v} = \frac{1}{{1 + \dfrac{{Kn\left( {1 + \dfrac{1}{\varepsilon }} \right)}}{2}}},\end{gather}
where Kn is the Knudsen number, and 
$\varepsilon $ is the aspect ratio for the cross-section shapes. The value of 
$\varepsilon $ is 1 for the capillary tube in this paper (Pang et al. Reference Pang, Fan and Chen2021).
The Knudsen number is a dimensionless number that quantifies the importance of molecular transport processes in comparison with continuum transport processes, which can be expressed as
\begin{equation}Kn = \frac{\lambda }{d},\end{equation}where λ is the mean free path of gas molecules. There are several methods to calculate the mean free path (Liehui et al. Reference Liehui, Baochao, Yulong and Zhaoli2019), and the following equation for dense gas is used in this paper (Corral-Casas et al. Reference Corral-Casas, Li, Borg and Gibelli2022):
\begin{equation}\lambda = \frac{{16\eta }}{{5{\rm \pi} p}}\sqrt {\frac{{{\rm \pi} RT}}{{2M}}} .\end{equation}The gas mixture's viscosity and averaged molecule weight are used in the above equation to calculate the gas mixture's mean free path.
Appendix C. Inner boundary condition
A material balance equation (MBE) for the high-permeability porous medium is solved to obtain the fluid composition within the high-permeability porous medium at a given time. We define the following MBE for the total moles of component i in the high-permeability porous medium:
\begin{gather}{n_{HP,i}}({t^{n + 1}}) = {n_{HP,i}}({t^n}) + [{N_{Inter,i}}({t^{n + 1}}) - {N_{q,i}}({t^{n + 1}})]\Delta t,\end{gather}
\begin{gather}\Delta t = ({t^{n + 1}} - {t^n}),\end{gather}
where 
${n_{HP,i}} = {V_{HP}}{c_{HP}}{z_{HP,i}}$ is the total moles of component i in the high permeability porous medium; 
${c_{HP}}$ is the overall molar concentration of the fluid in the high permeability porous medium; 
${V_{HP}}$ is the high-permeability porous medium pore volume; 
$\Delta t$ is a time step; 
${t^n}$ represent the value at the previous time step; 
${t^{n + 1}}$ represent the value at the current time step; 
${N_{Inter,i}}$ is the molar flux of component i through the interface between the high permeability and nanoporous media, which can be obtained from BFM. 
${N_{qi}}$ is the total molar flux for component i removed from (+ sign) or added to (− sign) the high-permeability. The 
${N_{qi}}$ value can be obtained by solving a volume balance equation for the high-permeability porous medium,
\begin{gather}{q_s}({t^{n
+ 1}})\Delta t = \frac{\displaystyle{\sum_{i = 1}^{{n_c}}
{{n_{HP,i}}({t^{n + 1}})} }}{{{c_{HP}}({t^{n + 1}})}} -
{V_{HP}},\end{gather}
\begin{gather}{N_{qi}} = {\lambda _l}{c_{lq}}{q_s}{x_{qi}} + {\lambda _v}{c_{vq}}{q_s}{y_{qi}},\end{gather}
where 
${q_s}$ is volumetric flux; 
${c_{lq}}$ and 
${c_{vq}}$ are the overall molar concentration of liquid and vapour phase removed from/added to the high permeability porous medium, respectively; 
${x_{qi}}$ and 
${y_{qi}}$ are the molar composition of the liquid and vapour fluid removed from/added to the high permeability porous medium, respectively; 
${\lambda _l}$ and 
${\lambda _v}$ are mobility of liquid phase and vapour phase. Each phase's molar concentration and composition can be calculated based on PR-EOS with the given pressure and fluid overall composition. The mobility can be calculated based on relative permeability,
\begin{gather}{\lambda _l} = ({k_{rl}}/{\eta _l})/({k_{rl}}/{\eta _l} + {k_{rv}}/{\eta _v}),\end{gather}
\begin{gather}{\lambda _v} = ({k_{rv}}/{\eta _v})/({k_{rl}}/{\eta _l} + {k_{rv}}/{\eta _v}).\end{gather}Appendix D. Dimensionless number groups
For the multicomponent system, the overall velocity of component i (
${v_i}$) is dependent on average velocity (
$\tilde{v}$), fluid molar density (
${C_t}$), diffusion coefficient (
${D_{ij}}$), partial viscosity (
${\eta _i}$), pore length (L) and pore radius (
${R_p}$). The dimensional analysis then shows that they should be expressed in the form
\begin{gather}S{h_i} = f(P{e_i},{\varphi _i},x_{ij}^\ast ,\eta _{ij}^\ast ,D_{ij}^\ast )\quad \textrm{for}\;i = 1,2, \ldots ,{n_c},\end{gather}
\begin{gather}S{h_i} = \frac{{{v_i}L}}{\displaystyle{\sum_{j = 1}^{n - 1} {{{D^{\prime}}_{ij}}\Delta {x_j}} }},\end{gather}
\begin{gather}P{e_i} = \frac{{\tilde{v}L}}{\displaystyle{\sum_{j = 1}^{n - 1} {{{D^{\prime}}_{ij}}\Delta {x_j}} }},\end{gather}
\begin{gather}{\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{gather}
\begin{gather}x_{ij}^\ast = \frac{{{x_j}}}{{{x_i}}}\quad \textrm{for}\;i = 1,2, \ldots ,{n_c};\;j \ne i,\end{gather}
\begin{gather}\eta _{ij}^\ast = \frac{{{\eta _j}}}{{{\eta _i}}}\quad \textrm{for}\;i = 1,2, \ldots ,{n_c};\;j \ne i,\end{gather}
\begin{gather}D_{ij}^\ast = \frac{{{D_{ij}}}}{{{D_{in}}}}\quad \textrm{for}\;i = 1,2, \ldots ,{n_c} - 1;\;j \ne i,\end{gather}
where 
$S{h_i}$ is i component Sherwood number, reflecting the ratio of overall to diffusive transport rate; 
$P{e_i}$ is i component Péclet numbers, reflecting the ratio of advective to diffusive transport rate; 
${\varphi _i}$ is fluid–solid friction modulus, reflecting the shear stress's impact on diffusive transport; 
$x_{ij}^\ast $ is mole fraction ratio; 
$\eta _{ij}^\ast $ is partial viscosity ratio and 
$D_{ij}^\ast $ is diffusion coefficient ratio.
The fluid mole fraction ratio, partial viscosity ratio and diffusion coefficient ratio are fluid characteristics, and they interplay with each other. We designate them all inclusively by the fluid modulus (
$\mathcal{F}$), (D1) can then be expressed as
\begin{equation}S{h_i} = f(P{e_i},{\varphi _i},\mathcal{F})\quad \textrm{for}\;i = 1,2, \ldots ,{n_c}.\end{equation}Along with the fluid–solid friction modulus, 
$\mathcal{F}$ also influences the relationship between Sh and Pe. However, predicting its specific impact on Sh and Pe is a complex task due to the intricate interplay among various fluid characteristics parameters. This complexity increases as the number of species increases in multicomponent fluids. When the fluid–solid friction modulus is kept constant, the effect of the fluid modulus on Sh and Pe is restricted to the transition system and minimal in systems where diffusion and advection are dominant. For the purposes of this study, we assume that we already have knowledge of the fluid properties and proceed to explore the relationship between the multicomponent Sherwood number (Sh), Péclet number (Pe) and fluid–solid friction modulus (φ). Therefore, the implicit function between Sh and Pe for any component in multicomponent systems can be expressed as (3.3).
