A many-body dissipative particle dynamics model, namely, MDPD, is applied for simulation of pore-scale, multi-component, multi-phase fluid flows in fine-grained, nanoporous shales. Since this model is able to simultaneously capture the discrete features of fluid molecules in nanometer size pores and continuum fluid dynamics in larger pores, and is relatively easy to parameterize, it has been recognized as being particularly suitable for simulating complex fluid flow in multi-length-scale nanopore networks of shales. A remarkable feature of this work is the integration of a high-resolution FIB-SEM (focused ion beam scanning electron microscopy) digital imaging technique to the MDPD model for providing 3D voxel data that contain the invaluable geometrical and compositional information of shale samples. This is the first time that FIB-SEM is seamlessly linked to a Lagrangian model like MDPD for fluid flow simulation, which offers a robust approach to bridging gaps between the molecular- and continuum-scales, since the relevant spatial and temporal scales are too big for molecular dynamics, and too small for computational fluid dynamics with known constitutive models. Simulations ranging from a number of benchmark problems to a forced two-fluid flow in a Woodford shale sample are presented. Results indicate that this model can be used to deliver reasonable simulations for multi-component, multi-phase fluid flows in arbitrarily complex pore networks in shales.

## I. INTRODUCTION

Approximately 75% of the sedimentary rocks on earth are clastic nanoporous rocks, including compacted and diagenetically altered sands (tight sands) and a variety of rocks formed from fine grained sediments such as mudstones, siltstones, claystones, and marlstones, which are often referred to as “s*hale*s.” These shales contain most of the world’s oil and possibly most of the natural gas resources as well, according to a recent estimation (Kuuskraa *et al.*, 2013). However, only a small fraction can be recovered using the latest technology—typically about 5% for oil and about 20% for gas (Kuuskraa *et al.*, 2013). This low recovery is at least in part due to the lack of the knowledge and understanding of the relevant fundamental physics that ultimately control the dynamics of fluids organic-rich shales with extremely low permeabilities in the micro- to nano-Darcy range and average pore sizes from a few nanometers (10^{−9} m) to a few micrometers (10^{−6} m). Filling these knowledge gaps could potentially facilitate the development and deployment of more effective recovery strategies from tight shales.

Most of the theoretical understanding of fluid flow in porous geomaterials (and the predictive models built upon such understanding) has been based on the concepts of classical continuum fluid dynamics and a rigid porous or fractured solid porous matrix, which determines the “ideal” non-slip boundary conditions (BCs) for fluid flow and transport (Bear, 1972). For most of the part, those concepts and models have proven to be adequate for developing a fundamental understanding of a single- and multi-phase fluid flow in permeable porous media such as aquifers, soils, and conventional oil and gas reservoirs (Bertoncello and Honarpour, 2013). Many pore-scale models have been successfully developed in either Eulerian or Lagrangian frame, based on such continuum computational fluid dynamics (CFD). To list a number of recent works, Pan *et al*. (2001; 2004) applied a lattice Boltzmann (LB) and pore-network modeling approach to simulate a single- and two-phase flow in complex pore networks; Tartakovsky and Meakin (2005; 2006) introduced a smoothed particle hydrodynamics (SPH) model for immiscible and miscible flow simulations in fractured porous media, and Huang *et al*. (2005a; 2005b) applied a finite volume method based on the volume-of-fluid (VOF) model for multiphase flow simulations in fracture aperture and fracture junction with complex geometries. However, the behavior of fluids in organic-rich, fine-grained, nanoporous shales is very different from that in coarse-grained porous media. For example, in an idealized spherical pore with a diameter of 100 nm, ∼6% of the fluid is within a distance of 1 nm from the solid surface, and in a pore with diameter of 10 nm, ∼49% is within a distance of 1 nm; the properties of fluids in this thin surface layer are expected to be physically and chemically different from those of bulk fluids. Thus in organic-rich, fine-grained, nanoporous shales, the discreteness of atoms and molecules starts to become important for flow and transport processes at smaller scales; solid organic materials play an important role as mechanical components, sorbents and sources of fluids, and large specific surface areas make surface reactions and surface transport more important. Therefore, understanding large-scale fluid flow and transport behaviors within organic-rich, fine-grained, nanoporous shales requires robust and accurate multiscale computational models that can bridge both the spatial and temporal scales between fluid molecular dynamics and pore-scale flow.

