A Molecular Dynamics Investigation on Methane Flow and Water Droplets Sliding in Organic Shale Pores with Nano-structured Roughness

Roughness of surfaces significantly influences how methane and water flow in shale nanopores. We perform molecular dynamics simulations to investigate the influence of surface roughness on pore-scale transport of pure methane as well as of two-phase methane–water systems with the water sliding as droplets over the pore surface. For single-phase methane flow, surface roughness shows a limited influence on bulk methane density, while it significantly reduces the methane flow capacity. In methane–water systems, the mobility of water is a strong function of surface roughness including a clear transition between immobile and mobile water droplets. For cases with mobile water, droplet sliding speeds were correlated with pressure gradient and surface roughness. Sliding water droplets hardly deform, i.e., there is little difference between their advancing and receding contact angle with structured roughness.


Introduction
Shale gas has recently reshaped the oil and gas production landscape and has been recognized as the potential alternative for the increasingly exhausted conventional geoenergy resources (Vedachalam et al. 2015). Shale gas is mainly located in micropores and mesopores with pore diameters ranging from 1 to 50 nm (Cao et al. 2015;Clarkson et al. 2013).
Shale rocks consist of inorganic minerals, such as quartz and clays, and organic matter. The main deposition of the heavy insoluble organic matter is kerogen (Loucks et al. 2012). It has extremely low permeability and porosity. For a shale sample that is primarily formed of kerogen, a significant fraction of its porosity lies in 2-7 nm pores and the characteristic pore throat size of 6 nm accounts for the largest fraction (Sakhaee-Pour and Bryant 2012).

