Pore-Scale Simulation of Confined Phase Behavior with Pore Size Distribution and Its Effects on Shale Oil Production

Confined phase behavior plays a critical role in predicting production from shale reservoirs. In this work, a pseudo-potential lattice Boltzmann method is applied to directly model the phase equilibrium of fluids in nanopores. First, vapor-liquid equilibrium is simulated by capturing the sudden jump on simulated adsorption isotherms in a capillary tube. In addition, effect of pore size distribution on phase equilibrium is evaluated by using a bundle of capillary tubes of various sizes. Simulated coexistence curves indicate that an effective pore size can be used to account for the effects of pore size distribution on confined phase behavior. With simulated coexistence curves from pore-scale simulation, a modified equation of state is built and applied to model the thermodynamic phase diagram of shale oil. Shifted critical properties and suppressed bubble points are observed when effects of confinement is considered. The compositional simulation shows that both predicted oil and gas production will be higher if the modified equation of state is implemented. Results are compared with those using methods of capillary pressure and critical shift.


Introduction
Shale oil production experienced a rapid increase in recent years and has been a major contributor to oil consumption in North America [1][2][3]. However, modeling and prediction of production from shale reservoirs are still challenging. One reason is that the state and motion of hydrocarbon molecules in nanopores is not thoroughly studied [4]. It is known that pore size in shale ranges from several to hundreds of nanometers [5]. The strong interactions between fluid molecules and rock surface lead to complex flow, transport, and phase behavior at the microscale. Also, the unique physics of fluids in nanopores has nontrival effects on predicted production and ultimate recovery from shale reservoirs [6][7][8]. Thus, it is essential to understand the mechanisms of fluids in confined space and include their effects when modeling shale reservoirs.
Owing to the abunt nanopores in the shale matrix, the phase equilibrium of oil and gas deviates significantly from that at the bulk condition [4]. Conventional equation of state (EOS) is not suitable to model phase change of fluids for shale reservoir [4]. The study of confined phase behavior is thus important and popular in recent years. Deviated phase equilibrium has been found experimentally. Qiao et al. [9] and Russo et al. [10] reported decreased saturation pressures and critical temperatures in nanopores from measured adsorption isotherms. Morishige et al. [11] measured critical temperatures of several gases in different pore sizes of nanopores and found they are all below the critical temperature at the bulk condition. By using the technology of nanofluidics, Wang et al. [12] and other researchers [13][14][15][16] observed that the critical temperatures and bubble point pressures of hydrocarbons are lower in nanopores compared to those in large pores.
Besides, confined phase behavior was also studied numerically by using molecular dynamic (MD), Monte Carlo (MC) method, and Lattice Boltzmann (LB) method. Sedghi and Piri [17] studied capillary condensation in nanopores using MD and claimed that there is a critical pore size below which capillary condensation cannot occur. Ambrose et al. [18] modeled the density distribution of methane in slit-like pores. They reported that the adsorption thickness is around 0.38 nm which is close to the molecular diameter of methane. Compared to MD, MC is more widely used to study confined phase behavior in nanopores. For instance, Peterson and Gubbins [19] studied the phase change of a Lennard-Jones fluid in nanopores usng MC and found the critical temperature decreases in nanopores compared to bulk. Similar work was also been done by Hamada et al. [20]. Singh et al. [21,22] reported that the critical temperature decreases linearly with the inverse of slit width and critical density fluctuates with a decrease in the slit pore. Besides, Vishnyakov et al. [23] conducted numercial study on the phase behavior of methane in slit pores. They found that a lower critical temperature and higher critical density of methane in slit pores and the shift of critical point depends on pore size and strength of fluid-solid force. More recently, the Lattice Boltzmann method has been applied to study phase behavior under effects of confinement. Huang et al. [24] evaluated the effect of large capillary pressure on the vapor-liquid equilibrium of several gases. It is found that the phase equilibrium of nanoscale droplets follows Laplace law and Kelvin equation. In another work of Huang et al. [25], the pseudo-potential LB model is applied to study the confined phase behavior of methane in a slit pore and a complex medium. A suppressed critical temperature and increased critical density are observed which is consistent with molecular simulation. Also, density function theory (DFT) is an efficient tool for theoretical studies of confined phase behavior. Findings from studies using DFT [26][27][28] were often consistent with molecular simulations.
Numerical simulation of confined fluid can produce a coexistence curve of liquid and vapor in nanopores. However, it is not straightforward to directly use it in thermodynamic modeling. Thus, efforts have been made to modify the EOS. To model altered phase behavior in nanopores, some researchers choose to couple Peng-Robinson (PR)-EOS with capillary pressure [6,[29][30][31]. It is found that high capillary pressure will suppress bubble point pressure and increase upper dew point pressure. Another popular way for modifying EOS is using the critical shift method [32][33][34]. The phase diagram is found shrink in nanopores and bubble point and dew point are both suppressed with a shifted critical point. Other researchers tried to propose a new form of EOS by considering molecule-wall interaction. In the work of Luo et al. [35], adsorbed and bulk fluids in nanopores were treated separately and a pore-size-dependent EOS was proposed. Travalloni et al. [36] extended PR-EOS by analytically analyzing interactions between fluid molecules and walls. Wang and Aryana [37] extended PR-EOS guided by results from molecular simulation. The shift of critical point was captured and a temperature-dependent parameter was introduced to account for the effects of capillary pressure. With the EOS proposed for confined fluids, it can be incorporated into the flash calculation and compositional simulation to model oil and gas production from shale. It is reported that the predicted oil production is higher and predicted gas production is lower if the effect of oil-gas capillary pressure is considered in reservoir simulation [6,30,38]. Implementing the critical shift method usually leads to both higher oil and gas production [34,39,40]. Other forms of EOS can also be applied in compositional simulation and not much work has been conducted.
In this work, we aim to present the process of thermodynamic modeling of confined phase behavior using pore-scale simulation and applying the knowledge learned into the reservoir-scale simulation. First, the pseudo-potential LB model is validated in modeling phase equilibrium in a single pore. Then, the effect of pore size distribution on confined phase behavior is studied using a bundle of capillary tubes. An effective pore size is used to represent distributed pores in the real rock sample and a modified EOS is proposed. Finally, the modified EOS is applied to model phase diagram in shale matrix and predict oil and gas production. Results are compared with those using methods of capillary pressure and critical shift.

Methodologies
In this part, we briefly introduce the pore-scale and large-scale methods applied in this work. Specifically, a multiphase lattice Boltzmann method is first presented for simulation of confined phase behavior. Then, a modified flash calculation and compositional simulation model will be introduced.

Pseudopotential Lattice Boltzmann Method
In this work, a lattice Boltzmann (LB) method is applied to model phase behavior of confined fluids. As a mesoscopic method, LB method simulates velocity distribution of fluid molecules. The general Equation of the LB method is given by: where n i is the velocity distribution function of fluid molecules at position x and time t, e i is the discrete velocity, subscript i denotes the index of lattice velocities. ∆t is the time step which is usually set as one in LB simulation. F i (x, t) is the external forcing term. Ω i represents the collision operator that accounts for the effects of collision among fluid molecules. In this study, the multi-relaxation time (MRT) collision operator is applied because of its superiority of stability at low temperatures [24]. Acoording to [41], the MRT collision process can be generally expressed by: where M is the transformation matrix and S is the diagonal collision matrix in the moment space and n eq j is the distribution function at an equilibrium state. To model phase change in a more meaningful scale, a three-dimensional, nineteen velocity model (D3Q19) is applied. In this model, discrete velocities e i are given as: Transformation matrix M is defined as the work of [41]. The diagonal collision matrix S ≡ diag 0, s e , s ε , 0, s q , 0, s q , 0, s q , s ν , s π , s ν , s π , s ν , s ν , s ν , s t , s t , s t . Relaxation parameters are chosen as: s q = s t = 8(2 − 1/τ)/(8 − 1/τ), s e = s ε = s ν = s π = 1/τ. τ is the relaxation time. In Equation (2), the distribution function at equilibrium state is written as: where c s = √ 1/3 is speed of sound in LB, ω i is the weighting parameter that is ω 1 = 1/3, ω 2−7 = 1/18, ω 8−19 = 1/36. ρ and u are the macroscopic density and velocity that can be obtained from the distribution function by ρ = ∑ i n i , u = (∑ i n i e i )/ρ.
To model the phase behavior of non-ideal gas, a pseudo-potential LB model is applied. In this model, the pseudo-potential is introduced to account for the intermolecular interaction. The intermolecular interactive force is given by [42]: where ψ(x) is the pseudopotential, c 0 is a constant depending on a lattice structure and c 0 = 6 when using the D3Q19 model. G f controls the strength of the interaction. In the pseudopotential model, the pressure of non-ideal gas is given by: The form of ψ(x) determines Equation of state. In this work, the Peng-Robinson Equation of state is used: where T is temperature, R is gas constant, a = 0.45724R 2 T 2 c /P c , b = 0.0778RT c /P c , T c and P c are the critical temperature and critical pressure respectively. Similar to the study in [43], we have set a = 2/49, b = 2/21, and R =1 in our simulation. Combing Equation (6) and Equation (7), the pseudopotential is in the following form: Also, the adhesive force between solid and fluid molecules is incorporated by creating an analog to intermolecular interaction [44]: where s(x + c i ) is a switching function: s = 1 or 0 when x + c i falls into the solid phase or fluid phase. G w is the wetting parameter that can be used to tune the strength of adhesive force and contact angles. To incorporate intermolecular and adhesive forces, the exact difference method is applied [45]. The current LB model is limited for single component and model for gas mixtures need to be developed in the future.

Thermodynamic Modeling of Vapor-Liquid Equilibrium
The vapor-liquid equilibrium (VLE) method is employed to model thermodynamic phase behavior. At thermodynamic equilibrium, the fugacity of components in liquid and vapor phases are equal: where f l i and f v i are the fugacity of component i in liquid and vapor phase respectively, P l and P v are the pressure of liquid and vapor phase respectively, x i and y i are the mole fraction of component i in the liquid and vapor phase respectively. The fugacity of a single component in a gas mixture can be given as: where φ l i and φ v i are fugacity coefficient of component i in liquid and vapor phase respectively. The component equilibrium ratio K i (K i = y i /x i ) is updated by using the successive substitution method: An equilibrium state is reached until: Energies 2021, 14, 1315 where ε is the tolerance that is smaller than 10 −14 .
When the pore size reaches nanometers, pressure difference between liquid and vapor phases could be very large. According to the Young-Laplace Equation, pressure difference or capillary pressure can be calculated by: where r is pore radius, θ is contact angle and σ is the interfacial tension given by: where χ i is the parachor of component i, ρ l and ρ v are the molar density of liquid and vapor phase respectively. To account for the effect of capillary pressure on VLE, capillary pressure is updated in parallel with fugacity by: util P old,n c − P new,n c < ε c . The threshold for iteration of capillary pressre should be small enough that the ultimate VLE will not be affected. In this work, ε c is set as 0.5 psi. Apart from effects of capillary pressure, it is found that the critical point of confined fluids are also altered by the strong interaction between solid surface and fluid molecules. One commonly used method is from Zarragoicoechea and Kuz [32] that is given by: where T c and P c are critical temperature and pressure at bulk condition, and T cp and P cp are the critical temperature and pressure in nanopores, respectively. σ ij = 0.244 3 √ T c /P c .

Simulation of Production from Shale Reservoir
In this work, the compositional model is used to predict oil and gas production from shale reservoir. An in-house compositional simulator has been developed considering multiple physics [46]. The mass balance Equation can be generally written as [31]: where φ is porosity, V j is the bulk volume of cell j, S o and S g are the oil and gas saturation, λ ro and λ rg are the relative mobility of oil and gas, Φ o and Φ g are the potential of oil and gas, q p o and q p g are the production rates of oil and gas respectively. In a compositional simulation, fluid property including density and viscosity will be updated by VLE every time step. T is the transmissibility of connection.
To model fractures in shale reservoirs, the embedded discrete fracture model (EDFM) [47,48] is used. The mass transfer between hydraulic/natural fractures and the matrix is calculated by: where A nnc is the contact area, k nnc is the harmonic average of fracture permeability and matrix permeability, d nnc is the distance of connection can be determined by: where dist(x, y, z) is the distance from location (x, y, z) to the fracture. The detailed model can be referred to our previous work [31].

Results and Discussion
At the beginning, we present pore-scale simulation of phase equilibrium in a capillary tube with a uniform pore size. Then, effect of pore size distribution on confined phase behavior is evaluated by using a bundle of capillary tubes with different sizes. With simulated coexistence curves, a modified EOS is established. Then, the new EOS is implemented into thermodynamic models for predicting the phase change of hydrocarbon mixtures. Finally, oil and gas production from a fractured shale reservoir is conducted using the modifed EOS, and results are compared with those using the methods of shifting critical points and implementing capillary pressure.

Phase Equilibrium in a Single Pore
In this part, we first evaluate the phase separation of methane at the bulk condition. Simulation is conducted in a 100 3 domain. The periodic boundary condition is applied at x, y, z-direction so that methane can be regarded as bulk gas. Initially, the computational domain is filled with methane with equal density everywhere and the temperature is set below the critical temperature. By a small perturbation applied, phase separation of methane can be achieved due to the intermolecular interactions. In this way, liquid and vapor densities of methane at the temperature are obtained from simulation. Repeating the same process in different temperatures, we can construct the coexistence curve using the LB method. In Figure 1, the coexistence curve from the simulation is compared with the theoretical one obtained using Maxwell construction on PR-EOS. The good match indicates the reliability of this method for phase behavior study.
To model fractures in shale reservoirs, the embedded discrete fracture model (EDFM) [47,48] is used. The mass transfer between hydraulic/natural fractures and the matrix is calculated by: where is the contact area, is the harmonic average of fracture permeability and matrix permeability, is the distance of connection can be determined by: where ( , , ) is the distance from location (x, y, z) to the fracture. The detailed model can be referred to our previous work [31].

Results and Discussion
At the beginning, we present pore-scale simulation of phase equilibrium in a capillary tube with a uniform pore size. Then, effect of pore size distribution on confined phase behavior is evaluated by using a bundle of capillary tubes with different sizes. With simulated coexistence curves, a modified EOS is established. Then, the new EOS is implemented into thermodynamic models for predicting the phase change of hydrocarbon mixtures. Finally, oil and gas production from a fractured shale reservoir is conducted using the modifed EOS, and results are compared with those using the methods of shifting critical points and implementing capillary pressure.

Phase Equilibrium in a Single Pore
In this part, we first evaluate the phase separation of methane at the bulk condition. Simulation is conducted in a 100 3 domain. The periodic boundary condition is applied at x, y, z-direction so that methane can be regarded as bulk gas. Initially, the computational domain is filled with methane with equal density everywhere and the temperature is set below the critical temperature. By a small perturbation applied, phase separation of methane can be achieved due to the intermolecular interactions. In this way, liquid and vapor densities of methane at the temperature are obtained from simulation. Repeating the same process in different temperatures, we can construct the coexistence curve using the LB method. In Figure 1, the coexistence curve from the simulation is compared with the theoretical one obtained using Maxwell construction on PR-EOS. The good match indicates the reliability of this method for phase behavior study.  After simulating phase equilibrium at bulk condition, phase change in a slit pore is studied. The simulation setup is shown in Figure 2. A three-dimensional capillary tube with a rectangular cross-section is used. The solid boundary condition is applied in the y and z-direction. In x-direction, constant pressure is held at the entrances of the tube. Similar to the case of bulk gas, the tube is initially filled with methane with equal Energies 2021, 14, 1315 7 of 17 density everywhere and the temperature is below the critical temperature. At inlet and outlet boundaries, the pressure/density of the gas is fixed. By gradually increasing the density of gas at the boundary, the change of average density inside the tube is recorded. Figure 3 presents examples of average densities as functions of densities at the boundary at three temperatures. The length of the side of the tube is 10 nm. Clearly, the average density undergoes a sudden jump, from a vapor-like density to liquid-like density. This sudden change of average density inside the tube indicates the occurrence of capillary condensation. The average densities right before and immediately after the sudden jump are regarded as the vapor and liquid densities at a specific temperature. By following the same procedure using different temperatures, we can construct the coexistence curve in the capillary tube. After simulating phase equilibrium at bulk condition, phase change in a slit pore is studied. The simulation setup is shown in Figure 2. A three-dimensional capillary tube with a rectangular cross-section is used. The solid boundary condition is applied in the y and zdirection. In x-direction, constant pressure is held at the entrances of the tube. Similar to the case of bulk gas, the tube is initially filled with methane with equal density everywhere and the temperature is below the critical temperature. At inlet and outlet boundaries, the pressure/density of the gas is fixed. By gradually increasing the density of gas at the boundary, the change of average density inside the tube is recorded. Figure 3 presents examples of average densities as functions of densities at the boundary at three temperatures. The length of the side of the tube is 10 nm. Clearly, the average density undergoes a sudden jump, from a vapor-like density to liquid-like density. This sudden change of average density inside the tube indicates the occurrence of capillary condensation. The average densities right before and immediately after the sudden jump are regarded as the vapor and liquid densities at a specific temperature. By following the same procedure using different temperatures, we can construct the coexistence curve in the capillary tube.  In Figure 4, the coexistence curves of methane in different sizes of nanotubes are compared with that at the bulk condition. The liquid-gas contact angle is assumed to be 30 o . A significant shift of the coexistence curve is observed in the nanotube, especially when the size of the tube is small. Compared to bulk condition, liquid density is decreased and vapor density is increased. Also, a clear shift of critical point is observed in nanopores. Specifically, the critical temperature decreases, and critical density increases with After simulating phase equilibrium at bulk condition, phase change in a slit pore is studied. The simulation setup is shown in Figure 2. A three-dimensional capillary tube with a rectangular cross-section is used. The solid boundary condition is applied in the y and zdirection. In x-direction, constant pressure is held at the entrances of the tube. Similar to the case of bulk gas, the tube is initially filled with methane with equal density everywhere and the temperature is below the critical temperature. At inlet and outlet boundaries, the pressure/density of the gas is fixed. By gradually increasing the density of gas at the boundary, the change of average density inside the tube is recorded. Figure 3 presents examples of average densities as functions of densities at the boundary at three temperatures. The length of the side of the tube is 10 nm. Clearly, the average density undergoes a sudden jump, from a vapor-like density to liquid-like density. This sudden change of average density inside the tube indicates the occurrence of capillary condensation. The average densities right before and immediately after the sudden jump are regarded as the vapor and liquid densities at a specific temperature. By following the same procedure using different temperatures, we can construct the coexistence curve in the capillary tube.  In Figure 4, the coexistence curves of methane in different sizes of nanotubes are compared with that at the bulk condition. The liquid-gas contact angle is assumed to be 30 o . A significant shift of the coexistence curve is observed in the nanotube, especially when the size of the tube is small. Compared to bulk condition, liquid density is decreased and vapor density is increased. Also, a clear shift of critical point is observed in nanopores. Specifically, the critical temperature decreases, and critical density increases with In Figure 4, the coexistence curves of methane in different sizes of nanotubes are compared with that at the bulk condition. The liquid-gas contact angle is assumed to be 30 • . A significant shift of the coexistence curve is observed in the nanotube, especially when the size of the tube is small. Compared to bulk condition, liquid density is decreased and vapor density is increased. Also, a clear shift of critical point is observed in nanopores. Specifically, the critical temperature decreases, and critical density increases with decreasing pore sizes. The shifted critical properties are consistent with results using molecular simulation [21,22] and density function theory [26,27]. decreasing pore sizes. The shifted critical properties are consistent with results using molecular simulation [21,22] and density function theory [26,27].

Phase Equilibrium with the Pore Size Distribution
In Section 3.1, we have shown the calculation of phase equilibrium in a single capillary tube. However, it is known that pore size in shale always yields a distribution [4]. Thus, it is important to study phase behavior in a medium with distributed pore sizes. As shown in Figure 5, normalized pore size data of a shale sample from Eagle Ford is used [49]. The measured data is from Hg intrusion [49]. Four representative pore sizes are used to represent the distributed pores. Note that the choice is not unique. The chosen representative pore size should have relatively high pore volume and significantly deviated phase behavior. As illustrated in Figure 5, pore sizes of 4 nm, 6 nm, 8 nm, and 14 nm are chosen to cover the pores with sizes of 3-5 nm, 5-7 nm, 7-9 nm, 9-30 nm respectively. The pore volume of each size of pores equals the volume of pores covered. Following the work of Jin et al. [50], pores larger than 30 nm are ignored where the confinement effect on phase behavior can be ignored.  [50]) and discrete pore size distribution used in LB simulation. PSD represents pore size distribution.

