Fluid transport through heterogeneous pore matrices: Multiscale simulation approaches in a with slip

Fluids confined in nanopores exhibit several unique structural and dynamical characteristics that affect a number of applications in industry as well as natural phenomena. Understanding and predicting the complex fluid behavior under nano-confinement is therefore of key importance, and both experimental and computational approaches have been employed toward this goal. It is now feasible to employ both simulations and theoretical methods, the results of which can be validated by cutting-edge experimental quantification. Nevertheless, predicting fluid transport through heterogeneous pore networks at a scale large enough to be relevant for practical applications remains elusive because one should account for a variety of fluid–rock interactions, a wide range of confined fluid states, as well as pore-edge effects and the existence of preferential pathways, which, together with many other phenomena, affect the results. The aim of this Review is to overview the significance of molecular phenomena on fluid transport in nanoporous media, the capability and shortcomings of both molecular and continuum fluid modeling approaches, and recent progress in multiscale modeling of fluid transport. In our interpretation, a multiscale approach couples a molecular picture for fluid interactions with solid surfaces at the single nanopore level with hierarchical transport analysis through realistic heterogeneous pore networks to balance physical accuracy with computational expense. When possible, comparison against experiments is provided as a guiding roadmap for selecting the appropriate computational methods. The appropriateness of an approach is certainly related to the final application of interest, as different sectors will require different levels of precision in the predictions.


INTRODUCTION
Fluid transport through heterogeneous pore matrices over broadly different time and length scales is at the heart of many technological and environmental processes. [1][2][3][4][5][6][7][8][9][10] Achieving a molecularlevel description of fluid-solid interactions and the correlated interfacial dynamics as well as fluid diffusion through heterogeneous pore matrices would help to build quantitative functionalitymorphology-transport relationships to enable the manufacture of more efficient and selective porous materials. [11][12][13][14][15][16] Molecular simulations can provide molecular structure and dynamics within the interfacial layers, where complex interactions between the porous matrices and fluids occur. [17][18][19][20][21][22][23][24][25] Such simulations, however, tend to be limited to single nano-to meso-pores of simple geometry, which falls short of accounting for the porous material complexity at larger length scales. On the other hand, approaches accounting for the latter must rely on a simplified picture of the confined fluid and its transport properties. [26][27][28] Achieving a balance between the two approaches remains a challenge, especially because different applications require different levels of precision in the predictions, while in most cases, it is preferable to obtain results at low to moderate computational costs.
In what follows, we present the selected simulation results to highlight the significant impact of the fluid structure on its dynamical behavior in a single pore at the nanoscale and compare them against experimental data at various length scales to assess the reliability of various approaches, when possible. Convective motion and molecular diffusion are the two main mechanisms responsible for fluid dynamics in porous matrices. Molecular diffusion becomes the dominant mechanism when fluids are confined in pores of size comparable to the size of the confined fluid molecules (i.e., in nanopores). The geometry of the nano-/meso-pores and their pore apertures as well as pore connectivity play an important role in regulating the amount of fluids that cross or enter natural porous

REVIEW
scitation.org/journal/phf media. [29][30][31] Because many pores in the subsurface are in this size regime, 8 it is crucial to understand the mechanisms by which the molecular diffusion occurs in nanopores, and when surfaces and pore defects strongly affect the transport mechanisms. To illustrate these concepts, we summarize the selected studies for fluid transport in slit-shaped single nanopores primarily from our body of work, produced by employing atomistic molecular dynamics (MD) simulations. It should be recognized that many groups contributed to these investigations. The interested reader is referred to a few recent reviews on this topic. 18,[32][33][34] Bridging the gaps between the results obtained within a single pore, a few pores, a few hundreds of pores, and ultimately within complex pore networks requires the development of simulation approaches with high computational efficiency. In principle, atomistic MD is able to describe molecular phenomena even when they occur in large systems. However, it requires high computational resources, and therefore, it can become inefficient and sometimes impractical to study fluid transport through heterogeneous porous matrices at length and time scales larger than a few nm and several hundreds of ns, respectively. 35 When practical applications require modeling such conditions, continuum fluid models, such as the no-slip Hagen-Poiseuille flow equation or Darcy's law, are often employed, although they do not correctly characterize flow within complex nanoporous materials because in such systems, non-continuum flow phenomena arise due to the well-ordered molecular structure near solid-liquid interfaces, inaccurate characterization of the local viscosity, and interfacial slip. 35 In Fig. 1, we   FIG. 1. Time and length scales in which atomistic MD simulations, lattice Boltzmann, kinetic Monte Carlo, macroscopic continuum fluid model, and multiple modeling approaches are used. The approach identified as "multiscale" seeks to reconcile the results from the different methods toward describing fluid transport both precisely (e.g., taking into consideration molecular phenomena) and effectively (e.g., achieving length scales relevant for the applications). Reproduced with permission from Vlachos, "A review of multiscale Analysis: Examples from systems biology, materials engineering, and other fluid-surface interacting systems," Adv. Chem. Eng. 30, 1-61 (2004). 36 Copyright 2005 Elsevier B.V. summarize time and length scales in which atomistic MD simulations as well as Monte Carlo, lattice Boltzmann (LB), Kinetic Monte Carlo (KMC), macroscopic continuum fluid model, and multiple modeling approaches are typically used. Developing a computational framework to couple all these techniques remains challenging but promises the ability to capture phenomena that occur at the single-pore level (e.g., edge effects) with larger-scale phenomena (e.g., preferential flow pathways) toward truly predicting fluid transport through heterogeneous porous matrices. Note that in the shortest time and length scales, one could implement quantum mechanics methods to elucidate interactions between fluids and pore surfaces, which could lead to reactivity.
We review in what follows some promising advances for the development and implementation of multi-scale approaches achieved by implementing non-equilibrium MD simulations, coupling MD simulations with LB calculations, stochastic approaches based on KMC, and improved continuum macroscopic flow models. We attempt to identify conditions under which each multiscale approach is preferable, based on comparison against the available experimental studies. The experimental tools considered include pulsed field gradient nuclear magnetic resonance (NMR), 37-39 quasielastic neutron scattering (QENS), [38][39][40] and microimaging, [41][42][43] among others.
The remainder of this Review is organized as follows: In the section titled Fluid transport in single nanopores, we present the selected MD results for the fluid structure and transport behavior in single nanopores and discuss the inaccuracy of the conventional continuum fluid model owing to the dominance of non-continuum flow phenomena. Then, in the section titled Fluid transport in hierarchical pores, we present some multiscale approaches that address the aforementioned challenges. We conclude presenting open questions and some emerging topics we consider promising testbeds for future research.