3
Kerogen is considered as the origin of gas in shale reservoirs and its organic nanopore serves as primary methane storage space. Therefore, kerogen determines not only the production of shale gas but also the gas migration mechanism (Tesson and Firoozabadi 2018;Wu et al. 2019a, b).
Kerogen exhibits rough porous structure connected by abundant nanometer-sized pores (less than 10 nm) (Loucks et al. 2009). Javadpour et al. (2015) measured nanoscale roughness in kerogen through atomic force microscopy and pointed out that surface roughness significantly influences shale gas transport. Therefore, a precise description of the transport characteristics of shale gas in rough nanopores is relevant for exploiting shale gas resources.
There is a significant amount of literature on the influence of surface roughness on fluid flow. Through a delicately designed set of experiments, Ou et al. (2014) have demonstrated the existence of laminar drag reduction for the flow of water through microchannels influenced by hydrophobic surfaces with well-defined surface roughness. The sensitivity of molecular-scale surface roughness on nanoscale fluid flow was probed in Pit et al.'s experimental work (Pit et al. 2000). From a numerical modeling perspective, Javadpour (2009) proposed an analytical gas flow model that combines the effects of slip flow and Knudsen diffusion in an organic nanopore. Based on this analytical model, the effect of surface adsorption and pressure gradients on the gas-transport behavior in smooth nanopores has been investigated. Some examples show that the contribution of the adsorbed molecules to the gas flux in organic shale represents more than a half of the total mass flux (Kazemi and Takbiri-Borujeni 2016b), while such contribution can be neglected in inorganic pores (Wu and Zhang 2016).  derived an analytical model for the influence of multiple transport mechanisms-that include continuum flow, transition flow and slip flow-on gas transport in shale nanopores. Despite extensive research on gas transport in nanopores, there is as of yet no clear description of the effect of surface roughness on how gas moves through such porous structures.
In addition to methane, also liquid water (connate water and injected water/retained fracturing fluid) is present in shale nanopores (Li et al. 2018;Xu et al. 2017). When it comes to multi-phase transport in porous media, wettability-usually characterized by contact angle and surface tension-is a significant factor controlling fluid distribution and fluid flow (Zhou et al. 2016). There is, however, limited research on the interactions of the various components of multi-phase transport mechanisms in nanoscale pores with roughness.
It is beyond the scope of this paper to review the experimental literature. Spatial resolution of experimental methods varies from sub-micrometer scale to centimeter scale. They include a gas-liquid two-phase flow patterns study in a long microchannel (Triplett et al. 1999) and the centimeter scale core-flooding experiments probing relative permeability in two-phase porous media flow (Perrin et al. 2009). Experimental measurements of fluid flow in a nanochannel are rare. Among these studies (Holt et al. 2006;Lee et al. 2012;Majumder et al. 2005;Whitby et al. 2008), the majority deals with carbon nanotubes and reported flow enhancement and water slippage over their hydrophobic surface. The order of magnitude of flow enhancement-the ratio of the real flux to the value calculated from continuum hydrodynamics-ranges from 10 to 10 5 depending on the channel size (Holt et al. 2006;Majumder et al. 2005;Whitby et al. 2008). Lee et al. (2012) conducted an experimental demonstration of water slippage on hydrophilic alumina surfaces and found that the water permeability is more than double the theoretical prediction using the Hagen-Poiseuille equation for nanochannels with diameters up to 20 nm. For polymeric liquids, De Gennes (2003) predicted a large interfacial slip for entangled polymer melts on ideal surfaces using the Navier boundary condition together with polymer reptation and this Navier-de Gennes model has been confirmed using de-wetting experiments. In most of their studies, the channel surface is smooth and therefore no effects of roughness have been involved.
Limited by the device resolution, experiments are very challenging to realistically replicate the conditions in nanopores of actual formations, let alone to accurately capture fluid flow and interfacial phenomena at the nanoscale. Simulation of pore-scale processes-if properly verified and validated-is an alternative for experimentation. Experimental evidence suggests that continuum simulations are not able to accurately represent the essential physics at the nanoscale (Holt et al. 2006). Given that the molecular nature of matter clearly manifests itself at length scales we are interested in, we have performed molecular dynamics (MD) simulations to study single-and multi-phase transport at the nanoporous scale.
Previous MD-based research includes nanoscopic simulations of fluid flow, droplet wetting and droplet interactions with smooth and rough solid surfaces. The effect of roughness on the behavior of water nanodroplets (Yaghoubi and Foroutan 2018), wetting transition (Khan and Singh 2014;Koishi et al. 2009) and gas flow characteristics (He et al. 2019;Yu et al. 2020) have been reported. Fang et al. (2019) performed MD simulation to elaborate on the mechanisms of extracting hydrocarbon using CO 2 in a nanoscale groove constructed by silicon. Thompson and Troian (1997) did a prominent work in MD. They developed a general boundary condition for liquid flow at solid surface, which allows to relate the degree of slip to the underlying static properties and dynamic interactions of the walls and the fluid. Priezjev and Troian (2006) used three different approaches-Stokes flow calculations, molecular dynamics (MD) simulations and statistical mechanical model-to study the effect of surface roughness on the slip behavior of a Newtonian liquid as a result of planar shear. The authors demonstrated a detailed investigation on the influence of periodic roughness, ranging from molecular to macroscopic dimensions, on the slip length. Stukan et al. (2012) applied MD simulations to probe the mechanism of spontaneous imbibition of a classic aqueous surfactant solution for enhanced oil recovery. Previous studies have shown how permeability enhancement in smooth nanopores depends on contact angle (Wu et al. 2017;Wu et al. 2019a, b). It is an open question if and how such enhancement works out on a rough surface.
The recent work by Ruecker et al. (2020) discussed the relationship between wettability and capillary pressure in a crude oil/brine/rock system; in their work, they simulated (using disjoining pressure curve) and measured (using atomic force microscopy) the water film stability, which controls the wettability alteration, during oil drainage process. The microscale rock surface roughness influences the wettability alternation (contact angle changes) process due to water film collapse and eventually influences on flow regime. In this work, the influence of nanoscale roughness on water droplet sliding and wettability in organic shale pores; this is completed with MD simulation and the organic-rich pore is constructed by carbon atom. The nanoscale roughness is much finer than that in carbonate rock microscale roughness as discussed in the work by Ruecker et al. (2020).
For AFM measurement on actual shale rock, Chalmers et al. (2012) reported that the pore size is mainly < 30 nm, and Zhao et al. (2019) indicated the roughness in shale pore, Ra, reaches 16.9 nm. As the best knowledge of the authors, this is the finest length scale accessed via AFM. However, it is very challenging to measure the shale roughness whose length scale is less than 5 nm. However, molecular dynamics simulations have been used extensively as appropriate tools to understand the fundamentals of rock-fluid interaction in organic shale pore, for example (He et al. 2019;Yu et al. 2020) have studied the roughness (length scale is at 2 nanometers) effect on shale gas transport. Besides, structured surface (with close length scale in our simulation) has also widely been applied in related roughness studies via MD (Khan and Singh 2014;Koishi et al. 2009;Yaghoubi and Foroutan 2018).
In this paper, we report MD simulations that probe single-phase (methane) and twophase (methane/water) transport in shale nanopores as a function of surface roughness. The two-phase systems probe equilibrium and non-equilibrium contact angles on surfaces with finite roughness. It requires to generate extremely larger simulation domain to consider more complicated system that consisting of entities such as organic matter (kerogen), shale grains, shale pores and minerals, and thus it is computational expensive to run MD simulations. For shale rock, we focus on its main deposition of the heavy insoluble organic matter, the kerogen (Loucks et al. 2012), and thus the shale rock system studied in this paper only consists of carbon atoms and is represented by graphite sheets. For kerogen, it has amorphous nano-structure connected by abundant nanometer-sized pores. Given at its complex pore structure at the tiny length scale, it is very challenging to fully reproduce a real kerogen pore. Here in this work, we adopted the recent MD work on kerogen study (He et al. 2019;Kazemi and Takbiri-Borujeni 2016a), in which the graphite slit nanopores with structured roughness have been used as an idealized representation of a kerogen nanopore system.
The remaining sections of this work are organized in the following manner. In the next section, the model system, simulation methodology and the force field models as used in the MD simulations are described. In the Results section, we first present single-phase methane transport in both smooth and rough pores. The effects of surface roughness on the methane volumetric flux and the relative contributions of surface adsorption and bulk transportation to the overall gas flux are quantified. We then move to two-phase methane/ water flow. This includes deformations and sliding speeds of liquid water droplets in a methane environment. The final section reiterates the conclusions of the current work and suggests future research directions.