Originally introduced by Hoogerbrugge and Koelman (1992) for microscopic hydrodynamics, the dissipative particle dynamics (DPD) models constitute a relatively new class of mesoscale models that can be used to simulate a single- and a multi-phase fluid flow, e.g., see Tiwari and Abraham (2006), Heldele *et al.* (2006), Visser *et al.* (2006), and Liu *et al.* (2006). A theoretical basis for DPD based on statistical mechanics was discussed in Español and Warren (1995) and Marsh (1998). A recent review of DPD theory, methodology, recent developments, and applications can be found in Moeendarbary *et al.* (2009) and Liu *et al.* (2015). Using DPD, a complex system can be simulated with a set of interacting particles, where each particle represents a small cluster of atoms or molecules instead of a single molecule (Figure 1). The particle-particle interaction force in a DPD embodiment consists of a “conservative” (non-dissipative) component, a random (fluctuating) component, and a dissipative component that represents the effect of viscosity. The inclusion of random particle-particle interactions in addition to the conservative particle-particle interactions is the most remarkable feature that distinguishes DPD from the SPH method. The random force generates fluctuations that correspond to thermally driven fluctuation in real fluid molecules. Such fluctuations are only important on microscopic physical scales, e.g., pores of sizes in nanometer ranges. Conversely, the DPD fluid can reproduce the continuum Navier-Stokes equations only on large scales (scales much greater than the particle size) with the effect of thermal fluctuations being negligible. This makes DPD essentially a mesoscale method between the molecular and continuum hydrodynamic scales, and facilitates simulations of complex fluid systems scales spanning a wide range of length and time scales. DPD rigorously conserves the number of particles (equivalently the mass) and momentum of the system. In addition, much larger time steps are possible in DPD simulations than in MD simulations. This is because the DPD particle-particle interactions are much softer than those (e.g., Lennard-Jones potential) experienced in MD. For these reasons, DPD is more efficient than MD for mesoscale flow simulations in terms of computational costs. Among the many existing DPD models, the “many-body” DPD model by Warren (2003), namely, MDPD, has been recognized as being particularly suitable for multi-fluid systems and relatively simple to parameterize for specific types of fluids. This MDPD model has been applied in simulations of liquid-vapor interface, surface tension, and multi-component fluid flows in micro-scale channels, e.g., Pan (2010), Chen *et al.* (2011), Ghoufi and Malfreyt (2011), Chen *et al.* (2013), and Chen *et al.* (2014). More recently, Liu *et al.* (2007a; 2007b) successfully applied a DPD model using an SPH-like conservative interaction force to simulate unsaturated flows in porous media. However in their works, only pair-wise particle-particle interaction potential and single phase were considered.

The major effort of the present work is the application of a robust MDPD model for simulations of pore-scale, multi-component, multi-phase fluid flows within fine-grained, nanoporous shales. This is a relatively unexplored area so far. The novelties of the work lie in the following two aspects. First, this model is able to and has been applied to capture simultaneously the discrete features of fluid molecules in nanometer size pores and continuum fluid dynamics in larger pores. An example in Figure 2 demonstrates the multi-length-scale modeling capability of MDPD, where the same particle-particle interaction force and parameters have been used to model fluid-fluid interfaces in two cases: (a) at the micrometer-scale where the interface is relatively stable; and (b) at the nanometer scale where the interface is diffused due to the Brownian motion of fluid particles. Second, an FIB-SEM (focused ion beam scanning electron microscopy) high-resolution digital imaging technique (Goral and Miskovic, 2015) has been utilized to provide the MDPD model with 3D voxel data that contain the invaluable micro-geometrical and compositional information of shale samples. The process of the digital imaging process is briefly discussed in the Appendix for clarity. The digitized microstructural information is used to construct pore networks with solid boundaries that enclose pores. To the best of the authors’ knowledge, this is the first time that FIB-SEM has been seamlessly integrated to a Lagrangian particle transport model like MDPD for modeling fluid flows in real multi-length-scale nanoporous shale samples, as compared with earlier DPD fluid flow simulations where only manufactured (or simulated) pore networks were used (Liu *et al.*, 2007b). It is worth of mentioning that although the FIB-SEM machinery has become a popular approach for characterizing shale rock samples since the works by Curtis *et al.* 2010; 2011; 2012a; 2012b, most of the earlier flow and transport simulations in shales were based on continuum CFD models. For example, Dewers *et al.* (2012) reported a finite element simulation of shale gas transport in Hayneville Formation mudstones by assuming sufficient interconnected large pores in the samples and that Darcy’s Law was valid. Unfortunately, those models are not appropriate for heterogeneous nanoporous shale rocks any more. In comparison, the integrated research method to be presented in this work offers a robust approach to bridging gaps between the molecular- and continuum-scales, as the relevant spatial and temporal scales are too big for molecular dynamics, and too small for continuum CFD with known constitutive models.

The rest of the paper is organized as follows. First, the MDPD model is presented in Section II. Then a number of validation test cases are reported in Section III. A pioneering flow simulation in a pore network enclosed in an organic-rich, nanoporous shale is presented in Section IV. Finally, the concluding remarks are given in Section V.

## II. MDPD PORE-SCALE FLUID FLOW MODELING METHODOLOGY

### A. Mathematical description of DPD and MDPD models