FLUID TRANSPORT IN SINGLE NANOPORES
Fluid structure governs diffusion enhancement/hindrance Mounting evidence, albeit mostly from computational studies, suggests that confinement significantly affects almost all the fluid properties related to interfacial interactions (adsorption and solubility), thermophysical (phase behavior), dynamical (translational, rotational, and vibrational diffusions), and ultimately transport properties of the confined species. 8,47 Because understanding the structural and dynamical properties of fluid species such as carbon dioxide (CO 2 ) and methane (CH 4 ) confined in nano-and mesopores in sedimentary rocks is of increasing interest, MD simulations at reservoir composition 48,49 have been used by Loganathan et al. 44 to quantify the partitioning of CO 2 and methane between bulk and nano-and meso-pores bounded by the montmorillonite clay mineral. MD has been proven useful to study multi-component systems described at the atomistic level. In MD, interactions between atoms and molecules are described explicitly via appropriate atomic force fields. The motion of such entities is described by solving numerically Newton's equations of motion. The resultant molecular trajectories, integrated over time, yield the macroscopic properties of such systems via statistical mechanics analysis. 50 For example, in Fig. 2 (left), the results show the preferential adsorption of CO 2 on clay surfaces, suggesting that CO 2 can effectively remove methane from such pores, leading to enhanced oil recovery and shale gas production. It is worth pointing out that the accumulation of fluids near pore surfaces leads to more frequent fluid-fluid and fluid-pore collisions, which reduce the mean free path (MFP) and the diffusion coefficient, and that this effect is more and more pronounced as the pores become narrower and narrower, e.g., micro-pores. By employing the Direct Simulation Monte Carlo (DSMC) method and MD simulations, Xie and his co-workers 51 showed that the MFP is spatially dependent when both the intermolecular interactions between gas molecules and collisions between gas molecules and wall atoms are taken into consideration. The behavior of fluid confined within larger meso-pores, on the other hand, may exhibit fluid-substrate interfacial and bulk-like properties simultaneously, with the latter becoming dominant as the pore width increases. 45 To support simulation results such as those in Fig. 2, one could employ QENS. For example, conducting QENS for ethane in mesoporous silica pores, Patankar et al. found that the self-diffusion coefficient of ethane confined in 11.1 nm and 41.5 nm pores was ∼4 times slower than that observed in the bulk under otherwise similar conditions. 45 If these results are due to the enhanced fluid density near the pore walls, it might be feasible to facilitate the fluid diffusion by adding a second fluid that preferentially adsorbs on the solid substrate, as discussed previously. Indeed, one QENS experimental study 45 reported that the ethane diffusion in mesoporous silica increased in the presence of CO 2 (27 mol. % CO 2 ) by a factor of ∼2 compared to the data obtained for pure ethane [see Fig. 2(a), right]. Le et al. 23 conducted MD simulations for systems containing butane and CO 2 in a ∼2 nm wide slit-shaped silica pore. The results also showed that CO 2 preferentially absorbed near the pore surfaces, consistent with the results reported in Fig. 2(a) (left), 44 enhanced diffusion of butane occupied in the middle of the pore, via a "molecular lubrication" mechanism, which is in good agreement with the prior experimental study. 23 Others also observed that CO 2 enhanced the diffusion of hydrocarbons through different porous materials. [52][53][54] The preferential adsorption of CO 2 to the pore surfaces can enhance hydrocarbon diffusion in nanopores, but one can question whether it is possible to achieve similar results if another fluid other than CO 2 , also preferentially attracted or even more strongly