Phase Equilibrium with the Pore Size Distribution
In Section 3.1, we have shown the calculation of phase equilibrium in a single capillary tube. However, it is known that pore size in shale always yields a distribution [4]. Thus, it is important to study phase behavior in a medium with distributed pore sizes. As shown in Figure 5, normalized pore size data of a shale sample from Eagle Ford is used [49]. The measured data is from Hg intrusion [49]. Four representative pore sizes are used to represent the distributed pores. Note that the choice is not unique. The chosen representative pore size should have relatively high pore volume and significantly deviated phase behavior. As illustrated in Figure 5, pore sizes of 4 nm, 6 nm, 8 nm, and 14 nm are chosen to cover the pores with sizes of 3-5 nm, 5-7 nm, 7-9 nm, 9-30 nm respectively. The pore volume of each size of pores equals the volume of pores covered. Following the work of Jin et al. [50], pores larger than 30 nm are ignored where the confinement effect on phase behavior can be ignored. decreasing pore sizes. The shifted critical properties are consistent with results using molecular simulation [21,22] and density function theory [26,27].

Phase Equilibrium with the Pore Size Distribution
In Section 3.1, we have shown the calculation of phase equilibrium in a single capillary tube. However, it is known that pore size in shale always yields a distribution [4]. Thus, it is important to study phase behavior in a medium with distributed pore sizes. As shown in Figure 5, normalized pore size data of a shale sample from Eagle Ford is used [49]. The measured data is from Hg intrusion [49]. Four representative pore sizes are used to represent the distributed pores. Note that the choice is not unique. The chosen representative pore size should have relatively high pore volume and significantly deviated phase behavior. As illustrated in Figure 5, pore sizes of 4 nm, 6 nm, 8 nm, and 14 nm are chosen to cover the pores with sizes of 3-5 nm, 5-7 nm, 7-9 nm, 9-30 nm respectively. The pore volume of each size of pores equals the volume of pores covered. Following the work of Jin et al. [50], pores larger than 30 nm are ignored where the confinement effect on phase behavior can be ignored. Figure 5. Normalized pore size distribution of shale rock obtained from experiments (data obtained from [50]) and discrete pore size distribution used in LB simulation. PSD represents pore size distribution.  [50]) and discrete pore size distribution used in LB simulation. PSD represents pore size distribution.
To study the effect of pore size distribution on phase behavior using LB, a bundle of capillary tubes is used. The simuation setup is shown in Figure 6. Four capillary tubes are Energies 2021, 14, 1315 9 of 17 arranged in parallel with the side length of 4 nm, 6 nm, 8 nm, and 14 nm respectively. The length of tubes in x-direction is determined by the volume fraction occupied by each pore. Similar to the study in a single capillary tube, the solid boundary condition is applied in the y and z direction and the pressure boundary condition is applied in the x-direction. Then, the tubes are filled by increasing the gas density at the boundary. Since the tubes have different sizes and volumes, the volume of the bundle is the summation of the volume of each tube. We consider the average density of the bundle over the whole volume. With different sizes of tubes, condensation may happen first in small tubes. Due to the shared boundary, fluids in small tubes may exchange mass with fluids in other tubes. To study the effect of pore size distribution on phase behavior using LB, a bundle of capillary tubes is used. The simuation setup is shown in Figure 6. Four capillary tubes are arranged in parallel with the side length of 4 nm, 6 nm, 8 nm, and 14 nm respectively. The length of tubes in x-direction is determined by the volume fraction occupied by each pore. Similar to the study in a single capillary tube, the solid boundary condition is applied in the y and z direction and the pressure boundary condition is applied in the x-direction. Then, the tubes are filled by increasing the gas density at the boundary. Since the tubes have different sizes and volumes, the volume of the bundle is the summation of the volume of each tube. We consider the average density of the bundle over the whole volume. With different sizes of tubes, condensation may happen first in small tubes. Due to the shared boundary, fluids in small tubes may exchange mass with fluids in other tubes. Following the procedure in Section 3.1, the coexistence curve of methane in the bundle of capillary tubes can be established which is shown in Figure 7. Since it is not straightforward to directly apply pore size distribution into modified EOS or reservoir simulation, we choose to use one pore size to represent the effect of pore size distribution on confined phase behavior. As illustrated in Figure 7, the coexistence in the capillary bundle is close to the one in the 8 nm capillary tube. Thus, we consider 8 nm as the 'effective' pore size for this rock when calculating confined phase behavior.  Following the procedure in Section 3.1, the coexistence curve of methane in the bundle of capillary tubes can be established which is shown in Figure 7. Since it is not straightforward to directly apply pore size distribution into modified EOS or reservoir simulation, we choose to use one pore size to represent the effect of pore size distribution on confined phase behavior. As illustrated in Figure 7, the coexistence in the capillary bundle is close to the one in the 8 nm capillary tube. Thus, we consider 8 nm as the 'effective' pore size for this rock when calculating confined phase behavior. To study the effect of pore size distribution on phase behavior using LB, a bundle of capillary tubes is used. The simuation setup is shown in Figure 6. Four capillary tubes are arranged in parallel with the side length of 4 nm, 6 nm, 8 nm, and 14 nm respectively. The length of tubes in x-direction is determined by the volume fraction occupied by each pore. Similar to the study in a single capillary tube, the solid boundary condition is applied in the y and z direction and the pressure boundary condition is applied in the x-direction. Then, the tubes are filled by increasing the gas density at the boundary. Since the tubes have different sizes and volumes, the volume of the bundle is the summation of the volume of each tube. We consider the average density of the bundle over the whole volume. With different sizes of tubes, condensation may happen first in small tubes. Due to the shared boundary, fluids in small tubes may exchange mass with fluids in other tubes. Following the procedure in Section 3.1, the coexistence curve of methane in the bundle of capillary tubes can be established which is shown in Figure 7. Since it is not straightforward to directly apply pore size distribution into modified EOS or reservoir simulation, we choose to use one pore size to represent the effect of pore size distribution on confined phase behavior. As illustrated in Figure 7, the coexistence in the capillary bundle is close to the one in the 8 nm capillary tube. Thus, we consider 8 nm as the 'effective' pore size for this rock when calculating confined phase behavior.