In a generic DPD model, particles interact via pairwise central forces, i.e., $\mathbf{F}ij=\mathbf{F}ijR+\mathbf{F}ijD+\mathbf{F}ijC$, where $\mathbf{F}ijR$ represents a random force, $\mathbf{F}ijD$ a dissipative force, and $\mathbf{F}ijC$ a conservative force between particle *i* and *j*, respectively. If **r**_{i} and **v**_{i} are used to denote the position and velocity of particle *i*, respectively, the random force, $\mathbf{F}ijR$, and the dissipative force, $\mathbf{F}ijD$, can be expressed as $\mathbf{F}ijR=\sigma wR(rij)\xi ij\mathbf{r}^ij$, and $\mathbf{F}ijD=\u2212\gamma wD(rij)(\mathbf{r}^ij\u22c5\mathbf{v}ij)\mathbf{r}^ij$, where **r**_{ij} = **r**_{i}−**r**_{j}, $rij=|\mathbf{r}ij|$, $\mathbf{r}^=\mathbf{r}ij/rij,$ **v**_{ij} = **v**_{i} − **v**_{j}. These forces constitute a thermostat if the amplitude, $\sigma $, of the random variable, $\xi ij$, and the viscous dissipation coefficient, $\gamma $, satisfy a fluctuation-dissipation theorem: $\sigma 2=2\gamma kBT$ and $wD(r)=(wR(rij))2$, where *k*_{B}*T* denotes the desired temperature in the unit of Boltzmann’s constant *k*_{B}. In the original DPD model, the conservative force, $\mathbf{F}ijC$, is defined as $\mathbf{F}ijC=aijwC(rij)\mathbf{r}^ij$, where *a*_{ij} denotes the magnitude of the force, and the weight function *w*^{C}(*r*) vanishes when the inter-particle distance *r* is larger than a cutoff range *r*_{c}. $\mathbf{F}ijC$ is usually derived from a soft and unspecific weight function, *w*^{C}(*r*_{ij}), thus allowing for a fairly large integration time step. Different weight functions describe different material properties. A common choice for *w*^{C}(*r*_{ij}) is *w*^{C}(*r*_{ij}) = 1 − *r*_{ij}/*r*_{c}, and *w*^{R}=*w*^{C}. The standard velocity Verlet algorithm (Revenga *et al.*, 1998) is employed to integrate the resulting equations of motion. A quadratic equation of state (EOS) was obtained with respect to the average particle density *ρ*, as shown in Figure 3(a). However, the original DPD model is not sufficient to model multi-species or multi-phase fluid flow phenomena such as those associated with liquid-vapor interfaces, liquid-liquid interfaces, and free capillary surfaces. A more complex EOS needs to be represented with the DPD method. To achieve this, a long-range attractive, short-range repulsive conservative force, **F**^{C}, is required. The approach employed in the present work is the so-called many-body DPD method (Warren, 2003), namely, MDPD. The $\mathbf{F}ijC$ in the standard DPD method is augmented by density-dependent contributions, and the resulting method includes the van der Waals loop in the EOS, as shown in Figure 3(b). In the MDPD model, the conservative force is expressed as

which consists of a long-range attractive part that is density-independent, and a short-range repulsive part that depends on a weighted average of the local particle density. The attractive component, *A*_{ij}*w*^{C}(*r*_{ij})$\mathbf{r}^ij$, can be obtained by simply turning the sign of the original force parameter *a*_{ij}, i.e., *A*_{ij} < 0, with a cutoff range *r*_{C} = 1, and $Bij(\rho \xafi+\rho \xafj)wd(rij)\mathbf{r}^ij$ is a many-body repulsive component with *B*_{ij} > 0, and shorter cutoff *w*_{d}(*r*_{ij}) = 1 − *r*/*r*_{d}, where *r*_{d} < *r*_{C}. The averaged local density, $\rho \xafi$, at the position of particle *i* can be computed as $\rho \xafi=\u2211j\u2260iw\rho (rij)$, where the normalized weight function, $w\rho $, needs to satisfy $\u222b0\u221e4\pi r2w\rho (r)dr=1$. $w\rho $ is defined as follows for a three-dimensional computational domain: $wd(r)=152\pi rd3(1\u2212rrd)2$.

### B. Treatment of solid-fluid interface

Because of the soft interaction between DPD particles, fluid particles may penetrate through a solid matrix when DPD and MDPD models are used for pore-scale fluid flow simulations. Such penetration is not physical and must be avoided. Most of the previous efforts to properly model solid-fluid interfaces in DPD fluid flow simulations were focused on imposing rigorous macroscopic boundary conditions (BCs), e.g., a non-slip BC at sharply defined impenetrable solid surfaces. That motivation originated from a strict mesoscale interpretation of DPD models, where a single DPD fluid particle represents a cluster of fluid molecules on scales well above the atomistic scale (Groot and Warren, 1997). In order to make the non-slip BC possible, additional forces must be exerted on fluid particles at the vicinity of solid-fluid interfaces. Parameterization must be carried out very carefully in order to avoid undesired spurious behavior, e.g., artificial slip (Pivkin and Karniadakis, 2005), oscillations in temperature (Revenga *et al.*, 1999), and particle layering (Pivkin and Karniadakis, 2006). The disadvantages of those specialized boundary models is that they are only applicable to mostly flat, or spherically curved solid surfaces, e.g., Liu *et al.* (2007b). To overcome the difficulty of the rigorous non-slip BC, Henrich *et al.* (2007) devised a boundary model, where a weak external repelling force is imposed on those fluid particles whenever they penetrate into a thin layer in the solid boundary. This prevents fluid particles from further diffusing into solid matrix indefinitely.