REVIEW
scitation.org/journal/phf adsorbed to the pore surface, is used. Recently, Loganathan et al. 46 employed Grand Canonical MD (GCMD) simulations to investigate the structure and dynamics of CO 2 -water mixtures in the interlayers of the smectite clay, Na-hectorite, at geological temperatures and pressures. The GCMD algorithm is a hybrid MD-Monte Carlo technique, which was introduced to maintain constant temperature and chemical potential during the virtual experiment. The approach assures that the number of particles and system energy varies corresponding to the standard grand canonical probability distribution. 55 In the GCMD approach, it is possible to use the grand canonical Monte Carlo algorithm to first establish and subsequently maintain a chemical potential gradient across a pore, while MD allows one to study the transport of fluids across the pore in response to the chemical potential gradient. Details of the GCMD simulations can be found in the work of Boinepalli and Attard. 55 In Fig. 2 showed that water is more strongly adsorbed on the clay surfaces than CO 2 , displacing CO 2 to the midplane of clay interlayers for the monolayer hydration structure (left) and closer to one of the basal surfaces for the bilayer hydration structure (right). Coming back to the question whether diffusion enhancement is observed if another fluid other than CO 2 , e.g., water, is mixed with hydrocarbon, Gautam et al. 56 used QENS experiments to examine the diffusion of propane in silica-based, cylindrical pores of diameter ∼1.5 nm in the presence of water. The experimental results showed that water hampered, rather than facilitated, propane diffusion. To provide a molecular understanding of these observations, Le et al. 25 conducted MD simulations for a system resembling the experimental one. 56 The results showed that although water molecules are strongly adsorbed on the pore surface due to such a narrow pore diameter as well as to the strong waterwater interactions, water molecules can form molecular bridges across the pore volume, which impedes propane diffusion. Le et al. 25 quantitatively examined the decrease in propane diffusion upon the increase in water loading inside the pore [as shown in Fig. 2(b), right]. The approaches implemented to conduct these investigations could be relevant for studying systems with extremely heterogeneous pores of size comparable to that of the fluid molecules. When the two or more fluids confined in a narrow pore are mixable, however, the phenomena of interest could be different compared to those discussed so far, and the diffusion of a fluid dissolved within another fluid, the latter confined in a pore, needs to be investigated.