Nanoslit Pore Simulation Setup
The simulation domains have periodic boundary conditions in all three Cartesian coordinate directions. The systems consist of solid graphite surfaces mimicking a kerogen pore system in shale. For single-phase methane flow simulations, we have two types of pores with smooth and rough surfaces, respectively, see Fig. 1. Pore widths are varied by setting the distance W-the distance between sheet at the top of domain and bottommost sheetto 120 Å, 160 Å and 200 Å, respectively. The rough pores have two additional layers at the top and at the bottom consisting of rectangular roughness elements, see Fig. 1. Periodic arrays of rectangular roughness elements have been used in the literature to mimic nanoscale roughness in MD simulations (Koishi et al. 2009;Lundgren et al. 2007).
For the two-phase methane/water simulations, the graphene sheet at the top of domain has been removed, see

Simulation Methods
The combination of the Lennard-Jones potential and the electrostatic ("Coulombic") potential is used for pairwise interactions: with ij and ij the strength and the length scale of the LJ interaction, respectively, q i and q j the charges of sites i and j, and 0 the dielectric permittivity of the vacuum. The LJ interactions have a cutoff distance of 10 Å. The electrostatic interactions are calculated directly if they are within the same cutoff range and are determined through the particle-particle-particle-mesh (PPPM) method if they are outside the 10 Å cutoff (Darden et al. 1993). Each atom type has been given its own size and strength . The cross-interaction LJ parameters between atoms of different types ( and ) are deduced from the Lorentz-Berthelot mixing rules (Hudson and McCoubrey 1960): For the intramolecular interactions, we have been using the SPC/E model for water (Wu et al. 2006), the OPLS models for CH 4 (Aimoli et al. 2014) and a force field proposed by (Stuart et al. 2000) for graphene. The reasons for choosing these models and details of the complete set of interaction parameters can be found in our previous paper (Yong et al. 2020). Previously, we demonstrated that the thermal motion of the carbon atoms in graphene has negligible impact on the dynamics of the fluid molecules so that also in this paper the carbon atoms in graphene have fixed locations. All MD simulations were performed in a NVT ensemble at a fixed temperature T = 300 K using the open-source molecular dynamic simulation code LAMMPS (Plimpton 1995). The time step, Δt, in this work has been set to 2 fs.
We have applied equilibrium MD (EMD) simulation and non-equilibrium MD (NEMD) simulation in this study. The latter imply applying a constant force on the atoms to generate flow. In the single-phase methane flow simulations, the EMD simulation is first conducted for 1 ns in the absence of a body force to equilibrate the system. This is followed by a 2 ns NEMD simulation. In the latter simulation, a constant external force in the X direction-mimicking a pressure gradient-is added to every atom in each methane molecule. In order for a molecule to accelerate uniformly in the X direction, the total force exerted on the molecule in X direction is distributed over the constituting atoms by the mass of the atoms.
In the methane/water simulation, a 2.4 ns and 10 ns evolution is performed for EMD and NEMD simulations, respectively. In these simulations, the two phases are segregated. A droplet of liquid water on the lower substrate is surrounded by methane. In order to have a uniform force density (and thus pressure gradient), the ratio of the external force exerted on a methane molecule and the force on a water molecule is inversely proportional to the ratio of number density of methane molecules and the number density to water molecules. The trajectories of methane and water are collected every 1 ps for data processing in all simulations.