However, in nanometer-sized pores, surface slip might be an important transport mechanism as well. Thus in order to account for the multi-length-scale nature of pores in shales, a boundary model is desired to model the continuum no-slip boundary in large pores, the non-continuum discrete motions of fluid particles in small pores, and the transition between the non-continuum and continuum fluid motions.

In this work, a simple boundary model is devised for modeling solid-fluid interfaces in the MDPD model. In this model, a dense layer of solid particles is defined with the MDPD particle attributes, and used to enclose pore space. Distribution of the solid particles is configured to be sufficiently dense to prevent penetration of the MDPD fluid particles that flow in pore space. This is equivalent to the application of a weak external repelling force approach (Pivkin and Karniadakis, 2006). The coordinates of solid particles are determined using 3D voxel data that are obtained through a digital rock imaging process, and these coordinates are fixed during simulations. For clarity, an example of data conversion from a 2D pixel map to coordinates of solid particle is illustrated in Figure 4. Later in Section IV, 3D voxel data of a shale sample obtained from the FIB-SEM digital rock imaging process will be used. For MDPD simulations, Eq. (1) is applied for (a) fluid-fluid, (b) solid-solid, and (c) solid-fluid particle interactions. However, the coefficients in Eq. (1) must be parameterized for each of the three interaction modes, so that fluid properties at the pore scale under consideration can be properly represented, e.g., liquid-vapor interface, miscibility between two types of fluids, and fluid wettability against the solid surface. Notice that this is not a completely new approach for modeling solid boundaries in DPD-type methods; for example, in their work, Pivkin and Karniadakis (2006) mentioned a similar idea, with which though they did not demonstrate a test case. Using the present approach, however, there is a potential artifact of undesired minor temperature oscillations in the vicinity of solid-fluid interfaces. Nevertheless, this simple boundary model is probably so far the only viable way for MDPD to model arbitrarily complex networks of enclosed pore space in nanoporous shales. As already shown in Figure 2, this boundary model is able to render a diffusive solid-fluid interface in a nanometer-scale pore, and also a stable solid-fluid interface in a micrometer-scale pore.

## III. MODEL VALIDATION

The MDPD model described above has been programmed in an in-house code based on the LAMMPS framework (Plimpton, 1995). Before the code could be reliably applied for problems of practical interest, code validation was performed in this section to validate the implementation and key capabilities of the model for simulating multi-component, multi-phase fluid flows. Three test cases are presented, with the particle mass, *m*, temperature, *k*_{B}*T* and many-body attraction cutoff range, *r*_{C}, chosen to be the reduced units of mass, energy, and length, respectively, i.e., *m* = (*k*_{B}*T*) = *r*_{C} = 1. In the simulations, a constant time step *dt* = 0.01 in the reduced time unit along with the constant-energy, constant-volume NVE ensembles was used. The gravity force was neglected.

### A. Liquid-vacuum interface

The objective of this test case is to validate if the implemented MDPD code accurately calculates properties of a specific type of fluid. We followed Ghoufi and Malfreyt (2011) to simulate a water liquid-vacuum interface, where the water density and surface tension were calculated. A 3D computational domain was defined, extending from −10*r*_{c} to 10*r*_{c} in the x-direction, from 0 to 5*r*_{c} in the y-direction, and from 0 to 5*r*_{c} in the z-direction. A periodic BC was used in all the directions. 1000 DPD particles were then randomly placed in *x* = [−5*r*_{c}, 5*r*_{c}]. The parameters in Eq. (1) are set as *A*_{ij} = −50, *B*_{ij} = 25, and *r*_{d} = 0.75*r*_{C} for modeling the water (Ghoufi *et al.*, 2010). In this case, one DPD particle approximately represents a cluster of three water molecules, i.e., *N*_{m} = 3, and the size of one DPD particle is about 90 Å^{3}. More details of the correspondence between the DPD reduced units and physical values can be found in Ghoufi and Malfreyt (2011). To calculate the thermodynamic properties, 1 × 10^{6} time steps were carried out to equilibrate the system. Figure 5 displays an instantaneous snapshot of the system after equilibrium. Finally, 1 × 10^{6} additional time steps were then conducted to calculate the time-averaged properties.