Confined fluids regulate diffusion
The transport behavior of guest molecules inside narrow pores in the presence of confined fluids is strongly dependent on the chemical properties of mineral solid substrates, as well as on the chemical nature of the fluid molecules in question. Bui et al. 19 reported that the self-diffusion coefficients estimated for methane confined in water-filled pores carved out of common minerals can vary by a factor of ∼4, although the simulations were conducted under the same conditions (i.e., of the same pore width and temperature). Very few experimental studies on the properties of methane in confined water within nanopores have been conducted. However, Hu et al. 57 recently employed NMR to study the behavior of water and methane within ZSM-22, MCM-41, and SBA-15 nanopores. These materials contain mostly the two most common elements on the Earth (O and Si) and possess uniform one-dimensional pore sizes of ∼0.5 nm, 3 nm, and 6 nm, respectively, which enables the quantification of pore size effects on transport. Hu et al. found that the self-diffusion coefficients of methane in water-filled ZSM-22 pores with pore sizes of ∼0.5 nm have a similar order of magnitude to those found in our previous simulation study 21 (∼10 −10 m 2 /s). Additional insights into multicomponent interactions under nanoconfinement were recently achieved by simulations. 58-60 For example, Ho et al. 58 showed that a thin supercritical CO 2 layer formed at water-solid interfaces could help accelerate water flow through a rough hydrophilic nanochannel, acting as a molecular lubricant that converts a stick-to-slip flow in nanopores.
While equilibrium MD simulations can be used to quantify the structure of confined fluids and their diffusion properties, nonequilibrium MD (NEMD) can be implemented to assess the transport properties of fluids confined in narrow pores, although usually in single pores because of the computational requirements. [61][62][63][64] In NEMD algorithms, the equations of motion are integrated numerically as is the case of equilibrium MD, but an external field is applied to impose, for example, a pressure gradient across a pore. More details have been described in some previous works. 61,65,66 Considering as an example the recent results from our group, Phan et al. 21 studied the transport properties of a gas mixture [ethane, methane, and hydrogen sulfide (H 2 S)] through hydrated nanopores. The results showed that H 2 S could permeate the hydrated pores much faster than the other species and that ethane is too large to transverse effectively the water-filled nanopores. This could have significant implications in several applications, such as the design of gas separation membranes and the shale gas sustainable deployment.
One might ask whether a similar transport behavior in narrow pores would be achieved when a fluid other than water, also strongly adsorbed on the pore surface, was present. One such fluid could be benzene, which could be a first approximation for aromatic hydrocarbons trapped in sedimentary rocks. Employing the boundarydriven NEMD (BD-NEMD) algorithm, 61 Phan et al. 67 studied the transport properties of two fluid mixtures (CO 2 -methane and H 2 Smethane) through amorphous silica nanopores filled with benzene. The BD-NEMD simulations are implemented by applying a constant force along the x axis, acting in the direction of the arrows (see Fig. 3, top) to all CO 2 , H 2 S, and methane molecules located in a thin slab (shaded region) of width dext = 20 Å within the permeate region of the simulation box. This external field establishes and maintains a constant pressure difference across the pore network and hence leads to a steady molar flux across the hierarchical porous media. For mixtures, the transport diffusivity of each component is estimated by dividing its permeability, estimated based on the linear relationship between the permeated number of molecules and time, by its solubility within the pore. The system considered in Fig. 3 was built as an approximate representation of organic-rich shale caprocks that contain a significant amount of organic carbon (>11.7 wt. %). Snapshots representing the simulated systems are shown in Fig. 3 (top). 68 Analysis of the simulation results shows that both CO 2 and H 2 S are favorably adsorbed inside the organic-filled pore, in part, displacing benzene. Surprisingly, CO 2 /H 2 S adsorption facilitates methane transport, as illustrated in Fig. 3 (bottom), because both fluids play a role as mobile carriers and possibly trigger benzene swelling, creating favored traveling paths within the pore networks. 67 Results such as those just summarized could help to determine unanticipated transport mechanisms and building up engineering approaches for CO 2 capture and storage in caprocks, albeit their effect on macroscopic transport needs to be quantified. Related to shale gas, it remains of great interest to quantify how hydrocarbons can be produced from the source rocks. Indeed, many investigated how hydrocarbons escape from kerogen. [2][3][4]8,18 Breakdown of continuum fluid models at the nanoscale Falk et al. employed MD to calculate hydrocarbon transport in kerogen [see Fig. 4(a)] and statistical mechanics to analyze the results. 69 They reported that the continuum Darcy's law, derived from the Navier-Stokes equation, 70 a fundamental principle for computational fluid dynamics (CFD) studies, fails to predict hydrocarbon permeability in shale nanoporous matrices-i.e., the kerogen material. In particular, the permeability k (k = K × η with permeance K and viscosity η), an elemental material characteristic, was observed to depend not only on the adsorbed amount but also on the fluid type [see Fig. 4(b)]. 69 The failure of the continuum flow description is likely due to adsorption. Because of the drastic confinement in the nanopores present in kerogen, the physical state of confined alkanes is remarkably different from their bulk counterpart at the same temperature and pressure, changing from the gaseous phase in the bulk to a condensed liquid-like phase under confinement. These results demonstrate that calculating flow properties in the nanopores using bulk viscosity can lead to inappropriate conclusions. Falk et al. 69 showed that also using the viscosities correspondent to the confined alkane density would not yield permeability predictions, using Darcy's law, in agreement with the MD simulation results [see the inset of Fig. 4

Physics of Fluids
To alleviate the shortcomings of Darcy's approach, many corrections have been suggested, for example, by invoking slippage in gas flow via the Klinkenberg effect. 69 Some have attempted to use MD data to correct the Hagen-Poiseuille (H-P) flow equation, introducing a Navier slip boundary condition (BC) for simulating flows in nanotubes with size of a few nm. 71 For example, Walther et al. 72 modified the H-P equation with the pressure correction proposed by Weissberg 73 to take into consideration the membrane-end losses. Nevertheless, these empirical corrections cannot completely account for the entirety of the complex behavior of hydrocarbons in heterogeneous nano-porous materials.
Atomistic MD is probably the most appropriate approach for accounting for non-continuum flow inside nano-porous materials. [2][3][4]8,18 However, this technique requires enormous computational resources when it is desired to study flow through pores longer and wider than a few tens of nanometers. 74 On the other hand, noncontinuum effects ranging from molecular ordering to velocity slip can inhibit the reliability of continuum fluid models such as CFD. 74 Is it possible to achieve more reliable and computationally efficient predictions toward gas transport in heterogeneous pore networks at a larger scale? The challenge lies in upscaling the simulation results achieved in a single pore, toward predicting flow transport through a pore matrix, as different phenomena involved in such multiscale porous matrices can affect the results.

FLUID TRANSPORT IN HIERARCHICAL PORES
Brute force non-equilibrium MD simulations with predictive empirical relations Fluid transport through heterogeneous pore networks is usually characterized through the tortuosity parameter τ (resistivity to transport). 75,76 This parameter relates the flux J to porosity ϕ and flux in the absence of medium J0 through the following expression: J = ϕ/τJ 0 . While porosity ϕ is estimated using adsorption experiments, permeability/transport and NMR experiments can be used to measure both J and J0 (and, therefore, extract τ). 18 From a numerical perspective, tortuosity can be assessed from x-ray tomography experiments 77 using, for example, random-walk (Brownian motion) simulations where the tortuosity is estimated as the square ratio between the walking path length and the length along the longitudinal coordinate.
To replicate such experiments via a direct atomistic simulation approach, Phan et al. 17 employed the BD-NEMD algorithm 61 (as mentioned above) to examine methane flow through hierarchical amorphous silica porous materials. Note that for single component systems, e.g., only methane through porous media, transport diffusivity can be estimated from the pressure gradient and molar flux, which were quantified by fitting the steady-state simulation results to Fick's first law. The pore networks considered in Fig. 5 were meant to represent the simplified structure for some of those often found in shale formations. By conducting simulations on "synthetic" model porous media, the transport properties (e.g., permeability and diffusivity) are collected and correlated with pore structure information. In Fig. 5(a), the results show a strong dependence of the methane permeability through the pores on the frameworks. Indeed, a linear correlation was observed between methane permeability and a characteristic parameter that accounts for porosity, constriction factor, and tortuosity [see Fig. 5(b)]. The constriction factor, a geometric parameter that can become important when the fluid molecules have size comparable to the pore size, was calculated as a function of cross-sectional area A(x) of a porous medium along the direction of transport and its length L. The simulation results suggest that the correlation identified by Phan et al. applies for hierarchical porous materials consisting of both micropores and mesopores. Ultimately, Phan et al. found that the transport properties of methane through porous media scale as a power-law function of porosity and constriction factor [as shown in Fig. 5(c), left], which are descriptors that can be determined experimentally. Phan et al. identified as exceptions those porous systems that containing extremely strong blockages [system 3 (green) or 8 (orange)]. These results suggest that when characterization data are accessible, it is possible to forecast transport properties for engineering and natural materials. For example, a power law is also obtained for Fontainebleau sandstone [see Fig. 5(c), right]. 78,79 The good agreement proposes that the method employed is capable of quantifying molecular effects on fluid transport as well as reliably predicting fluid transport properties through shale rocks using macroscopic pore characterization data as a sole input.
While important new insights can be achieved from the above method, modeling and upscaling approaches for describing fluid transport within realistic and more complex porous materials are still missing, as NEMD cannot be implemented for systems much larger than those considered in Fig. 5. To overcome the computational barriers, numerous computational and theoretical approaches have been attempted, including coarse-grained MD simulations, 80,81 LB calculations, 82-84 stochastic KMC algorithms, 85-87 random-walk particle-tracking (RWPT) method, and modified macroscopic continuum flow model. Each method displays advantages and suffers from limitations, relying upon the simulation time scale and the size of the samples considered. The goal remains to access further

Coupling MD and LB simulations
As already mentioned, in principle, MD simulations could describe microscopic transport phenomena accurately, but their applicability to large length scales (e.g., >100 nm) is drastically restricted by the availability of computing resources. 88,89 As an alternative, since the 1980s, the LB method has become attractive to overcome the computational limitations that MD simulations entail. 90 Compared to MD simulations, the virtue of the LB method is its capability of mesoscale characterization of fluid transport and fluid-solid interactions in porous media with complex boundary geometries, which renders it a computationally affordable and reasonably accurate alternative to model microscopic transport. [91][92][93] The LB method considers fluid flow as a collection of hypothetical

REVIEW
scitation.org/journal/phf particles described by appropriate particle velocity distribution functions. These particles, located on the lattice fluid nodes, exchange and propagate velocity distribution values after their collision with other particles, as well as with the walls, which are represented as solid nodes. 94 A typical LB implementation for microflow requires the input of initial fluid density and momentum and the prescription of BCs (e.g., fluid-solid interface and periodic boundaries), which can be obtained from the correspondent MD models. 84 During the past few decades, the LB approach has been extensively developed and applied to study flow in confinement, [95][96][97][98][99] from transport through T-/Y-junctions, 100 near nanotubes, 101 corrugated surfaces, 102,103 nano-/micro-channels, 95,104-106 and nano-/micro-porous complex media. [107][108][109] Of interest for the present overview is (1) how accurate the LB results can be at the microscale, e.g., compared to their counterpart MD models, and (2) how MD data can be mapped into LB implementations to improve the LB predictions. In fact, MD simulation results are often used to benchmark and identify the applicability range of LB models: 95,[110][111][112][113] Myriad studies have shown that LB algorithms are capable of reproducing MD results, except at the nano-scale interfaces 114 or under conditions of extreme confinement. 115 These findings provide best practice approaches to couple MD onto LB to yield reliable predictions.
Numerous coupling approaches have been developed to further improve the LB accuracy. 96,116,117 One such approach is the two-way coupling of data (e.g., velocity) exchange at the overlapping regions of the MD and LB domains during the simulation, which applies to fluid-substrate interactions for confined single-phase transport as well as to fluid-fluid interactions for multiphase transport. 96,101 This coupling usually applies to the cases where nano-scale simulations are needed regionally only, such as for the description of fluid transport around a nanotube 101 (where LB and MD are applied to the continuum domain and around the nanotube, respectively) or that of fluid transport through a nanochannel 118 (where the bulk region inside the nanochannel is described by the LB method, while the interfacial phenomena, such as fluid slip, are described by MD simulations).
The key to the success of this approach is the agreement between LB and MD results in the overlapping domain or at the interfacial BC between the two subdomains. An example dataexchange protocol is the one in which the molecule velocity at the interface in the MD model is reset to match a given distribution (e.g., Maxwellian) with mean and variance dictated by the LB results; in contrast, the particle velocity distribution in the LB model is reconstructed to match the homogenized MD results near the interface. 96 While this approach is very promising, its computational efficiency is limited by the need of securing convergence between MD and LB methods in the interfacial region.
An alternative coupling approach, which is less computationally demanding because it does not require synchronization between MD and LB datasets, is the "one-way mapping" of the final MD results into the implementations of larger-scale LB models. Successful examples involve the mapping of slip length, 119,120 fluid viscosity, 119 interfacial tension, 121,122 adsorption layers, 123 and velocity profiles 84 as obtained from MD simulations into the correspondent LB models. An exemplar of the latter approach is the use of MD velocity data to identify suitable BCs to be implemented in the correspondent LB models. 84,124,125 In this approach, velocity data from MD simulations are often normalized by the bulk velocity. 107,126,127 Once the normalized LB velocity results match the correspondent MD data, the LB model is regarded as valid and representative of the physics at the solid-gas interface. However, recent results 84,128 suggest that this approach may not be always accurate when using MD or DSMC data. For example, one study 128 showed that two LB models, which implemented (1) the combined scheme of diffuse reflection and bouncedback BCs and (2) discrete Maxwellian BCs, respectively, successfully predicted the normalized velocity from the DSMC; however, only the former scheme reproduced the dimensionless flux, while the latter overestimated the dimensionless flux in transition flow regimes.
Similar observations were achieved in a comparative LB vs MD study from our group. 84 As shown in Fig. 6(a), four LB models implementing different BCs (denoted as BC1-4, details of BCs are referred to as in the caption of Fig. 6) predict similarly agreeable velocity profiles when compared to MD data; however, not all LB models accurately predict the gas permeability correction factor as observed from core-flooding experiments [see Fig. 6(e): overestimation by implementing BC1 and BC3], which implies that the normalized velocity in Fig. 6(a) may not be the proper method to identify correct BCs for the LB calculations. 84 It was instead suggested that using the non-normalized axial velocity 113 data from the correspondent MD models [ Fig. 6(b)] can identify the correct BCs [i.e., see the good agreement of velocity predictions from the MD and LB results by implementing BC4 in Figs. 6(b)-6(d)]. It was shown that by implementing the newly identified BCs in the LB calculations, it is possible to reproduce the permeability correction factor measured experimentally [see Fig. 6(e)].

Slip boundary conditions in LB simulations
As indicated in Fig. 6, the choice of the BC strongly impacts the LB predictions for properties such as permeability. [135][136][137] When using MD data to benchmark LB models, it is, therefore, crucial to understand the applicability of different slip boundaries in the LB method. MD has been used to identify slip BCs in a variety of conditions, sometimes achieving the unexpected results. 138 In LB calculations, the shape of the boundary segments embedded in each lattice site is considered either planar or curved (nonplanar), reflecting the physical structure of porous media. The physical boundary is often approximated as zigzag segments, i.e., the planar boundary. The drawback of this approximation is that it decreases the resolution of images while depicting porous media, and therefore, it can impact the accuracy of LB results; on the other hand, if the image is converted into high-resolution pixeled planar boundaries, the computational time of the subsequent LB simulations can be long, as well as the efforts required to build the model for the LB calculations themselves. The choice of specific BCs for describing gas slip over planar surfaces has been extensively reviewed in the literature, 94,99,130,139 including diffuse reflection, 140,141 specular reflection, 142 the combined bounce-back and specular reflection, 143 and the combined diffuse reflection and specular reflection, 144 In contrast to planar boundaries, the curved boundaries can better maintain the physical structure of the porous media and can overcome the resolution limitations imposed by the planar approximation. However, accounting for the slip over curved, and in general, non-planar boundaries has been a challenging task at the LB level. 145,146 The common solution is to interpolate the local velocity distribution functions according to the relative distance between fluid node, solid node, and wall. 125,[147][148][149] The relationships between the incoming and outgoing velocity distribution functions are then stipulated by the specific slip BCs, e.g., the interpolated bounce-back, 147 interpolated combined bounce-back and diffuse reflection, 125 interpolated combined bounce-back and Maxwellian diffuse reflection, 145,148 and multi-reflection scheme, 149 which are adapted from planar BCs. Alternative solutions to describe slip on curved boundaries may resort to changing the coordinate system, e.g., to adopt a body fitted curvilinear coordinate system when using the combined bounce-back and specular reflection scheme, 150 or to reconstruct velocity BCs, e.g., using Navier's slip length to derive the velocity at the wall, 151 or adopting a method that counterextrapolates slip velocity and non-equilibrium velocity distribution function. 152 Using LB, augmented when needed by MD results, it is now possible to predict macroscopic transport properties (e.g., permeability) for complex pore networks. The approach has been able to predict quantities that can be measured experimentally and also to predict the hydrocarbon productivity of a given formation or the amount of CO 2 that could be sequestered in a reservoir. When rock formations are particularly heterogeneous, as is the case for many shale formations, it might, however, be impractical to explicitly simulate extremely large rock samples, as for LB approaches to succeed, it is necessary to explicitly simulate fluid transport through pore networks large enough to represent the elementary pore volume.

Coupling deterministic MD and stochastic KMC simulations
When it is desired to estimate the transport properties of heterogeneous materials, it might be appropriate to implement stochastic methods and estimate the properties of interest for many representations of the pore media, provided that they share some wellidentified observables. For example, when porosity, pore size distributions, and mineralogy of a formation are known, one could generate several pore networks with such given characteristics, estimate the transport properties in each representation, and then use appropriate averages for predicting the over-all property of the formation. For such an approach to be feasible, a computationally efficient approach must be available. It is within this scenario that stochastic KMC methods become potentially useful. KMC simulations are often used to stochastically describe transport phenomena controlling various engineering and natural processes. [153][154][155][156][157] They can access lengthy time scales (from ms to h) and wide spatial dimensions (from nm to μm) at a relatively low computational cost. 158,159 Fundamentally, KMC takes into consideration state-tostate transition rates to compute trajectories for a system that "wanders" stochastically along the phase space. 159 The KMC implementation requires rate constants that seize the probability per unit Published under license by AIP Publishing

REVIEW scitation.org/journal/phf
time of state-to-state transitions, which can be obtained from other methods, including MD and LB simulations, and experiments. The order of such transitions creates a sample path or trajectory whose statistics follows the so-called Master equation that determines the dynamics of the system. 159 While the approach has been tested successfully in catalysis, 154,[160][161][162] its implementation to study transport, including models and algorithms, has been described at length in prior contributions. 85 KMC allows huge computational savings, compared to MD, as it coarse-grains the molecular trajectories. The energy barrier the system has to overcome in order to get from one energy basin to another determines the waiting time before one event occurs. 159 In a multiscale approach such as represented in Fig. 1, the information of energy barrier could be obtained via atomistic MD studies. For example, to study fluid transport through a porous medium, Apostolopoulou et al. 85 developed a lattice-based 1D KMC model, which was able to reproduce the atomistic MD simulations conducted by Phan et al. 21 to quantify methane transport through hydrated slit-shaped pores carved out of different solid substrates. Quantitative agreement was achieved between KMC and MD results. Because the KMC model is much more computationally efficient than MD, once the necessary input parameters obtained from MD were appropriately used as a KMC input, it was possible to quantify the contribution of various pore network characteristics to methane transport, including a systematic study on pore length, pore connectivity, and pore width. To quantify the effect of pore width, Apostolopoulou et al. 87 extended the KMC model to 3D and predicted fluid diffusivity in mesoscale pores carved out of minerals commonly found in shale rocks. The analytical solution of the diffusion equation as well as against a set of atomistic MD results is used to validate the 3D KMC model. Subsequently, the stochastic 3D KMC model was used to estimate the diffusion of pure methane in pores of varying width. Albeit some, rather small deviations were observed between MD and KMC results, the latter method allowed for a systematic quantification of the effect of pore width on the diffusion coefficient of methane, showing that the mineralogy of slit-shaped pores can affect the transport properties of supercritical methane only when the pores are smaller than ∼3 nm. When the slit pores were wider than ∼5 nm, the bulk-like diffusion coefficients are predicted for methane diffusing through all the materials considered by Apostolopoulou et al. 87 The lattice KMC model developed by Apostolopoulou et al. 85 can be considered as a bottom-up approach for multi-scale studies. Any fabricated or natural networks can be examined, as long as kinetic (diffusion constants) and thermodynamic (interfacial barriers) properties are available. The latter can, in general, be obtained accurately from MD simulations. It is likely that KMC approaches such as the one just described will provide better insights into the diffusion of multiple components in shale formations and possibly aid the formulation of approaches to enhance the natural gas or oil recovery through shale reservoirs, as well as the permeation of CO 2 in geological formations for permanent sequestration.

REVIEW scitation.org/journal/phf
To demonstrate the ability of KMC to predict macroscopic quantities, Apostolopoulou et al. 86 used a 2D version of the KMC algorithm to evaluate the permeability of heterogeneous pore networks meant to replicate the properties of shale samples. They compared the KMC approach to the results from the effective medium theory (EMT) 164 and simplified renormalization methods, 163 which are deterministic approaches often used in practical settings because, in part, of their simplicity. For example, Bhatia and co-workers 165,166 recently used the EMT to predict effective diffusivities in hierarchical mesoporous zeolite materials. The EMT induces the famous Bruggeman's combination rule, 167 which combines the self-diffusivities for the different domains within the materials in a non-linear fashion. Although the EMT is broadly used, it encounters a few drawbacks, in particular, when many pore blockages exist, and the percolation threshold is approached. 18 Apostolopoulou et al. 86 showed that KMC can accurately account for transitions from low-to highpermeability domains within a material. The most valued characteristic of KMC, compared to the two deterministic methods considered, is its ability to capture anisotropy. 86 It is therefore expected that KMC could be applicable to low-connectivity networks and could evaluate the impact of small-scale heterogeneities in which low connectivity exists locally.
When the KMC was used to estimate the permeability of an Eagle Ford shale sample for which data on pore size distributions and one SEM image were available 163 [see Fig. 7(a)], a reasonable agreement with experimental permeability was achieved using ten stochastic realizations for the inorganic, organic, and dualpermeability networks considered representative of the shale rock [see Fig. 7(b)]. [168][169][170] While, as shown in Fig. 7, the KMC approach yields results that are broadly consistent with those from EMT and renormalization approaches, it is likely that more pronounced differences will emerge between the three methods when samples more heterogeneous than the one used by Apostolopoulou et al. are considered. It is worth noting that the precision of the KMC approach can be enhanced by extending the analysis to more 2D SEM images of a plug sample so that local heterogeneities, anisotropy, and porosity are properly accounted for. It will be challenging to efficiently employ the KMC approach to analyze the permeability of an entire plug sample while taking into account realistic molecular phenomena such as adsorption. The extension of the KMC model to 3D will

REVIEW
scitation.org/journal/phf allow for the systematic quantification of the effect of a number of features on the predicted permeability. 171 It should be noted, however that it is difficult to obtain some important information on complex pore matrices, such as the connectivity of the 3D pores and some macroscopic properties from 2D SEM images. Newer techniques, e.g., focused ion-beam SEM (FIB-SEM) have been applied to overcome these challenges. For example, using tomography and FIB-SEM experiments, Botan et al. 172 developed a multiscale approach by employing statistical mechanics to study adsorption/transport in porous media while accounting for adsorption in various degrees of porosity. First, a lattice model is built for a given porous solid based on 3D FIB-SEM images with each site representing a porosity domain type, e.g., micro-, meso-, macro-, and non-porous. Subsequently, the

REVIEW scitation.org/journal/phf
lattice is sketched onto 3D data and the multiscale simulation is conducted to quantify fluid transport as a function of pressure gradient applied across the sample. The advantage of this multiscale model is to enable changes in adsorption and transport upon modifying simulation parameters. In addition, this approach does not require inferring any types of adsorption or flow so that the various phenomena can be taken into consideration. Nevertheless, the extraordinarily small size of samples and the costs for the use of the FIB-SEM technique hinder the broad applicability of the FIB-SEM in characterizing the porous media structure. 173 Recently, Tahmasebi et al. 174 suggested an effective stochastic approach that uses a few 2D images and creates 3D models or realizations of a porous network. It appears that the approach developed by Tahmasebi et al. 174 can yield insights into the complexity of shale reservoirs, as well as to other natural porous media with hierarchical pore structures, which could help in improving the accuracy of the KMC method. However, we point out that the accuracy of the stochastic method proposed by Tahmesebi et al. depends strongly on the availability of experimental data to reduce the uncertainty in the realizations. Once a representative 3D model is developed, both stochastic and deterministic approaches could be implemented, as discussed in the section titled Coupling MD and other deterministic simulation techniques.

Coupling MD and other deterministic simulation techniques
Tallarek et al. 175 presented an attractive bottom-up multiscale simulation approach that is deterministic in nature. Similarly to the work of Apostolopoulou et al., 87 Tallarek et al. 175 derived the effective bed diffusion coefficient (D bed ) taking into account for the solute dynamics as well as features such as solute properties, mobile-phase elution strength, and surface chemistry obtained from atomistic MD simulations. This information is incorporated into random-walk particle-tracking (RWPT) 176,177 -Brownian dynamics (BD) simulations used to extract the effective diffusion coefficient (see Fig. 8, top). In BD simulations, only trajectories and interactions between key particles are computed, while other components of the system are not included explicitly in the simulations but affect indirectly the dynamics of the Brownian particles via a random force. This reduces the dimensionality of the problem, making BD more computationally effective than the corresponding MD simulations. 178 To validate the accuracy of the RWPT approach, Tallarek et al. 175 compared the results of benzene diffusivity and density in a slit-pore model with a mixture of water/acetonitrile achieved from MD simulations, showing high levels of consistency. This suggests that the approach is applicable for estimating diffusion in systems with spatially dependent mobility [see Figs Other approaches exist. For example, Borg et al. proposed a sequential multiscale simulation method to accurately quantify macroscopic water flows through laboratory-scale carbon nanotube (CNT) membranes at moderate computational costs. 35 The multiscale simulation flow approach combines MD simulations with a continuum fluid model. Figure 9(a) shows an ideal laboratoryscale membrane of thickness L containing many, potentially billions, parallel CNTs. Fluid flow through the membrane can be described as follows: Δp = KmM, where K is the flow resistance through an arbitrary CNT and mM is the steady state mass flow rate through such CNT. The flow model through one CNT is decomposed into three components: one perfect long channel, one entrance, and one exit, as illustrated in Fig. 9(b). The flow resistance K is equivalent to the sum of the three correspondent flow resistances [ Fig. 9(c)]. Borg et al. conducted one MD simulation for the flow through the combined entrance/exit regions (subdomain 1) and another for the long developed-flow CNT region of length L ′ (subdomain 2) via applying periodic BCs, as shown in Fig. 9(d). The MD results were then fed as input parameters to a continuum flow model based on the Hagen-Poiseuille-Weissberg (H-P-W) equation. 72 Subsequently, the generated data are used to correct the H-P-W equation. This approach is easy to implement but still incorporates non-continuum fluid effects. Borg et al. 35 compared the results of the nanotube flow enhancement factor E as predicted from their H-P-W description against a range of experimental data, as shown in Fig. 9(e). The comparison shows that predictions of Borg et al. agree pretty well with some, [180][181][182] but not all the literature results. 183,184 Further studies are needed to resolve the discrepancy. The multiscale computational approach proposed by Borg et al. 35 appears to be applicable to complex membrane frameworks because it overcomes the drawbacks of the conventional H-P-W formalism, known to be not appropriate for non-ideal configurations.

SUMMARY, CONCLUSIONS, AND PERSPECTIVES
This Review focuses on multiscale fluid transport through heterogeneous porous matrices, with porosity ranging from the molecular (nm) to the macroscopic scale (>micrometer). The increased importance of such complex and multiscale porous materials has raised novel challenges regarding the fundamental understanding of fluid transport through them. In particular, the thermodynamics and dynamics of fluids confined in hierarchical porous media remain somewhat unclear. Among questions left unexplained are the role of interfaces between porosity domains, the breakdown of hydrodynamics at the nanoscale, and the interplay between adsorption and transport mechanisms at each pore scale. Addressing these questions will enable the rational design of hierarchical porous solids for many applications. We have discussed the applicability of each multiscale approach as well with its limitations for modeling fluid transport through heterogeneous porous matrices: • Atomistic equilibrium/non-equilibrium molecular dynamics (MD) simulations are an efficient technique that reveals a wealth of molecular-detailed properties for confined fluids. However, predicting the transport properties of fluids in hierarchical porous materials consisting of thousands of pores and ultimately within complex pore networks is at present unattainable with MD approaches. • Fluid transport properties in hierarchical porous materials can be measured by employing appropriate

Physics of Fluids
REVIEW scitation.org/journal/phf experimental tools such as pulsed field gradient nuclear magnetic resonance techniques and microimaging; notwithstanding, describing transport in these solids with a general framework remains challenging. Several models such as effective medium theory and simplified renormalization methods have been developed. However, despite significant progress, crucial questions remain to be addressed such as the role of surface energy barriers and the role of gas/liquid interfaces. Such effects can lead to complex activated transport mechanisms that remain poorly understood. When transport in different porosity scales is well-understood and robust homogenization techniques are developed, prediction of flow through macroscopic pore networks will be achieved using only simple parameters describing the structural features of the porous solid derived using limited datasets and others representing the fluid properties. Such understanding will launch applications in many sectors from energy to catalysis and environmental remediation. It should be stressed that the approach most fruitful for a given application strongly depends on the application itself, as different levels of accuracy are needed to predict, for example, the separation of multicomponent mixtures through a porous membrane or the gas permeability through a heterogeneous rock formation.

DATA AVAILABILITY
Data sharing is not applicable to this article as no new data were created or analyzed in this study.