Modified EOS from Pore-Scale Simulation
With shifted coexistence curve, a modified EOS can be built to account for the effects of confinement. Different forms of EOS have been proposed in the literature [35][36][37]. In this work, a modified EOS in Equation (23) is used to match with the simulation data [2].
The modified a and b account for the changed impulsive and repulsive forces in confined space. By carefully tuning the value of the above parameters, the coexistence curve from the modified EOS can match with simulation data. As shown in Figure 8, the modified EOS can capture the shift of the coexistence curve correctly. Compared to results from LB simulation, vapor densities from modified EOS have a good match with simulation. Liquid density is well predicted at high temperatures but overestimated at low temperatures. Note that the EOS used for matching the simulation data is not unique. Readers can use any form of EOS and a better match could be obtained.

Modified EOS from Pore-Scale Simulation
With shifted coexistence curve, a modified EOS can be built to account for the effects of confinement. Different forms of EOS have been proposed in the literature [35][36][37]. In this work, a modified EOS in Equation (23) is used to match with the simulation data [2].
The modified and account for the changed impulsive and repulsive forces in confined space. By carefully tuning the value of the above parameters, the coexistence curve from the modified EOS can match with simulation data. As shown in Figure 8, the modified EOS can capture the shift of the coexistence curve correctly. Compared to results from LB simulation, vapor densities from modified EOS have a good match with simulation. Liquid density is well predicted at high temperatures but overestimated at low temperatures. Note that the EOS used for matching the simulation data is not unique. Readers can use any form of EOS and a better match could be obtained. The modified EOS is then extended to model the confined phase behavior of shale oil. The fluid model used is shown in Tables 1 and 2. The P-T phase diagram of shale oil using the modified EOS is presented in Figure 9. The pore diameter used for thermodynamic modeling is 8 nm. In the literature, implementing capillary pressure and critical shift method are the two commonly used methods for confined phase behavior. To make a comparison, we also plot the phase diagrams using the two methods as well as that at the bulk condition in Figure 9. Compared to bulk, all three methods generate shifted phase diagrams, and suppressed bubble point pressures are observed. However, the method of implementing capillary pressure cannot capture the changed critical point. The modified EOS in this study and critical shift method using Equations (18) and (19) can both capture the shift of critical point. Qualitatively, the shift is more significant when using the modified EOS. In Figure 10, the saturation pressures in different sizes of pores are compared when T = 240 °F. It is seen that saturation pressures decrease linearly with the increasing inverse of pore size. Also, the suppression of saturation pressure is more significant in all sizes of pores when using the modified EOS. The modified EOS is then extended to model the confined phase behavior of shale oil. The fluid model used is shown in Tables 1 and 2. The P-T phase diagram of shale oil using the modified EOS is presented in Figure 9. The pore diameter used for thermodynamic modeling is 8 nm. In the literature, implementing capillary pressure and critical shift method are the two commonly used methods for confined phase behavior. To make a comparison, we also plot the phase diagrams using the two methods as well as that at the bulk condition in Figure 9. Compared to bulk, all three methods generate shifted phase diagrams, and suppressed bubble point pressures are observed. However, the method of implementing capillary pressure cannot capture the changed critical point. The modified EOS in this study and critical shift method using Equations (18) and (19) can both capture the shift of critical point. Qualitatively, the shift is more significant when using the modified EOS. In Figure 10, the saturation pressures in different sizes of pores are compared when T = 240 • F. It is seen that saturation pressures decrease linearly with the increasing inverse of pore size. Also, the suppression of saturation pressure is more significant in all sizes of pores when using the modified EOS.  Table 2. Binary interaction coefficient of studied fluid [34].   Table 2. Binary interaction coefficient of studied fluid [34].