Figure 6 shows a time-averaged density profile as a function of the x-axis. The density near *x* = 0 (the center of the water drop) is 6.88, matching the value reported in an MD study (Ghoufi and Malfreyt, 2011). Due to the simple geometry in this problem, the water-vacuum interfacial tension, $\gamma WV$, can be calculated by subtracting the mean tangential stresses, $\sigma yy$ and $\sigma zz$, from the normal stress, $\sigma xx$:$\gamma WV=Lx\u27e8\sigma xx\u221212(\sigma yy+\sigma zz)\u27e9$. The calculated $\gamma WV$ is 12.4, matching the value reported in Ghoufi and Malfreyt (2011). The values for water density and water-vacuum interfacial tension can be converted into the physical units with the equations: $rc=rc*(\rho *NmV)13\xc5$, $\rho =\rho *NmMNarc3kgm\u22123$, and $\gamma =\gamma *kBTrc2Nm\u22121$, where the superscript * denote values in the DPD reduced unit, *V* is the volume of one water molecule (30 $\xc5$), $M$ is the molar weight of a water molecule (18 g mol^{−1}), *N*_{a} is Avogadro’s number, and *k*_{B} is Boltzmann’s constant, and *T* is 298 K. In the physical units, $\rho =994kgm\u22123$ and $\gamma =70.6\xd710\u22126Nm\u22121$, which agree well with the MD results in Ghoufi and Malfreyt (2011). This indicates that the implemented MDPD code delivers accurate predictions of thermodynamic properties for fluids of interest.

### B. Interfacial tension between two fluids

Interface stability and interfacial tension between two types of fluids are crucial properties for characterizing fluid-fluid miscibility. To model and simulate multi-fluid flows, it is desirable to have a simple approach for independently calculating the interfacial tension between liquid-liquid pairs. In this numerical experiment, we set up the simulation problem that is shown in Figure 7. It contains an 11 × 10 × 10 non-wetting liquid slab with 4000 particles, and an 11 × 10 × 10 wetting liquid slab with 4000 particles. These slabs are initially attached to each other. The total tension of the system, $\gamma total$, consists of three components: the wetting liquid surface tension, $\gamma WV$, the non-wetting liquid surface tension, $\gamma NV$, and the interfacial tension between the wetting and non-wetting fluids, $\gamma WN$. For this nomenclature, the subscripts “W,” “N,” and “V” denote the wetting fluid, non-wetting fluid, and vacuum, respectively. $\gamma WN$ can be calculated by subtracting the mean tangential stresses, $\sigma yy$ and $\sigma zz$, from the normal stress, $\sigma xx$, i.e., $\gamma total=Lx\u27e8\sigma xx\u221212(\sigma yy+\sigma zz)\u27e9$, where *L*_{x} is the length of the simulation box in the x direction. Since $\gamma WV$ and $\gamma NV$ remain unchanged, $\gamma WN$ can then be calculated by subtracting the two surface tensions from the total tension, $\gamma total$, i.e., $\gamma WN=\gamma total\u2212\gamma WV\u2212\gamma NV$, where $\gamma WV=\gamma NV=7.67$ is obtained when there is only one fluid type in the system.

The MDPD parameters used in this simulation are as follows: *A*_{WW} = −40, *B*_{WW} = 25, and *r*_{d} = 0.75*r*_{C} between the wetting fluid particles; *A*_{NN} = −40, *B*_{NN} = 25, and *r*_{d} = 0.75*r*_{C} between the non-wetting fluid particles; and *B*_{WN} = 25 and *r*_{d} = 0.75*r*_{C} between the wetting and non-wetting fluid particles. By changing *A*_{WN}, different $\gamma WN$ values can be obtained, leading to different miscibility between the two liquids; see Figure 8. By increasing the *A*_{WN} value from −40 to −5, that is, turning the attractive interaction between the wetting and non-wetting liquids from strong to weak, the liquid–liquid interface changes accordingly from being completely miscible to entirely immiscible. Notice that before *A*_{WN} approaches 0, e.g., when *A*_{WN} = −5, the two liquids already tend to detach from each other since the inter-liquid attraction becomes too weak between the two liquids. A correlation between $\gamma WN$ and *A*_{WN} is shown in Figure 9. This matches the result in Chen *et al.* (2011). In general, successively altering the attraction parameter, *A*, while fixing the repulsion parameter, *B*, is an adequate approach for tuning the interaction of concern. Similarly, a solid–liquid interaction can also be tuned by adjusting the attraction parameter, *A*, from −40 to 0, which will result in an increase of solid–liquid interfacial tension as well as a decline of wettability of the liquid on the solid.

### C. Static fluid in nanometer-scale slit pore

In this test case, a static fluid residing in a nanometer-scale slit pore was simulated using MDPD. The objective is to demonstrate the ability of MDPD to model the static contact angle between a single fluid and bounding solid matter. This is an important requirement for characterizing the wetting properties of fluids in organic-rich, nanoporous shales. MDPD can model the interaction force between solids and liquids by choosing a proper set of parameters, *A*_{SL}, *B*_{SL}, and *r*_{d} in Eq. (1). The subscript “S” and “L” denote the solid and liquid, respectively. For simplicity, we fixed *B*_{SL} = 25 and *r*_{d} = 0.75*r*_{C}, and selected six *A*_{SL} values to study the dependence of the wettability of the liquid on the solid matter. To model liquid-liquid interactions, the MDPD parameters used in this simulation are *A*_{LL} = −50, *B*_{LL} = 25, and *r*_{d} = 0.75*r*_{C}. Notice that in a similar simulation (Pan, 2010), different choices of values for *r*_{C} and *r*_{d} in Eq. (1) were used. In the present test problem, the simulation domain extended from −30 to 30 in the x direction, from −5 to 5 in the y direction, and from −2.5 to 2.5 in the z direction. Periodic BCs were prescribed in the x and z directions.