Single-Phase Methane Flow in Nanopore
The temperature in the single-phase flow simulations involving methane is fixed at T = 300 K. There is a fixed number of methane molecules (4239) so that the q i q j 4 0 r ij volume-averaged density decreases when the width W is increased: from 3.73 mol/l for W = 120 Å to 2.24 mol/l for W = 200 Å. Pressures were determined by means of a virial expression (Yong et al. 2020) after the systems were equilibrated. For W = 120 Å, 160 Å and 200 Å, the pressures were 8.3 ± 0.16 MPa, 7.1 ± 0.13 MPa and 6.0 ± 0.08 MPa, respectively. The pressure, density and temperature remain the same (within 1.5%) in the rough pores as compared to the smooth ones.
In shale gas operations, the superficial gas flow velocities are of the order of 10 −2 -10 −1 m/s (Kenomore et al. 2018). The random thermal velocities at 300 K are of the order of 10 2 -10 3 m/s. Such disparity of organized (flow) velocities versus random velocities poses a problem for extracting average velocity profiles from MD simulations. If we denote the standard deviations related to the random fluctuations of the velocity as , then the standard deviation of the average velocity based on N statistically independent velocity samples is av = √ to approximately 10 which then means that N can be reduced by a factor 10 6 . With u av ≈ 10 the random velocities are still significantly larger than the average velocities which is important for a correct representation of the essential physics. The large relative velocities are achieved by setting a very large pressure gradient of 0.17 × 10 8 MPa/m. We will keep this in mind when interpreting the results of our simulations. Figure 3 shows the methane density distributions over the width of the pore for NEMD simulations. The EMD simulations details are omitted because the EMD density distribution shows no difference (within 1%) with NEMD simulations. This indicates even with a very high pressure gradient density profiles are not affected by the pressure gradient, and thus pressure keeps unchanged in EMD and NEMD simulations (pressure is determined by means of a virial expression from the equilibrium bulk region of which the density keep unchanged in EMD and NEMD simulations).
In the center region, the methane density is uniform (He et al. 2019;Yu et al. 2020) and the same (within 0.8%) in smooth and rough pores. The density peaks at the substrate indicate the formation of an adsorbed layer of methane. The adsorption peaks at the rough substrates are lower than at the smooth substrates. These phenomena agree with previous work (Yu et al. 2020) on single-phase methane flow in rough nanochannels at T = 298 K and P = 10 MPa. In Yu et al. (2020), it is shown that the adsorption peaks drop constantly as the roughness increases in an organic shale pore with fixed pore width of 5 nm, while the bulk density is invariable for different roughness. The adsorption density (simplified as the maximum density value of the adsorption peak) is about 10 times higher than that in the bulk region for ideally smooth pores.  Fig. 3 Left: methane density distribution along Z direction in a smooth pore in NEMD simulations; right: density distribution in a pore with roughness fraction f = 0.25 in NEMD simulations. Three widths W as indicated. Gray regions are defined as the adsorption layers. The pressure gradient equals 0.17 × 10 8 MPa/m; T = 300 K; average densities are 3.73 mol/l for W = 120 Å, 2.80 mol/l for W = 160 Å and 2.24 mol/l for W = 200 Å. The difference in density between smooth and rough pores is less than 0.5%. Uncertainties are smaller than the symbol size Fig. 4 Average velocity in X direction as a function of Z for three width W as indicated. Solid curves are least-square fits (according to Eq. 4). Left: smooth substrate; right rough substrate with f = 0.25. Pressure gradient, temperature and densities are the same as in Fig. 3. On data points without an error bar, the uncertainty is less than the symbol size Takbiri-Borujeni 2016b) it has been stated that between 20 and 40% of an overall gas flux is due to adsorbed layers, which is similar to our results. For rough pores we find Φ = 0.06, 0.039 and 0.033 for W = 120, 160 and 200 Å, respectively, which indicates that mobilizing adsorbed methane in rough pores is much harder than in smooth pores.