Compositional Simulation of Shale Reservoir Considering Confined Phase Bahavior
As illustrated in Section 3.3, the shifted phase diagram and saturation pressure of shale oil indicate that the properties of oil and gas should also be altered in nanopores, which in turn will affect the predicted oil and gas production. To evaluate the production of shale oil reservoirs using the different EOS, the fractured shale reservoir shown in Figure 11 is studied. As shown, a horizontal well is drilled in the middle of the 2D reservoir (shown as a bold black line). The red lines represent 17 stages of hydraulic fractures that are connected with the well. The properties of shale matrix and hydraulic fractures are summarized in Table 3. This small reservoir should be regarded as the SRV (stimulated reservoir volume). Apply no-flow boundary to the outer bound of this reservoir and production is constrained by constant bottom-hole pressure. As illustrated in Figure 8, we use 8 nm as the effective pore size for the shale matrix. The bubble point pressure of this oil sample in 8 nm pore is 2630 psi, 2345 psi, and 2105 psi when using methods of capillary pressure, critical shift, and modified EOS respectively. At bulk condition, the bubble point pressure is 2768 psi. At initial reservoir pressure (3700 psi), the fluid shale matrix is thus pure oil. At a later stage of production, gas starts to appear and three-phase flow will be involved. The three-phase relative permeability in Figure 12 is applied for flow in the shale matrix. In the fractures, assume that the relative permeabilities are a linear function with saturations.