To form the solid matter of the slit pore, 4500 solid particles were placed randomly in the two regions of y = [−5, −4] and y = [4, 5], respectively. These were run to thermodynamic equilibrium in these regions using the original DPD method, and then frozen for the rest of the simulations. The width of the slit pore was 8*r*_{C}, corresponding to 8.616 nm in physical units. Next, 4000 liquid particles were randomly placed in a region bounded by x = [−13, 13] and z = [−4, 4], and run for 10 000 time steps to thermodynamic equilibrium using the MDPD method. The bottom and top solid particle layers of the domain were dense enough to prevent the liquid particles from penetrating indefinitely.

The instantaneous particle distributions corresponding to *A*_{SL} ranging from −40 to −5 are shown in Figure 10. The larger is the absolute value of *A*_{SL}, the more wetting is the liquid in the slit pore. For example, the case of *A*_{SL} = −40 corresponds to a fully wetting effect with a contact angle of *θ* = 0°. This is shown in Figure 10(a). As the absolute value of *A*_{SL} was successively reduced, the liquid gradually turned from wetting to non-wetting; see Figures 10(b)–10(e). Finally, when *A*_{SL} = −5, a fully non-wetting effect with a contact angle of *θ* = 180° is observed in Figure 10(f). Overall, the simulation results presented in this test problem have demonstrated that the MDPD method can be used to conveniently model the contact between solid matter and liquids.

## IV. FLUID FLOW IN WOODFORD SHALE ORGANIC-MATTER-HOSTED PORE NETWORK

In this section, we assess the developed capability for problems of practical interest. A workflow that integrates high-resolution digital rock imaging to the MDPD model was utilized for simulating a two-component, two-phase fluid flow in heterogeneous, organic-rich, nanoporous shale. The choice for a region of interest in a shale rock sample and the corresponding fluid flow simulation results are discussed as follows.

### A. Region of interest

Figure 11 displays FIB-SEM reconstructed images of a (5000 nm)^{3} domain obtained in a Woodford shale sample, which shows three of the main compositions: organic matter, inorganic matter, and pores enclosed in organic matter. As seen in Figure 11(c), most of the pores enclosed in organic matter are relatively small, isolated, and not connected to others. This indicates that no fluid flow is possible between those pores. A fluid flow is only possible in those inter-connected larger pores, which are therefore the candidates of region of interest for simulations. However, a domain like that in Figure 11 is prohibitive for MDPD simulations in terms of storage and computing time. Instead, a smaller subdomain is considered to be representative if it contains sufficient details of interconnected pores with both large and small pore throats enclosed in the organic matter. Notice that it is not trivial to locate such a region in Woodford shale samples. In many individual samples, we found no satisfying region at all. The sample shown in Figure 11 was among the many that we inspected, and was the one we eventually chose.

We identified a 440 nm × 950 nm × 1600 nm subdomain as a region of interest from the domain in Figure 11, and then used it to prepare the simulations. This subdomain is shown in Figure 12(a), where surface geometries of an organic-matter-hosted pore network are represented by a finite element mesh generated based on the 3D voxel data. The largest pore was denoted as “primary pore” in the rest of the text. Rich details of interconnected large and small pore throats can be seen in the primary pore. It is perhaps only part of a larger pore network, but geometric details of this pore are complex enough for demonstrating fluid flows in nanometer-sized pores. In addition, there is also a large isolated pore in this region, through which we assumed no fluid flow. Figure 12(b) shows 590 608 MDPD solid particles that are used to enclose pore space in the region of interest. The particle coordinates were determined using the voxel map and converted into LAMMPS data format using the procedures described earlier. As one can see, the geometric details of pores are different at multiple locations in this region between Figures 12(a) and 12(b). For example, one can see much more small isolated pores in Figure 12(b), but not in Figure 12(a). This is because the finite elements may not preserve all of the geometric information from the voxel map, depending on the minimum element size specified during the mesh generation. This also indicates the possibility that some small pore throats connecting bigger pores could have been ignored in a finite element mesh, resulting in a false representation of the original geometry of pore networks. However, in the MDPD model, no geometric information was lost by conversion from the voxel map to MDPD particles.

### B. Simulation setup