Volumetric Flux of Single-Phase Methane Flow in Nanopore
We now analyze the velocity profiles in Fig. 4 in terms of the modified Hagen-Poiseuille (mHP) equation Firoozabadi 2015, 2016;Wang et al. 2016) where the modification accounts for wall-slip effects: with Γ the body force (force per unit volume) representing the pressure gradient driving the flow, the viscosity of the fluid, and L s the slip length. The superficial velocity (volumetric flow rate divided by cross-sectional area) follows from integrating Eq. 3 over the width of the pore where we assume that the viscosity is constant, i.e., not a function of z: The simulated velocity distributions (Fig. 4) have been fitted with parabolic curves: v(z) = az 2 + c with a and c the two least-squares fitting parameters. The slip length L s is related to the fitting parameters according to L s = − c Wa − W 4 . The superficial velocity can be expressed in terms of the fitting parameters as J f = 1 12 aW 2 + c. From fitting, the slip lengths for the smooth pores are L s = 57.6 Å, 80.2 Å and 88.9 Å for W = 120 Å, 160 Å and 200 Å, respectively. The slip lengths for the rough pores are L s = 3.15 Å, 3.69 Å and 9.32 Å which implies that slip can be neglected in rough pores.
Equation 4 with L s = 0 has been used to estimate the superficial velocity in rough pores. This requires the dynamic viscosity of methane at a temperature of 300 K and at the bulk density (i.e., the plateau densities in Fig. 3). For this, we used the NIST web-book (Shen et al. 2011). Viscosity values are = 13.12, 12.78, 12.47 μPa·s for the situations with W = 120 Å, 160 Å and 200 Å, respectively. Figure 5 shows that for the rough pores the estimates of the superficial velocity based on Eq. 4 agree well with the ones based on the curve fit. This good agreement, in hindsight, justifies the use of a uniform viscosity that can be estimated at bulk conditions.
For the smooth surfaces, we want to perform a similar exercise: comparing theoretical superficial velocities (Eq. 4) with MD simulation results. This would require a priori knowledge of the slip length which we do not have. The theoretical results (based on Eq. 4) in Fig. 5 that relate to smooth pores therefore use, in addition to the NIST viscosity data, the slip length estimates from the curve fits. Although this makes the theoretical results not fully predictive, we see good agreement between theory and MD. It implies that the contribution of the superficial velocity not related to slip effects is well captured by the MD simulations. Figure 6 shows the configuration as used for water/methane simulations. A liquid water droplet is deposited on a substrate (in the case of Fig. 6, this is a substrate containing roughness elements) and is immersed in methane. As described above, the width of the pore has been fixed to W = 120 Å. There are 2015 water molecules and 4239 methane