Compositional Simulation of Shale Reservoir Considering Confined Phase Bahavior
As illustrated in Section 3.3, the shifted phase diagram and saturation pressure of shale oil indicate that the properties of oil and gas should also be altered in nanopores, which in turn will affect the predicted oil and gas production. To evaluate the production of shale oil reservoirs using the different EOS, the fractured shale reservoir shown in Figure 11 is studied. As shown, a horizontal well is drilled in the middle of the 2D reservoir (shown as a bold black line). The red lines represent 17 stages of hydraulic fractures that are connected with the well. The properties of shale matrix and hydraulic fractures are summarized in Table 3. This small reservoir should be regarded as the SRV (stimulated reservoir volume). Apply no-flow boundary to the outer bound of this reservoir and production is constrained by constant bottom-hole pressure. As illustrated in Figure 8, we use 8 nm as the effective pore size for the shale matrix. The bubble point pressure of this oil sample in 8 nm pore is 2630 psi, 2345 psi, and 2105 psi when using methods of capillary pressure, critical shift, and modified EOS respectively. At bulk condition, the bubble point pressure is 2768 psi. At initial reservoir pressure (3700 psi), the fluid shale matrix is thus pure oil. At a later stage of production, gas starts to appear and three-phase flow will be involved. The three-phase relative permeability in Figure 12 is applied for flow in the shale matrix. In the fractures, assume that the relative permeabilities are a linear function with saturations.    (a) (b) Figure 12. Oil-water relative permeability (a) and oil-gas relative permeability (b) in shale matrix.
Simulated production from the REV using different confined models are presented in Figure 13 including reservoir pressure, cumulative oil production, and cumulative gas production. As shown in Figure 13a, compared to bulk condition, implementing confined phase behavior leads to a faster pressure decline. In terms of the rate of pressure decline, there is not much difference between the methods of capillary pressure, critical shift, and modified EOS used in this study. However, the oil and gas production is quite different when using different methods. As illustrated in Figure 13b, the cumulative oil production is higher if confined phase behavior is considered. The use of modified EOS in this study leads to the highest oil production. Oil production considering critical shift is second high followed by that using capillary pressure. The increased oil production rate is due to the suppressed bubble point in nanopores. It is apparent that more oil production can be expected if the suppression of bubble point is more significant. On the other hand, gas production is lower than that at bulk condition if only capillary pressure is considered. However, gas production is higher if the shift of critical points or modified EOS in this study is applied. Cumulative gas productions are close between the methods of critical shift and modified EOS. The comparison of cumulative gas productions using different method are presented in Figure 13c.
Since we only studied a simplified reservoir model in this work, a direct comparison with field data is not possible. The aim of this work is to provide an insight for simulation of shale reservoir in terms of altered phase behavior. As discussed above, the effects of different methods of confined phase behavior on oil/gas production from shale reservoirs cannot be ignored. The method of implementing capillary pressure fails to capture the altered critical points in nanopores. Considering capillary pressure alone is not sufficient to model confined phase behavior and both oil and gas production will be underestimated. Applying the critical shift method from Zarragoicoechea and Kuz [32] and modified EOS leads to higher oil and gas production. Besides, it is observed that the modified EOS has close performance with the critical shift method in predicting gas production but different in predicting oil production. The difference between these two methods requires more work in the future including comparison with experimental results or other methods. Simulated production from the REV using different confined models are presented in Figure 13 including reservoir pressure, cumulative oil production, and cumulative gas production. As shown in Figure 13a, compared to bulk condition, implementing confined phase behavior leads to a faster pressure decline. In terms of the rate of pressure decline, there is not much difference between the methods of capillary pressure, critical shift, and modified EOS used in this study. However, the oil and gas production is quite different when using different methods. As illustrated in Figure 13b, the cumulative oil production is higher if confined phase behavior is considered. The use of modified EOS in this study leads to the highest oil production. Oil production considering critical shift is second high followed by that using capillary pressure. The increased oil production rate is due to the suppressed bubble point in nanopores. It is apparent that more oil production can be expected if the suppression of bubble point is more significant. On the other hand, gas production is lower than that at bulk condition if only capillary pressure is considered. However, gas production is higher if the shift of critical points or modified EOS in this study is applied. Cumulative gas productions are close between the methods of critical shift and modified EOS. The comparison of cumulative gas productions using different method are presented in Figure 13c.
Since we only studied a simplified reservoir model in this work, a direct comparison with field data is not possible. The aim of this work is to provide an insight for simulation of shale reservoir in terms of altered phase behavior. As discussed above, the effects of different methods of confined phase behavior on oil/gas production from shale reservoirs cannot be ignored. The method of implementing capillary pressure fails to capture the altered critical points in nanopores. Considering capillary pressure alone is not sufficient to model confined phase behavior and both oil and gas production will be underestimated. Applying the critical shift method from Zarragoicoechea and Kuz [32] and modified EOS leads to higher oil and gas production. Besides, it is observed that the modified EOS has close performance with the critical shift method in predicting gas production but different in predicting oil production. The difference between these two methods requires more work in the future including comparison with experimental results or other methods.

