A Theoretical Analysis of Pore Size Distribution Effects on Shale Apparent Permeability

Apparent permeability is an important input parameter in the simulation of shale gas production. Most apparent permeability models assume a single pore size. In this study, we develop a theoretical model for quantifying the effect of pore size distribution on shale apparent permeability. The model accounts for the nonuniform distribution of pore sizes, the rarefaction effect, and gas characteristics. The model is validated against available experimental data. Theoretical calculations show that the larger the pore radius, the larger the apparent permeability. Moreover, the apparent permeability increases with an increase in the width of pore size distribution, with this effect being much more pronounced at low pressure than at high pressure.


Introduction
The unconventional gas production boom in North America has attracted increasing interest internationally [1].Unconventional gas production has been the topic of much research, but its accurate prediction is still a challenge [2][3][4][5][6], especially for shale gas.The simulation of shale gas extraction requires information about the apparent permeability of shale, which deviates its intrinsic permeability due to the complex fluid dynamics in tight porous media.Gas transport in shale involves different transport mechanisms, including viscous flow, slip flow, Knudsen diffusion, and surface diffusion [7][8][9].A variety of studies have been done on studying these different mechanisms [10][11][12][13][14].However, most of studies consider shale media as a bundle of uniform capillaries with an effective pore radius  e .The determination of  e is not a trivial work.Some researchers employ the mercury injection method to determine  e .For conventional rocks, Winland [15] proposed using the pore size at 35% mercury saturation as  e .Pittman [16] extended Winland's work and pointed out that the pore size at 25% mercury saturation controls the permeability.For tight gas sands, Rezaee et al. [17] indicated that the pore size at 10% mercury saturation is a good permeability predictor.For shale media, very high pressure would be required for mercury to enter small pores around 3.6 nm [18].At such high pressures, the deformation of the rock and the possibility destruction of pore structure will affect the measurement [18].For determination of  e , a more sophisticated method is the Effective Medium Approximation (EMA).Ghanbarian and Javadpour [19] combined EMA with the Javadpour model [7] to estimate the apparent permeability of shale.Their model is very sensitive to  e that is determined by EMA.As seen from the above literature review, the determination of  e is still an open issue.
Shale media possesses a complex pore architecture [20].Thus, the uniform capillary model does not adequately account for the heterogeneity of shale.Civan [21] indicated that the apparent permeability of tight rock is directly related to pore size distribution (PSD).But limited work has studied the effect of PSD on gas transport in shale.Villazon et al. [22] employed a log-normal density function to characterize the PSD of shale and proposed an apparent permeability model.The model is developed based on the Hagen-Poiseuille equation and the general slip boundary condition [23].Thus, the model only considers viscous flow and slip flow.But, at low pore pressure, Knudsen diffusion dominates gas transport.The fractal capillary model has received much attention over several decades [24][25][26][27].Recently, some researchers used the fractal theory to describe the PSD of shale and developed a series of apparent permeability models [28][29][30].These models assume that the PSD of shale can be modeled as a power law distribution.Apart from a power law distribution, the PSD of shale can be described by a nonsymmetrical distribution, such as a gamma distribution [31][32][33][34].
Fewer studies have been carried out to theoretically investigate the effect of PSD on shale apparent permeability.In this article, an apparent permeability model that accounts for the nonuniform distribution of pore sizes is developed and validated with experimental results.The model can describe unique flow behaviors in shale and can be used to study the effect of PSD on shale apparent permeability.