Droplet Sliding Over Nanopore in Methane Background
molecules in a simulation. The system is at a temperature of 300 K. The pressure under equilibrium conditions the pressure has been determined as 8.3 ± 0.16 MPa. After equilibration we apply a force in the x direction on each atom in each molecule with the force per atom type weighted in such a way that (1) each molecule is accelerated uniformly and (2) the body force is uniform throughout the volume between the substrates. The force distributions mimics a pressure gradient p x = 0.18 × 10 8 MPa/m. This represents a pressure difference of L x p x ≈ 0.2 MPa over the length of the domain which is much smaller than the overall absolute pressure of approximately 8 MPa.  In what follows, we will mainly study the effect of the surface roughness-as quantified by the surface fraction f-on the deformation and motion of the water droplet.
Droplet motion is monitored by plotting its footprint as a function of time, see Fig. 7. For this we consider contours of the short-time-average concentration of water droplets in a xy plane closely (4 Å) above the rough substrate. Time averaging is done over 0.2 ns to reduce thermal noise and have a footprint with a clear outline. Figure 7 shows a gradual displacement of the footprint in the positive x direction as a result of the pressure gradient that drives the flow of methane and water. The vertical (xz) cross sections in Fig. 7 show a sliding water drop that hardly deforms under the shear stresses provided by the methane flowing over it.
Plotting the center position of the footprint x c as a function of time, as we do in Fig. 8, allows us to determine the sliding speed of the water drop. To a good approximation linear behavior of x c versus t is observed so that, from the moment it starts sliding, the drop has a  Fig. 9 that shows droplet sliding speed as a function of f.
The pressure gradient on droplet sliding was investigated further. In Fig. 10, results for sliding speed as a function of pressure gradient for a number of roughness values f are presented. A smooth surface ( f = 1 ) is unique in the sense that sliding speed is proportional to the pressure gradients. This allows to make predictions for droplet sliding velocities at the much lower pressure gradients encountered in shale gas operations, where at 0.18 × 10 8 MPa/m, the sliding speed over a smooth surface is approximately 20 m/s, at a more common pressure gradient (Dai et al. 2017) of 0.14 MPa/m (we assume a single gas well drainage radius is 100 m) the sliding speed would be approximately 2 × 10 −7 m/s. Over rough surfaces, the relationship between sliding speed and pressure gradient is still linear, see Fig. 10. Rough surfaces, however, have a critical pressure gradient below which the drop gets stuck, i.e., does not slide. This even is the case for an almost smooth surface (with f = 0.97). For this case, the critical pressure gradient is approximately 0.02 × 10 8 MPa/m. Scaling this back to more practical pressure gradients as encountered in shale (0.14 MPa/m), this implies that already a minute level of surface roughness The results may help to explain water will become immobile and be retained in the matrix permanently once it is imbibed into the low-permeability micropores of the shale matrix (He 2011;Li et al. 2019). Average velocity profiles over the pore width of the methane/water system are given in Fig. 11 for four different levels of roughness. The trends in the methane profiles can be interpreted along the same lines as the single-phase methane results in the previous subsection: at the substrate with the highest roughness (f = 0.76), the average methane velocity at the substrate is almost zero (almost no-slip). Smoother substrates allow for more and more slip. The upper substrate of the pore at z = W = 120 A is smooth (see Fig. 2) so that methane has a considerable slip velocity there. The water velocity profiles hardly show a velocity gradient which implies that the fluid deformation in the drop is very small; we will come back to this when discussing dynamic contact angles. In all cases, the velocity of the water drop is lower than the slip velocity of methane at the lower solid surface.