The values of the attractive interaction parameter, *A*, for fluid-fluid and fluid-solid interactions are listed in Table I. The Type-1 fluid was defined as a wetting liquid that initially resided in the primary pore, and the Type-2 fluid was defined as a non-wetting liquid to be injected in the primary pore in order to extract the Type-1 fluid. The rest of the MDPD parameters used in this simulation problem are *B*_{ij} = 25 and $rd=0.75rC$ for all particle-particle interactions between the same and different types. As shown in Figure 12(a), the flow inlet boundary was defined at the open pore throats in the lower left corner of the primary pore, and the flow outlet boundary was defined at the open pore throats in the top right corner. In order to get an initial condition such that the primary pore was fully saturated with the Type-1 fluid, a total of 250 000 Type-1 fluid particles were injected into the empty primary pore with a reflection BC imposed at the flow inlet and flow outlet boundaries to prevent fluid particles from exiting the pore. Additionally 5000 time steps were then run to equilibrate the system. This resulted in a highly pressurized pore space with an averaged fluid pressure of 2.05 in the reduced pressure unit. The coordinates of Type-1 fluid and solid particles were then saved and used as initial conditions for the simulations of Type-2 fluid injection that follow. The injection of Type-2 fluid particles into the primary pore was designed to last for a total time of 30 000 in the reduced time unit, which is roughly equivalent to 20 $\mu $s. Simulations were conducted with three Type-2 particle injection rates:

Test 1, one particle every five time steps;

Test 2, one particle every two time steps;

Test 3, one particle every time step.

A
. | Solid . | Type-1 fluid . | Type-2 fluid . |
---|---|---|---|

Solid | $\cdots $ | −40 | −10 |

Type-1 fluid | −40 | −40 | −20 |

Type-2 fluid | −10 | −20 | −40 |

A
. | Solid . | Type-1 fluid . | Type-2 fluid . |
---|---|---|---|

Solid | $\cdots $ | −40 | −10 |

Type-1 fluid | −40 | −40 | −20 |

Type-2 fluid | −10 | −20 | −40 |

Upon an approximate conversion of units to the real physical units, these injection rates correspond to ∼2 × 10^{−12} g/s, ∼5 × 10^{−12} g/s, and ∼10 × 10^{−12} g/s, respectively. In this problem, an accurate calibration to match specific types of fluids like in Section III A is out of our scope. To ensure that fluid particles did not arbitrarily exit the pore, the reflection BC was maintained at the flow inlet and flow outlet boundaries. Meanwhile, an ejection BC was used to randomly delete any fluid particle that was very near the outlet boundary at a removal rate equal to the corresponding injection rate. This treatment guaranteed a constant total mass of fluids in the primary pore, and also a roughly constant average fluid pressure during the injection. Notice that this resembled a possible scenario of producing from nanoporous, kerogen-rich shale at reasonable flow rates.

A total of 240 CPU processors (2.66 GHz for each one) were used in parallel computing for each injection test, and the total wall time of each test was about 270 h. Unavoidably, the parallel scalability of the simulations was largely impaired by operations like particle injection at the inlet and particle removal at the outlet, as well as the highly irregular spatial distribution of the geometry. Regular load rebalancing algorithms by re-partitioning subdomains according to the updated particle distribution did not help much in this case. This MDPD system appears containing not a huge number of particles as compared to large-scale MD simulations. If we take the Type-II kerogen (C_{1512}H_{1764}N3_{6}O_{144}S_{18}) to be the main composition in the organic matter as an example, where six of such kerogen structures occupy a volume of (3.25 nm)^{3} (Pawar and Huang, 2016), then an equivalent MD system may contain 3 543 648 kerogen structures and 750 000 fluid molecules. However, the possible physical time evolution by such an MD simulation is severely constrained by the allowable time step size that is at the scale of femtosecond (10^{−15} s), and by the tremendous MD potential calculations in each step. Therefore, so far it is still practically impossible for MD simulations to model an equivalent system of such size.

### C. Simulation results

Figures 13–15 show a series of instantaneous fluid particle distribution simulated in the three injection tests. Type-1 and Type-2 fluid particles were colored yellow and blue, respectively, and the size of Type-2 fluid particles was set larger to highlight the flow path of the injected fluid. Four scenes are displayed in each figure, with each scene representing a typical moment when the same amount of Type-2 fluid particles had been injected between the three tests. For example, a moment when 1000 Type-2 fluid particles had been injected is depicted in Figures 13(a), 14(a), and 15(a) for each test, respectively. Distribution of the Type-2 fluid particles looks similar in those scenes. Likewise, the scenes in (b)–(d) of each figure look like their counterparts in the other two figures. The equivalence in those flow paths of the injected fluid indicates that the evolution of preferential flow path in the primary pore is relatively independent of the injection rate.