Model Development
When the characteristic size of pore channels is comparable to or smaller than the gas mean free path, the molecule-wall collisions predominate over the intermolecular collisions, which is known as the rarefaction effect.At this condition, the Navier-Stokes equations break down [37].Thus, the applicability of the conventional Darcy law is invalid.Shale is referred to as extremely tight porous media with pore sizes in the nanoscale [38].Due to the rarefaction effect, shale gas transport is a complex process involving viscous flow, slip flow, and Knudsen diffusion [12,14].The Knudsen number quantifies the rarefaction of fluid.The Knudsen number for producing shale gas reservoirs is less than unity, which corresponds to the Darcy flow, slip flow, and early transitionflow regimes [39].In this work, we consider slip flow as a part of Knudsen diffusion [40].Moreover, we do not consider the surface diffusion in this work.For a single capillary, the molar flow rate for viscous flow  v is given by where  is the molar density, mol⋅m -3 ,  is the viscosity, Pa⋅s,  is the pore pressure, Pa,  is the ideal gas constant, J⋅mol -1 ⋅K -1 ,  is the temperature, K,  is the compressibility factor,  is the porosity,  is the tortuosity,  0 is the length of capillary, m, and  is the pore radius, m.The compressibility factor can be estimated using the Soave−Benedict−Webb−Rubin equation of state (SBWR-EOS) due to its accuracy at reservoir condition [41].Gas viscosity data are taken from the NIST database [42].
For Knudsen diffusion, the molar flow rate  K is [43] where  is gas molar weight, kg⋅mol -1 .
The molar flow rate  in a single capillary is [40,44] Since shale is assumed to be composed of parallel capillaries with different sizes, the total molar flow rate through unit area  T can be obtained using the following integral: where  t is the total number of capillaries per unit area, m -2 ,  max is the radius of the maximum capillary, m,  min is the radius of the minimum capillary, m, and () is the probability density function for capillaries (pores).In this work, we take  min as the gas molecular radius.The value of  max is taken as 200 nm [45].Moreover, following Moghaddam and Jamiolahmady [33], we assume that the PSD of shale can be modeled by a gamma distribution with a mean,  m = /, and variance  2 = / 2 , where  is the pore radius, m,  is the shape parameter,  is the rate parameter, m -1 ,  m is the mean pore radius, m, and Γ() is the gamma function.The parameters  and  can be any value greater than zero.Note that ( 8) is limited to a relatively narrow monomodal distribution of pore sizes.The pore size distributions described by the gamma density function are shown in Figure 1 for several values of  and .
As shown in Figure 1(a), the most likely pore radius (the pore radius corresponding to the maximum value of the gamma density function) is smaller than  m .Moreover, the difference between  m and most likely pore radius increases with an increase in the width of the distribution (see Figure 1(a), solid lines).This behavior is due to the asymmetry in the gamma distribution, which is characterized by a relatively long tail at large .For most porous materials, the PSD is not symmetric, but skewed [46].Based on the Darcy equation, the total molar flow rate through unit area is Comparing ( 6) and ( 9), the apparent permeability is obtained as