Static and Dynamic Contact Angles at a Rough Substrate
We now discuss the shape of static and sliding droplets in terms of contact angles: the static contact angle for equilibrated drops and receding and advancing contact angles for drops feeling a pressure gradient. Contact angles have been determined by fitting circular arcs to the time-averaged outline of the xz midplane of a droplet, see Fig. 12. The angle with which the arc intersects with the top of the substrate is identified as the contact angle. For a static drop, we use one circular arc over the entire cross section; for sliding drop, we use two arcs, one at the advancing side and one at the receding side of the drop.
In terms of methods for contact angle calculation, for a given fluid/droplet configuration at pore scale, three approaches can be used to determine contact angle. The first is geometric-based contact angle estimated directly from the fluid-rock image (e.g., Andrew et al. 2014;Yang and Zhou 2020). The second is to use the energy balance to calculate the effective contact angle for a fluid configuration change that could be obtained from high-resolution imaging (Blunt et al. 2019). The concept of topology and Fig. 10 Left panel: the relationship between water droplet sliding speeds and pressure gradient for various f. An enlarged view of details within red square is shown in right panel. The solid lines are linear functions that best fit the nonzero data points. Dash lines denotes a transit region where water droplet starts sliding from a "stuck" state at the given conditions. P = 8.3 MPa and T = 300 K. On data points without error bars, the uncertainty is less than the symbol size integral geometry has been used to derive an effective macroscopic contact angle to study the relationship between intrinsic, advancing and receding contact angle based on fluid configuration imaged at pore scale (Sun et al. 2020a, b). In this work, the droplet configuration is determined based on the molecular fraction concentration that simulated based on MD simulation, and then the contact angle is estimated using a curve fitting approach. Other alternative approaches to estimate the effective contact angle  1 for panel a.,b.,c.,d, respectively. The pressure gradient is 0.18 × 10 8 MPa/m, P = 8.3 MPa and T = 300 K Fig. 12 Fitting procedure for determining the receding and advancing contact angles in the example of f = 0.4. a Snapshot of symmetry plane XZ cross section in the simulation system. b Time-averaged water droplet density contour (the averaging time is 5 ns). c Points with concentration 0.015/Å 3 and (least squares) fitted arcs for determining contact angles. The filled symbols in (c) do not take part in the fitting process. The horizontal reference line in (c) is the top of the substrate based on the MD simulated droplet configuration includes the above mentioned three approaches but they require other relevant information.
In Fig. 13, static and dynamic contact angles are reported as a function of roughness parameter f at a pressure gradient of 0.18 × 10 8 MPa/m. The incipient sliding motion at this pressure gradient is at approximately f = 0.7. Drop deformation is limited given that the difference between advancing and receding contact angle is less than 6°. This difference is hardly larger than the uncertainty with which contact angles have been determined through curve fitting. These uncertainties are likely the reason why sometimes the static contact angle does not lie between advancing and receding angle as would be expected. Limited drop deformation as reflected by limited contact angle hysteresis is consistent with the average velocity profiles in Fig. 11 that hardly show velocity gradient for the water phase. It is interesting to note that the contact angles as a function of f go through a maximum and sliding starts directly beyond the f associated with maximum contact angles.
We added more simulations to study the effect of roughness, such as vertical aspect ratio and patterns, on the dynamic contact angle (deformation of drop) and its pinning effect. As shown in Figure S1a, f = 0.25, a more shaper pillar surface is created by adding 5 layers graphite sheets pillar height. The interlayer distance keeps at 3.35 Å, and a water droplet is placed on the center of the substrate in the presence of methane. We then performed EMD and NEMD simulations with the same conduction as in the previous system where a onelayer pillar substrate is used. The quantitative comparison of the static and dynamic contact angles for one-layer pillar substrate and five-layer pillar substrate is shown in Figure S1b. It is shown the vertical aspect ratio roughness influences the water droplet contact angles, while there is still a small contact hysteresis effect. We further tested another configuration of surface roughness that has a "ridges" pattern with the same f, f = 0.25. The schematic of ridges substrate is shown in Figure S2a. A same simulation conduction has been conducted as in one-layer pillar substrate simulation, and the contact angles results are summarized in Figure S2b. We also noticed a limited contact angle hysteresis behavior for both substrates in this case. The simulated results indicate that water contact angle is dependent on roughness patterns and roughness vertical/lateral area fraction, while its hysteresis effect over graphite substrate is very limited and independent on the surface topography. But we do agree that the wetting behavior may change if more complicated surface roughness configuration is used, and we are currently performing more simulation work on more complicated surface that generated based on, e.g., fractal theory, which we hope to consider for another publication in the future.

Summary & Conclusions
In this work, we have presented a molecular-scale study on single-phase methane flow and two-phase methane/water flow in organic shale pores with surface roughness taken into account by means of full-atoms MD simulations. We investigated the effects of roughness on interactions between methane, water droplets and organic shale surface under a constant pressure gradient driven flow. It mainly involves studies on roughness influence on methane adsorption, methane and water/methane flow transport characteristic and water wettability. The conclusions of our simulation are: • Surface roughness significantly influences methane transportation and reduces permeability at organic nanopore. • Methane density distribution is not affected even with a very high pressure gradient. In agreement with the literature, the relative contribution of adsorption layers caused by slip flow to the overall gas flux is in a range of 19-28%. • The modified HP equation is used to analyze the methane velocity profiles recovered from MD simulations. Our results justify the use of a uniform viscosity that can be estimated at bulk conditions. The single-phase methane flow with a slip length in organic shale pores can be described by a combination of the modified HP question and MD results. • Water droplets exhibit a proportional relationship between its sliding speeds and pressure gradient over ideally smooth organic surface, in a practical shale pore where surface roughness immobilizes liquid water. • Our results indicate the static water contact angle is dependent on the surface topography, while there is always a limited contact angles hysteresis effect for water droplet over both smooth and rough organic substrate, and the sliding water droplets hardly deform over organic surface in the presence of methane.