Figure 16 displays the percentage of the Type-1 fluid extracted from the primary pore (ratio of recovery) versus the elapsed injection time. At the beginning, the ratio of recovery for the Type-1 fluid grew linearly with respect to time before the preferential flow path of Type-2 fluid was not fully formed. Once the preferential flow path of Type-2 fluid was developed with the ratio of recovery of the Type-1 fluid at about 60%, the recovery of the Type-1 fluid went to stagnation quickly regardless of the injection rate. In addition, Figure 17 shows the extraction speed of the Type-1 fluid (rate of recovery) versus the elapsed injection time. Notice that the rate of recovery for the Type-1 fluid dropped dramatically soon after the injection of the Type-2 fluid was started, despite the injection rates. The simulation results indicate that it may not be very challenging to reach a relatively low ratio of recovery for the fluid stored in the pores (60% or less). However, efficiency of the recovery is extremely low if 80% or more is pursued. This indicates that smarter geo-techniques have to be developed for further improvement of production. Moreover, the ratio of recovery estimated from the simulations is noticeably higher that those reported in oil and gas industry, e.g., 5% for oil, and about 20% for gas (Kuuskraa *et al.*, 2013). One also has to distinguish the difference that those numbers reported in literature are based on the reservoir-scale production, rather than a tiny representative pore network used in our simulations. Thus from the simulations as well as the data in literature, we could infer that the majority of oil or gas stored in the pores enclosed by organic matters was not recovered from shale reservoirs with existing geo-engineering processes because most of the large pore networks are likely to be stay disconnected and isolated from others in heterogeneous, nanoporous shales (as evidence is already shown in Figure 11(c)). Therefore, it would be no surprise to see those actual low ratios of recovery in industry, if no effective fracturing is performed to create sufficient connections between those large pore networks. This also indicates that more advanced fracturing techniques may be needed.

In conclusion, this simulation test problem has completely demonstrated the advantages and robustness of a coherent integration from a high-resolution digital rock imaging process to the MDPD flow model, which enables pore-scale, multi-component, multi-phase fluid flow simulations in real fine-grained, nanoporous shale samples.

## V. CONCLUSIONS

A many-body dissipative particle dynamics (MDPD) model has been applied for simulating multi-phase fluid flow in nanopore networks enclosed in fine-grained, nanoporous shales. A distinguished part of this work that sets it apart from the previous ones (Liu *et al.*, 2007a; Liu *et al.*, 2007b and Chen *et al.*, 2011) is the integration of high-resolution FIB-SEM digital imaging technique, which provides 3D voxel data that contain invaluable geometrical and compositional information of shale samples. Although applications of FIB-SEM to studying shales are not new, this is the first time that FIB-SEM has been seamlessly integrated to a Lagrangian particle transport model like MDPD for modeling fluid flows in multi-length-scale nanoporous shale samples, as compared with earlier DPD fluid flow simulations where only manufactured (or simulated) pore networks were used. Since the relevant spatial and temporal scales are too big for molecular dynamics, and too small for computational fluid dynamics with known constitutive models, this integrated research method offers a robust approach to bridging gaps between the molecular- and continuum-scales. Additional remarks worth of noting are as follows.

First, it is not possible for MDPD to properly simulate a fluid dynamics system that contains dissimilar fluids (Ghoufi and Malfreyt, 2011), e.g., water and supercritical CO_{2}, and methane and supercritical CO_{2}. However, those are important scenarios in fine-grained, nanoporous shales. New conservative force formulation must be explored for the MDPD to treat dissimilar fluids. Second, it is computationally prohibitive to use MDPD to simulate a shale sample when the characteristic length of the sample extends from tens of microns to several millimeters. Computer memory required for the digital rock structure reconstruction and data storage could exceed tens of terabytes for a large shale specimen (Goral and Miskovic, 2015). Despite the known limitations, the combination of FIB-SEM and MDPD model offers an unparalleled approach for investigating flow physics in nanopore networks enclosed in organic-rich, heterogeneous, nanoporous shales.

## ACKNOWLEDGMENTS

This work is supported by the Idaho National Laboratory (INL) Laboratory Directed Research and Development (LDRD) Program under DOE Idaho Operations Office Contract No. DE-AC07-05ID14517. The authors wish to acknowledge helpful reviews and comments from Ian Walton and John McLennan from Energy and Geoscience Institute. We also thank Jeff Gelb and Jack Kasahara from Carl Zeiss Microscopy for supporting much of this work. Additionally, the authors wish to acknowledge the help and support of Earl Mattson from INL. We are also grateful to U.S. Geological Survey (USGS) Core Research Center for providing access to core materials.

### APPENDIX: A WORKFLOW FROM DIGITAL ROCK IMAGING TO NUMERICAL MODELING

A digital-rock-physics (DRP) technique recently discussed in Goral and Miskovic (2015) is employed to enable a representative and uniquely informative perspective on the structures and chemical compositions of shales. The procedures are presented as follows.

Nondestructive micro and nano X-ray microscopy is applied in order to identify a region of interest (ROI) within a shale sample. The selected ROI is then imaged with destructive FIB-SEM serial sectioning.

The FIB-SEM data are processed, segmented, and reconstructed into digitized information of organic matter, nonorganic (mineral) matter, and pore space in either of the two phases and between.

The ROI of the organic-matter-hosted pore system within the shale rock matrix is identified and separated. The resulting pore network geometry is then gridded with a tetrahedron volume mesh to enable further alternative modeling and simulation studies.

An example of workflow from DRP to numerical modeling is displayed in Figure 18. Notice that image post-processing/segmentation is not a straightforward process, which highly relies on image interpretation. Improper image processing and segmentation can lead to, for example, misidentification of organic matter as pore space, resulting in overestimation of porosity and permeability.