Model Validation with Experimental Data.
To eliminate or minimize the effect of gas adsorption/desorption, shale permeability tests with helium are used to validate our model (see (10)).For fitting the model to the experimental data, the , , and  are treated as fitting parameters. is also treated as a fitting parameter, unless its value is known.The database used in this work consists of experiments from Aljamaan et al. [35] and Mathur et al. [36].Aljamaan et al. [35] carried out permeability measurements on two shale samples from Barnett and Eagle Ford.They conducted the tests with helium, methane, nitrogen, and carbon dioxide at temperature 296.95 K.Note that the permeability measurements with sorbing gases (methane, nitrogen, and carbon dioxide) are not used to validate our model.Figure 2 shows the comparison between the model results and the experimental data for the Eagle Ford shale sample.Table 1 lists the fitted parameters. min for helium is 0.26 nm.As can be seen from Figure 2, our model captures the trends of the experimental data very well.Figure 2 also shows that the apparent permeability is high initially, drops steeply with increasing pore pressure, and approaches a constant value at high pore pressure (>5 MPa).This is because, at low pore pressure, the dominant transport mechanism is Knudsen diffusion, which enhances gas transport in shale [40].With increasing pore pressure, the effect of Knudsen diffusion vanishes, which results in the decrease of the apparent permeability.Furthermore, according to (10), at high pore pressure, the apparent permeability of shale converges to its intrinsic permeability and is not sensitive to pore pressure.
The present model is also compared with the experimental data from Mathur et al. [36].This dataset includes two shale samples (XX86 and XX88) from the Wolfcamp  [35].The solid line is calculated from (10).The parameter set is presented in Table 1.
formation.These measurements are performed under steady state at temperature 298.15 K. Helium and nitrogen are used as permeating fluids.Once again, we only use the helium measurements to validate the present model.Figure 3 shows that the present model matches the experimental data well.Error bars for the experimental data of Aljamaan et al. [35] and Mathur et al. [36] are omitted because experimental uncertainties are not reported in those studies.The fitted parameters are listed in Table 1.The  values fall within the range of 2.3−11.9 reported by Katsube et al. [47].Moreover, the  m values range from 2.27 nm to 41 nm, which is consistent with the findings of Javadpour et al. [38].They found that pore radiuses in the range of 2 nm to 100 nm are the dominant flow channels in shale.Figures 3 and 2 show that the permeability values of the Eagle Ford shale sample are significantly higher than those of the Wolfcamp shale samples.Intuitively, one should expect permeability to vary from one sample to another.Moreover, we note that, in Mathur et al. 's [36] experiment, the measurements are conducted on "as-received" samples.But, in Aljamaan et al. 's [35] experiment, the measurements are conducted on dried samples.Ghanizadeh et al. [48] indicated that permeability values measured on dried samples are higher than those measured on "as-received" samples, which is due to the presence of pore fluid (bound/adsorbed water and hydrocarbon) in "as-received" samples [48].The bound/adsorbed water and hydrocarbon partially block pore channels, which results in the decrease of permeability [49].

Gas Transport Dynamics in Shale.
The flow regimes for producing shale gas reservoirs include continuum flow, slip flow, and transition flow [40].The present model incorporates viscous flow and Knudsen diffusion to cover these flow regimes.To enunciate the gas transport dynamics in shale, we define the permeability ratios for viscous flow  rv and Knudsen diffusion  rk .
Figure 4 shows the permeability ratios  rv and  rk as a function of pore pressure.The results shown in Figure 4 are obtained with  = 5.33 and  = 0.13.The other parameters are listed in Table 2.As can be seen from Figure 4, the permeability ratio  rv increases with increasing pore pressure, and the permeability ratio  rk decreases with increasing pore pressure, which is consistent with the molecular simulation results of Hari et al. [50].This behavior is due to the fact that when pore pressure is increased, the gas molecules collide with each other more frequently than with the pore wall, and thus viscous flow contributes more to the total flow.Figure 4 also shows that when 0 MPa <  < 6 MPa,  rv increases steeply, while for  > 6 MPa  rv increases slowly.Specifically,  rv increases sharply from 0.327 to 0.709 when the pore pressure increases from 1 to 5 MPa.Meanwhile,  rk decreases sharply from 0.673 to 0.291.When the pore pressure increases from 10 to 50 MPa,  rv shows a minor increase from 0.832 to 0.967, and  rk decreases smoothly from 0.168 to 0.033.This is because that, according to (11), with increasing pore pressure  rv gradually converges to 1. Thus, at high pore pressure (>30 MPa),  rv is not sensitive to pore pressure.Moreover, at high pore pressure (>30 MPa), Knudsen diffusion is negligible.Thus, at this condition,  rk is also not sensitive to pore pressure.2. As shown in Figure 5(a), the apparent permeability increases with an increase in  m .Specifically, at the pore pressure of 1.2 MPa, the apparent permeability increases from 0.771 to 1.152 uD when  m increases from 2 to 4 nm.This behavior is due to the fact that the apparent  [36].The solid and dashed lines are calculated from (10).The parameter set is presented in Table 1.