Conclusions
In this work, we conducted a pore-scale simulation of confined phase behavior using the pseudo-potential LB method. The model was validated with Peng-Robinson equation of state at bulk condition. Vapor-liquid equilibrium in nanopores was simulated by modelling adsorption isotherms in capillary tubes. A sudden jump in adsorption isotherms indicated the occurrence of capillary condensation in nanopores. In addition, effect of pore size distribution on phase equilibrium is evaluated by using a bundle of capillary tubes with varied sizes. Size and volume fraction of each tube is determined by matching with experimental data of pore size. To include the deviated phase equilibrium into a reservoir-scale simulation, a modified EOS is built guided by results from LB simulation. The following conclusions can be made in this study: • By comparing constructed coexistence curves, we propose that an effective pore size can be used to represent a real rock sample with distributed nanopores.

•
Shifted critical properties and suppressed bubble points were observed when using the modified EOS. • Compared to methods using capillary pressure and critical shift, the phase diagram using modified EOS shrinks more significantly. Figure 13. Comparison of reservoir performance from compositional simulation using methods of capillary pressure, critical shift, and modified EOS; (a) reservoir pressure, (b) cumulative oil production, and (c) cumulative gas production.

Conclusions
In this work, we conducted a pore-scale simulation of confined phase behavior using the pseudo-potential LB method. The model was validated with Peng-Robinson equation of state at bulk condition. Vapor-liquid equilibrium in nanopores was simulated by modelling adsorption isotherms in capillary tubes. A sudden jump in adsorption isotherms indicated the occurrence of capillary condensation in nanopores. In addition, effect of pore size distribution on phase equilibrium is evaluated by using a bundle of capillary tubes with varied sizes. Size and volume fraction of each tube is determined by matching with experimental data of pore size. To include the deviated phase equilibrium into a reservoirscale simulation, a modified EOS is built guided by results from LB simulation. The following conclusions can be made in this study: • By comparing constructed coexistence curves, we propose that an effective pore size can be used to represent a real rock sample with distributed nanopores. • Shifted critical properties and suppressed bubble points were observed when using the modified EOS. • Compared to methods using capillary pressure and critical shift, the phase diagram using modified EOS shrinks more significantly.
• Compositional simulation indicated that both oil and gas production will increase if modified EOS is implemented. Considering capillary pressure or shifted critical point only will underestimate oil and gas production and ultimate recovery from shale reservoir.