Effect of Pore Size Distribution on Shale Apparent Permeability.
permeability is positively associated with pore radius.The larger the pore radius, the larger the apparent permeability.Figure 5(a) also shows that the effect of  m on the apparent permeability is much larger at low pore pressure than at high pore pressure.Specifically, the difference between  a3 and  a1 decreases from 0.742 to 0.092 uD when the pore pressure increases from 0.6 to 6 MPa.At low pore pressure, the difference between  a3 and  a1 is contributed by both viscous flow and Knudsen diffusion.But, at high pore pressure, Knudsen diffusion could be ignored.Thus, the increase in pore pressure decreases the effect of   on the apparent permeability.2. For the three considered distributions, the apparent permeability increases with an increase in  2 .At the pore

Viscous flow
Figure 4: Permeability ratio as a function of pore pressure.The solid line is calculated from (11), and the dash dot line is calculated from (12).The parameter set is presented in Table 2.
pressure of 1.2 MPa, the apparent permeability increases from 0.863 to 1.315 uD when  2 increases from 0.75 to 3 nm 2 .This result is attributed to the width of the PSD being represented by  2 .The gamma distribution is a positive distribution.Thus, with increasing  2 , the PSD presents a lower pick at small pore sizes but a longer tail at large pore sizes, as shown in Figure 1(a).Moreover, the molar flow rates for viscous flow  v and for Knudsen diffusion  K have a strong dependence on pore radius. v is proportional to  4 (see (1)), and  K is proportional to  3 (see (3)), respectively.Thus, in the case of the same  m , increasing  2 results in the increase of the apparent permeability.Figure 5(b) also shows the effect of  2 on the apparent permeability is more pronounced under low pore pressure than under high pore pressure.This is because at high pore pressure Knudsen diffusion is negligible.The effect of PSD on gas transport in shale can also be seen from the behavior of the fractional flow rate  q , which is defined as In (13), the denominator characterizes the total flow rate, and the numerator characterizes the flow rate contributed by pores with radius  ≤   .Figure 6(a) shows the effect of  m on  q .The results shown in Figure 6(a) are obtained with  1 = 4,  1 = 2,  2 = 8,  2 = 8 1/2 ,  3 = 16, and  3 = 4.The corresponding values of mean and variance are given in the legend of Figure 6(a).The other parameters are listed in Table 2.As can be seen from Figure 6(a), in the case of the same  2 , increasing  m shifts the fractional flow curve toward the right, which means that the majority of the flow occurs through large pores.Specifically, at the pore pressure of 10 MPa,  q1 for the pore radius of 3 nm is 36.4%.At the same pore pressure,  q3 for the pore radius of 3 nm is only 3.3%, which means that more than 90% of the total flow was contributed by pores with radius  ≥ 3 nm.Moreover, at the conditions of pore pressure and pore radius of 10 MPa and Figure 5: Apparent permeability as a function of pore pressure.The lines are calculated from (10).The parameter set is presented in Table 2. (a) The effect of  m on the apparent permeability.(b) The effect of  2 on the apparent permeability.5 nm,  q1 is equal to 85.1%,  q2 is equal to 81.5%, and  q3 is equal to 59.5%.This is because, in the case of the same  2 , increasing  m means that more large pores are present.Thus, the main contribution to the total flow rate is from these large pores.Figure 6(a) also shows the effect of pore pressure on  q .As can be seen from Figure 6(a), decreasing the pore pressure shifts the fractional flow curve toward the left.When the pore pressure is 10 MPa and the pore radius is 4 nm,  q1 = 65.7%, q2 = 54.9%, and  q3 = 23.9%.When the pore pressure is 0.01 MPa and the pore radius is 4 nm,  q1 = 68.7%, q2 = 57.7%, and  q3 = 25.8%.This is because Knudsen diffusion dominates gas transport at low pore pressure, which enhances gas flow.Thus, at the same pore radius,  q at low pore pressure is larger than that at high pore pressure.
Figure 6(b) shows the sensitivity of  q to the variance of PSD  2 .The results shown in Figure 6(b) were obtained with  1 = 3,  1 = 1,  2 = 6,  2 = 2,  3 = 12, and  3 = 4.The corresponding values of  m and  2 are given in the legend of Figure 6(b).As shown in Figure 6(b), in the case of the same  m , increasing  2 shifts the fractional flow curve toward the right.Specifically, at the pore pressure of 10 MPa,  q1 for the pore radius of 4 nm is 18.2%.At the same pore pressure,  q3 for the pore radius of 4 nm is 61%.This behavior is explained by the fact that decreasing  2 improves the uniformity of the PSD; namely, the width of the PSD decreases with decreasing  2 ; see Figure 1(a) (solid lines).Figure 1(a) (solid lines) also shows that, with decreasing  2 , the most likely pore radius approaches the mean pore radius.When  2 equals 0, the pore size is uniform and equal to  m .At this condition,  q0 shows a step increase from 0 to 1 at  m , as shown in Figure 6(b).Thus, in the case of the same  m , with decreasing  2 ,  q increases more steeply at  m .Figure 6(b) also shows that  q at low pore pressure is larger than that at high pore pressure.It is worth noting that the difference between  q at low pore pressure and that at high pore pressure increases with increasing  2 .This is explained by the fact that a larger  2 corresponds to a higher quantity of extremely small pores ( ≤ 1 nm), as shown in Figure 1(b) (solid lines).At low pore pressure, the smaller the pore radius, the stronger the molecule−wall collisions and the stronger the rarefaction effect, which enhances gas transport.

Conclusion
In this paper, we have studied the effect of PSD on shale apparent permeability.We have derived an apparent permeability model with the assumption that shale media consists of a bundle of tortuous capillaries with a gamma distribution.The apparent permeability model matches experimental data well.Although the model is not completely realistic, it allows us to understand the effect of PSD on shale apparent permeability.Moreover, the proposed model possesses an analytical form and can be easily incorporated into existing reservoir simulators without large code change.The major findings are as follows.In the case of the same variance of PSD, increasing  m increases the apparent permeability.In the case of the same  m , the apparent permeability increases with an increase in the width of PSD.Moreover, the apparent permeability is more sensitive to PSD at low pore pressure than at high pore pressure, because Knudsen diffusion enhances gas flow at low pore pressure and is negligible at high pore pressure.It would be interesting to generalize the present work by considering the effect of surface diffusion and stress sensitivity.This will be done in future research.Figure 6: Fractional flow rate as a function of pore radius.The lines are calculated from (13).The parameter set is presented in Table 2.

(a)
The effect of  m on the fractional flow rate.(b) The effect of  2 on the fractional flow rate.

Figure 1 :
Figure 1: Representative plots of the gamma distribution: (a) probability density function, (b) cumulative distribution function.

Figure 2 :
Figure 2: Comparison of the model results with the experimental data of Aljamaan et al.[35].The solid line is calculated from(10).The parameter set is presented in Table1.

Figure 5 (
a) shows how the mean value of distribution affects shale apparent permeability.The results shown in Figure 5(a) were obtained with  1 = 4,  1 = 2,  2 = 8,  2 = 8 1/2 ,  3 = 16, and  3 = 4.The corresponding values of  m and  2 are given in the legend of Figure 5(a).The other parameters are listed in Table

Figure 3 :
Figure 3: Comparison of the model results with the experimental data of Mathur et al.[36].The solid and dashed lines are calculated from(10).The parameter set is presented in Table1.

Figure 5 (
b) shows how the variance of distribution affects the apparent permeability.The results shown in Figure 5(b) were obtained with  1 = 3,  1 = 1,  2 = 6,  2 = 2,  3 = 12, and  3 = 4.The corresponding values of  m and  2 are given in the legend of Figure 5(b).The other parameters are presented in Table

Table 1 :
Fitted parameters of the apparent permeability model.

Table 2 :
Parameters for the sensitivity analysis.