The Brooks and Corey Capillary Pressure Model Revisited from Pore Network Simulations of Capillarity-Controlled Invasion Percolation Process

Relating the macroscopic properties of porous media such as capillary pressure with saturation is an on-going problem in many fields, but examining their correlations with microstructural traits of the porous medium is a challenging task due to the heterogeneity of the solid matrix and the limitations of laboratory instruments. Considering a capillarity-controlled invasion percolation process, we examined the macroscopic properties as functions of matrix saturation and pore structure by applying the throat and pore network model. We obtained a relationship of the capillary pressure with the effective saturation from systematic pore network simulations. Then, we revisited and identified the microstructure parameters in the Brooks and Corey capillary pressure model. The wetting phase residual saturation is related to the ratio of standard deviation to the mean radius, the ratio of pore radius to the throat length, and pore connectivity. The size distribution index in the Brooks and Corey capillary pressure model should be more reasonably considered as a meniscus size distribution index rather than a pore size distribution index, relating this parameter with the invasion process and the structural properties. The size distribution index is associated with pore connectivity and the ratio of standard deviation to mean radius (σ0/r¯), increasing with the decline of σ0/r¯ but the same for networks with same σ0/r¯. The identified parameters of the Brooks and Corey model might be further utilized for correlations with other transport properties such as permeability.


Introduction
Invasion percolation (IP), introduced by Wilkinson and Willemsen [1], provides a realistic description of the slow biphasic fluid displacement processes in porous materials. IP is relevant to various applications (such as secondary oil recovery) and was generally used to model a drainage process during which a wetting fluid is displaced by a non-wetting fluid [2]. IP was also used to model imbibition, where a wetting fluid invades a porous medium originally saturated by a non-wetting fluid [3]. Moreover, the IP algorithm was adapted to describe immiscible displacement in drying porous media [4].
Invasion percolation is driven by the capillary pressure P c , which can be modeled with the normalized or effective wetting phase saturation as proposed by Brooks and Corey [5]: Processes 2020, 8 S e = P e P c λ (P c ≥ P e ) (2) where entry pressure P e and S e denote the entry pressure and effective saturation, respectively. The scaling parameter λ denotes the pore size distribution index, which must be determined experimentally for each porous medium. The entry pressure (P e ) describes the minimum air pressure that needs to act on the liquid phase prevailing in the medium to allow the air to invade the porous medium. S e is a linear function of saturation S and is related to the residual saturation S r (or, in the context of drying porous media, irreducible saturation), defined as: The irreducible saturation S r defines the amount of liquid that the medium with a given pore structure can retain. This liquid is held in the medium by capillary forces and cannot be removed by the capillary pressure.
To investigate whether the irreducible saturation is correlated with the microstructure of porous media, it may be worth making an analogy between the irreducible saturation in the capillarity-controlled invasion process and the mercury entrapment in porosimetry as these processes share common principles. In porosimetry, mercury cannot be retracted from a large pore unless a meniscus appears at one end of that pore so that mercury can penetrate to the surface. Thus, if a large pore is embedded in a network of fine pores, it remains full unless some fine pores evacuate at some lower extrusion pressure-this effect is referred to as shielding effect [6,7]. Portsmouth and Gladden [8] remarked that there is a characteristic entrapment of mercury for a network comprising Gaussian pore size distribution, which depends on the ratio of the standard deviation of the distribution to the mean pore radius-this ratio defines the shape of pore size distribution (PSD). Furthermore, Conner and Lane [9] performed analyses and have concluded that higher connectivity can yield lower entrapment of mercury. Other researchers have also disclosed the fact that topology of the medium and spatial correlations between pores and throats can influence porosimetry results [8,10,11]).
Nevertheless, the irreducible saturation does not have the same correlations with microstructural features as the mercury entrapment. Isolated mercury from the bulk of the network is forced to move if the Laplace equation is fulfilled along with the continuous decreasing of the non-wetting phase pressure. That is to say, the mercury in sites surrounded with narrow void spaces is trapped in the medium when the Laplace equation is not satisfied. However, this situation is different in capillarity-controlled invasion processes. For instance, in slow isothermal drying an isolated liquid cluster is considered as stationary when the liquid in this cluster cannot move into neighbor clusters by improving the capillary pressure that acts on the network, but the liquid can continuously be reduced by evaporation so that the moisture is transported through the vapor phase. Neglecting evaporation at the isolated clusters would mean that the liquid in the isolated clusters is not reducible because these clusters are surrounded by the gas phase saturated with vapor. In such situation, the residual saturation in the network represents the saturation for which liquid can be transferred by the capillary pressure out of the porous medium. Thus, mercury porosimetry can provide information on the geometry and topology within the porous medium, but the irreducible saturation indicates a transport property.
In addition to the irreducible saturation, the Brooks and Corey capillary pressure model also requires the entry pressure P e (or bubbling pressure). The entry pressure, as defined by White et al. [12] is the capillary pressure at the inflection point on the capillary pressure-saturation (P c -S) curve, where the value of dP c /dS is a minimum at high saturation. Early numerical studies conducted at the pore scale show that the entry pressure depends on the interplay of the pore geometry and wettability of the porous medium (see, e.g., [13] and references therein). Furthermore, on the basis of pore-scale simulations, Raeini et al. [14] suggested that the entry pressure increases as the pore body to pore throat contraction ratio increases. Ferrari and Lunati [15] assumed perfect water-wet conditions, similar to Raeini et al. [14] but contrary to Rabbani et al. [13], and conducted two-phase flow numerical It is assumed that all throats and pores are initially filled with free capillary water. A gas reservoir is located at the top of the network (inlet). At the bottom, the network is connected to a water reservoir (outlet) (Figure 2). Air invades from the top of the network by IP algorithm. Periodic boundary conditions are applied to the lateral faces of the network. The effects of gravity and liquid viscosity are neglected in this study. Thus, capillary forces completely control the fluid configurations. Network representation by the throat-pore model means that the invasion process happens at menisci throats and pores. Only the invasion process at throats is however considered in the throat- It is assumed that all throats and pores are initially filled with free capillary water. A gas reservoir is located at the top of the network (inlet). At the bottom, the network is connected to a water reservoir (outlet) (Figure 2). Air invades from the top of the network by IP algorithm. Periodic boundary conditions are applied to the lateral faces of the network. The effects of gravity and liquid viscosity are neglected in this study. Thus, capillary forces completely control the fluid configurations. It is assumed that all throats and pores are initially filled with free capillary water. A gas reservoir is located at the top of the network (inlet). At the bottom, the network is connected to a water reservoir (outlet) (Figure 2). Air invades from the top of the network by IP algorithm. Periodic boundary conditions are applied to the lateral faces of the network. The effects of gravity and liquid viscosity are neglected in this study. Thus, capillary forces completely control the fluid configurations. Network representation by the throat-pore model means that the invasion process happens at menisci throats and pores. Only the invasion process at throats is however considered in the throat-  Network representation by the throat-pore model means that the invasion process happens at menisci throats and pores. Only the invasion process at throats is however considered in the throat-node model, assuming perfect water-wet conditions (contact angle is 0 • ) [4,24,29,30]. The liquid transition between pore and throat is neglected in the current work. Once a throat meniscus reaches the entrance of a pore, perfect piston-like meniscus movement is also assumed in the pore invasion step, whereby thermodynamic equilibrium is considered. This is because Tsakiroglou and Payatakes [10] illustrated that mercury menisci could reach equilibrium not only at the entrance to the throat but also at the entrance to the pore. Here, we assumed that the pore radius is larger (or equal) than the radii of the neighbor throats, leading to smaller capillary pressure required for the air to invade into the pore. The Young-Laplace equation is traditionally used to estimate capillary pressure in menisci throats or pores: where θ is the equilibrium contact angle, σ (N/m) and r (m) stand for the liquid surface tension, the radius of the meniscus (where the throat radius is used for the meniscus throat, and the pore radius is assumed for the meniscus pore), respectively. The capillarity-controlled invasion percolation algorithm implemented here follows the original invasion percolation procedure proposed by Wilkinson and Willemsen [1]: (1) Initially, all the pores and throats are completely filled with water. Neglecting pressure difference due to liquid head, the liquid pressure in the network and the liquid reservoir are set to zero for simplicity, without affecting the final capillary-pressure relationship. The air pressure necessary for any of the interior throats to be displaced with air is computed. We use Equation (5) to find the pressure difference between air and water necessary for a meniscus to penetrate the throat: where P c (Pa) is the capillary pressure, P g (Pa) is the gas pressure and P l (Pa) denotes the liquid pressure. Thus, the capillary pressure is equal to the gas pressure.
The pores at the inlet reservoir do not contain liquid ( Figure 2). Thus, the network initially has a meniscus in a throat with radius equal to the radius of that throat. Once a single throat gets unfilled of liquid, it forms a meniscus pore. The radius of the meniscus in the pore is considered as the radius of this pore due to the assumed pore body shape. Therefore, if a throat is unfilled of water, the adjacent downstream pore is also unfilled completely at the same time, creating new menisci at the other throat entrances (filled with water) attached to that pore body. Accordingly, the water phase in the attached throats becomes disintegrated at that pore.
(2) For the menisci that have been formed in the previous step, the capillary threshold is calculated by Equation (4). The meniscus throats with the smallest capillary thresholds are detected. (3) The air pressure is increased to the smallest capillary threshold calculated in step 2 (called entry pressure P e ) to invade accessible throats with that threshold and create meniscus throats. (4) Once a meniscus pore is invaded, it forms the active capillary thresholds at neighbor throats, which still contain liquid. When a meniscus throat is invaded, it forms a new active meniscus pore. Hence, the newly formed menisci (meniscus throats/pores) are identified, along with calculating their capillary threshold pressures. (5) The network saturation is calculated by integrating all residual liquid volume in the network and dividing it by the total void space volume. The smallest capillary threshold pressure is then assigned to the gas pressure. This is one data point in the capillary pressure-saturation curve. (6) Due to the random radius distribution of throats/pores, the initially backbone liquid cluster may be divided into several liquid clusters (may be isolated from the bottom reservoir). We account for the number of liquid clusters by a variant of the Hoshen-Kopelman algorithm [31]. The invasion percolation happens simultaneously in each cluster, which is still connected to the liquid reservoir according to step 3. If the invading pressure is greater than the capillary threshold of any of them, they are also invaded. (7) This algorithm is repeated from steps 2 to 6 until the clusters get isolated from the liquid reservoir, which is connected to the network bottom. This means all the pores at the bottom layer are unfilled of water. The remaining liquid in the network is considered as the irreducible saturation.
The Monte Carlo method with 40 realizations is used to reduce the uncertainty of the reported parameters.

Wetting Phase Residual Saturation (S r )
The wetting phase residual saturation shares the physical characteristics of mercury entrapment. In this regard, studies are reported in Section 1, which indicates that the mercury entrapment varies with the width of the pore size distribution as well as with the connectivity of the network. With this background, we simulated networks with different ranges of the mean radius, various widths of the pore size distribution as well as different averaged connectivity of pores. Figure 3 shows the wetting phase residual saturation as function of structural characteristics of the pore network structure exposed to the capillarity-controlled invasion percolation (IP) process.
Processes 2020, 8, x FOR PEER REVIEW 6 of 17 (7) This algorithm is repeated from steps 2 to 6 until the clusters get isolated from the liquid reservoir, which is connected to the network bottom. This means all the pores at the bottom layer are unfilled of water. The remaining liquid in the network is considered as the irreducible saturation.
The Monte Carlo method with 40 realizations is used to reduce the uncertainty of the reported parameters.

Wetting Phase Residual Saturation ( )
The wetting phase residual saturation shares the physical characteristics of mercury entrapment. In this regard, studies are reported in Section 1, which indicates that the mercury entrapment varies with the width of the pore size distribution as well as with the connectivity of the network. With this background, we simulated networks with different ranges of the mean radius, various widths of the pore size distribution as well as different averaged connectivity of pores. Figure 3 shows the wetting phase residual saturation as function of structural characteristics of the pore network structure exposed to the capillarity-controlled invasion percolation (IP) process.  A similar level of is observed in Figure 3a for the networks with an identical value of 0 / ̅ = 0.1 and two different values of ̅ = 10 µ m and 100 µ m. This may infer that is independent of the pore size. By contrast, different morphologies at 0 / ̅ = 0.1 show that decreases with increase in the ratio of ̅ ̅ ⁄ (Figure 3a), the ratio ̅ ̅ ⁄ indicates the role of spherical pore volume in the capillary invasion process, corresponding to the ratio of the spherical pore volume over the volume of the cylindrical throat.
To study whether irreducible saturation has characteristics that are dependent on the shape of PSD, the value of the ratio of standard deviation over the mean radius ( 0 / ̅ ) is also systematically varied. Figure 3a shows that decreases with increasing 0 / ̅ . This result is contradictory with the A similar level of S r is observed in Figure 3a for the networks with an identical value of σ 0 /r = 0.1 and two different values of r = 10 µm and 100 µm. This may infer that S r is independent of the pore size. By contrast, different morphologies at σ 0 /r = 0.1 show that S r decreases with increase in the ratio of r p /L t (Figure 3a), the ratio r p /L t indicates the role of spherical pore volume in the capillary invasion process, corresponding to the ratio of the spherical pore volume over the volume of the cylindrical throat. To study whether irreducible saturation has characteristics that are dependent on the shape of PSD, the value of the ratio of standard deviation over the mean radius (σ 0 /r) is also systematically varied. Figure 3a shows that S r decreases with increasing σ 0 /r. This result is contradictory with the results obtained by Portsmouth and Gladden [8], who reported the mercury entrapment increases when expanding the TSD. One should keep in mind that S r in the current work corresponds to that of a drying process that takes place under capillarity-controlled IP. When gas invades the network, the liquid phase splits up into several regions with different sizes, the so-called liquid clusters that consist of at least one nonempty throat connected to one nonempty pore, or even forms single liquid throats/pores surrounded by air, namely single menisci. The isolated clusters and single menisci, which are separated from the bottom reservoir, cannot be drained to the bottom reservoir by invasion percolation. The remaining liquid in these isolated clusters and single menisci could be removed by evaporation when a vapor pressure gradient is present in the porous medium. However, in the capillarity-controlled period of drying, the porous medium is saturated with vapor, which leads to the liquid in isolated clusters and single menisci being trapped in the porous medium. Thus, the isolated menisci clusters and single menisci are not of interest in the current work. Comparing the liquid before and after invasion in Figure 4, it can be seen that the majority of remaining liquid is trapped in throats with radii smaller than the mean radius. The mercury extrusion algorithm is different from the IP algorithm, which can partially describe the drying process since the trapped mercury could be reshaped by refilling the neighboring throats/pores, meeting the demands of the Laplace/Washburn equation [8]. Trapped mercury in the porous medium could express the shielding effect, being surrounded by the smallest throats in the network. The pore radius is assumed larger than the radius of neighbor throats, so that the pore could be trapped by the smaller throats. Therefore, the probability of mercury entrapment is enhanced by expanding the TSD. This might explain why values of S r simulated for the capillarity-controlled invasion process are contradictory with observations on the mercury entrapment from porosimetry. Regarding drying, Lai et al. [32] tested the irreducible saturation for six rock samples, showing that it decreases with increasing spread of TSD/PSD. Their observation indirectly validates our simulations.
The wetting phase residual saturation shows a positive linear correlation with the averaged pore connectivity in the network (Figure 3b). In each network, throats are randomly plugged, and the averaged connectivity is calculated by averaging the connectivities of all pores in the network. The pore connectivity affects the porosity of the network, as the porosity (ε) increases with plugging more throats onto a pore, which means that the void space increases during this procedure. Here, porosity is calculated as ε = V void V total , where V void is the void space of the porous medium and V total is the total volume of the porous medium. Figure 4 depicts the volume distributions of liquid-filled throats and pores in a network at the initial stage (S = 1) and at the final stage (S = S r ) with different r p /L t . At r p /L t = 0.01, the pores contain a very small amount of liquid that can be neglected when compared with liquid stored in the throats (Figure 4a). The spherical pore volume contribution increases with increasing pore radius (from Figure 4a-d). By comparing the liquid-filled volumes before and after the invasion, the fraction of invaded pores is seen to be larger than that of invaded throats. This could demonstrate that the pores play a significant role when they store a large fraction of liquid volume in the porous medium. The effect of throat length has been reported by Sok et al. [33] based on stochastically generated networks with coordination number distribution (Z = 4) specified from Fontainebleau sandstone sample, showing that increasing the throat length leads to a contrary effect on the residual saturation.
Moreover, the porosity of network is calculated by dividing all void space volume by the total space volume that is assumed as constant, leading to an increase along with the pore connectivity. Such additional volume might overcome the effect of pore connectivity. Thus, normalized wetting phase residual saturation is proposed as the ratio of the wetting phase residual saturation to the unity of the porous porosity. Figure 5 exhibits the relation of wetting phase residual saturation within networks containing the same void space but possessing various pore connectivity. This is inversely proportional to the pore connectivity, which agrees with experimental observations by Ojha et al. [34]. Increasing pore connectivity provides more pathways for the liquid to expel the network. With increased pore connectivity, more throat menisci are formed in the network once one liquid pore has completely been invaded. As a result, it enhances the IP as more throat menisci provide more options for air invasion. Such observations can be supported by the percolation theory (see, e.g., [35,36]).   The residual saturation reported in the current work might be questioned for higher values than those reported by Sok et al. [33] and Ferreira et al. [37]. This discrepancy can partly be due to the fact that the present work does not account for the film effect. In fact, liquid films play a crucial role in the wetting phase residual saturation [37], because they reduce the resistance for liquid transport and thereby enhance the liquid hydraulic connectivity [38]. Another reason might be that the spread of size distributions in the present network simulations is smaller than in real porous media. Sok et al. [33] reported a pore volume range from 1 × 10 −5 mm 3 to 1 × 10 −2 mm 3 based on the structure of Fontainebleau sandstone. The enlarged TSD/PSD could decrease the irreducible saturation (see Figure 3).

The Entry Pressure (P e )
In this section, we show that TSD and PSD play critical roles in entry pressure. Figure 6 shows the relation of entry pressure P e and mean capillary pressure (P c,mean) with the porous medium properties. The entry pressure is computed using the algorithm described in Section 2. The mean capillary pressure is calculated by the Young-Laplace equation (Equation (4)), averaging the throat and pore radii in the network with perfect wetting assumption. Here, the mean capillary pressure P c,mean is a reference capillary pressure.
The networks in Figure 6a are generated for a constant standard deviation of the TSD/PSD and six different mean values of the TSD/PSD. As can be seen, P e is close to the P c,mean when the ratio 1/r is small. A small value of the ratio standard deviation (i.e., 5 µm) to the mean radius means that the standard deviation has a weak impact on the radius distribution so that the network can be viewed as morphologically uniform. Figure 6b clearly indicates that the entry pressure is inversely proportional to the standard deviation. Networks with a wider distribution are more easily invaded by the air since the probability of having a larger meniscus at the interface is higher. Figure 6c shows that the entry pressure is not related to the connectivity for network coordination numbers smaller than 6. The entry pressure increases with connectivity slightly up to Z = 6, but decreases drastically at Z = 8. The reason might be that network structure is completely altered from cubic (Z = 3 to 6) to the octahedral (Z = 8). For the cubic based networks (Z = 3 to 6), the vertical surface throats are not altered and surface pores are connected with only one throat. However, three throats are connected with each surface pore for the octahedron network. As a result, the air more easily invades the octahedrons network, leading to a lower value of entry pressure. As a final remark in this section, one can note that the mean capillary pressure varies linearly with the ratio 1/r but it remains uniform for various values of σ 0 and Z. Figures 7 and 8 show the variations of capillary pressure-effective saturation at various morphologies obtained from the pore network modeling (Section 2). Three microstructural parameters of the network are discussed: (1) mean throat/pore radius (r) (Figure 7), (2) standard deviation in PSD and TSD (σ 0 ) (Figure 8), (3) pore connectivity (Z) (Figure 8). They all affect the capillary-saturation profile of a porous medium in the capillarity-controlled IP.

The Entry Pressure (Pe)
In this section, we show that TSD and PSD play critical roles in entry pressure. Figure 6 shows the relation of entry pressure Pe and mean capillary pressure (Pc,mean) with the porous medium properties. The entry pressure is computed using the algorithm described in Section 2. The mean capillary pressure is calculated by the Young-Laplace equation (Equation (4)), averaging the throat and pore radii in the network with perfect wetting assumption. Here, the mean capillary pressure Pc,mean is a reference capillary pressure.  As can be seen in Figure 7, the capillary pressure strongly depends on the mean radius. With increasing mean radius, the value of capillary pressure decreases. This behavior is expected as the capillary pressure is inversely proportional to the mean radius according to the Young-Laplace equation (Equation (4)). Since the capillary pressure in the Brooks and Corey model is presented as a function of effective saturation, the saturation is converted into effective saturation (S e ) using Equation (3). One can see that the overall behavior of the capillary pressure-effective saturation profiles obtained in this work agrees with the experimental measurements by Hilpert and Miller [39] and simulation observations by Sweijen et al. [40].
The role of two other structural parameters in the behavior of capillary pressure-saturation is presented in Figure 8. As shown, the capillary pressure is almost stable against the effective saturation for the network with a narrow pore size distribution. This might be the reason that narrow TSP/PSD confines radii in a small range. It can also be seen that the capillary pressure decreases with the increase in standard deviation when the mean radius is fixed. Figure 8b shows that the entry pressure (again, the minimum value of dP c /dS e when S e = 1) is almost the same for the networks with connectivity from three to six. The capillary pressure decreases with increasing pore connectivity. The reason is that a pore connected with more throats offers more selections for air invasion, allowing the capillary pressure to decline. Figures 7 and 8 show the variations of capillary pressure-effective saturation at various morphologies obtained from the pore network modeling (Section 2). Three microstructural parameters of the network are discussed: (1) mean throat/pore radius ( ̅ ) (Figure 7), (2) standard deviation in PSD and TSD ( 0 ) (Figure 8), (3) pore connectivity (Z) (Figure 8). They all affect the capillary-saturation profile of a porous medium in the capillarity-controlled IP. As can be seen in Figure 7, the capillary pressure strongly depends on the mean radius. With increasing mean radius, the value of capillary pressure decreases. This behavior is expected as the capillary pressure is inversely proportional to the mean radius according to the Young-Laplace equation (Equation (4)). Since the capillary pressure in the Brooks and Corey model is presented as a function of effective saturation, the saturation is converted into effective saturation (Se) using

Evaluating the Pore Size Distribution Index (λ) Based on the Capillary Pressure-Effective Saturation
As we could obtain the capillary pressure-saturation profile, the irreducible saturation, and the entry pressure from the pore network simulations, we can evaluate the pore size distribution index by revisiting the Brooks and Corey capillary pressure model. Figure 9 shows that this evaluation results in a nonlinear profile rather than in a specific value. This infers that the pore size distribution index λ depends on the progress of the process (or the effective saturation). Now, it might be questioned whether the profile has a physical meaning. To address this question, first it should be noted that the pore size distribution index λ is a fitting parameter in the Brooks and Corey capillary pressure model. A specific value of this parameter is classically determined by fitting to experimental data of capillary pressure for minimal statistical deviation. Brooks and Corey used this value to express the heterogeneity of a porous medium. However, the primary factor in controlling the capillary pressure is menisci radii distribution. The distribution of the menisci throats/pores radii becomes narrower during, for example, a drying process that follows invasion percolation rules. Additionally, the Brooks and Corey model correlates the capillary pressure with effective saturation. Effective saturation should though be understood as saturation that is controlled by the capillary pressure. In other words, effective saturation in the current work should only indicate saturation in the major clusters that are connected with the bottom reservoir. The isolated clusters and the single menisci should thus not be included in the Brooks and Corey capillary pressure model. However, the smallest radius of throats/pores tends to be located in isolated clusters and single menisci ( Figure 4). This means that the pore size distribution index of the entire network might not be a reasonable quantity for correlating the effective saturation with the capillary pressure.
Since the capillary pressure relies on the radii of menisci, which are still active at a certain effective saturation, an index that would represent the distribution of such menisci would be more accurate for the Brooks and Corey capillary pressure model. Figure 10 is one example of the meniscus size distribution (the normalized meniscus size distribution) during the invasion process. At the initial moment, the standard deviation of menisci equals the standard deviation as assumed for the whole network. The standard deviation of menisci drops dramatically as the large throats are easily invaded, which corresponds to the strong rise of the capillary pressure in the high saturation range (S e ∼ 1). Then, the menisci size distribution becomes gradually narrow with decreasing effective saturation along the invasion process. Finally, a sharp decline in the spread of the meniscus size distribution is observed when the effective saturation approaches zero (S e ∼ 0). Kewen [16] found that the pore size distribution index increases with the decline in the fractal dimension (the heterogeneity) of the network. Figure 9 shows that the index is close to zero in the high saturation range, where the heterogeneity of the menisci is maximum. With ongoing invasion, the distribution of menisci becomes narrower, leading to a decline in the heterogeneity of the menisci. Thus, the index increases slowly during this period. The menisci radii are almost the same when S e ∼ 0, meaning that the heterogeneity of the menisci is then minimum. Figure 9 also demonstrates the index soaring at the lower limit of effective saturation. Therefore, it would be more reasonable to consider the index as a meniscus size distribution index in the Brooks and Corey capillary pressure model. Figure 9 reveals the correlation of the conventionally defined meniscus size distribution index with the structural properties of the porous medium. The meniscus size distribution index is a function of the ratio σ 0 /r and pore connectivity. Figure 9a shows that almost the same profiles are observed for networks with identical σ 0 /r and the same connectivity. The index λ escalates with the decline of σ 0 /r (Figure 9b), which also agrees with Kewen [16]. The index λ is not sensitive upon pore connectivity. It increases when widening the pore connectivity from Z = 3 to Z = 6, but drops again for a pore connectivity of 8. This is similar to the behavior of entry pressure shown in Figure 6. current work should only indicate saturation in the major clusters that are connected with the bottom reservoir. The isolated clusters and the single menisci should thus not be included in the Brooks and Corey capillary pressure model. However, the smallest radius of throats/pores tends to be located in isolated clusters and single menisci (Figure 4). This means that the pore size distribution index of the entire network might not be a reasonable quantity for correlating the effective saturation with the capillary pressure. Since the capillary pressure relies on the radii of menisci, which are still active at a certain effective saturation, an index that would represent the distribution of such menisci would be more accurate for the Brooks and Corey capillary pressure model. Figure 10 is one example of the meniscus size distribution (the normalized meniscus size distribution) during the invasion process. At the initial moment, the standard deviation of menisci equals the standard deviation as assumed for the whole network. The standard deviation of menisci drops dramatically as the large throats are easily invaded, which corresponds to the strong rise of the capillary pressure in the high saturation range (Se ~ 1). Then, the menisci size distribution becomes gradually narrow with decreasing effective saturation along the invasion process. Finally, a sharp decline in the spread of the meniscus size distribution is observed when the effective saturation approaches zero (Se ~ 0). Kewen [16] found that the pore size distribution index increases with the decline in the fractal dimension (the heterogeneity) of the network. Figure 9 shows that the index is close to zero in the high saturation range, where the heterogeneity of the menisci is maximum. With ongoing invasion, the distribution of menisci becomes narrower, leading to a decline in the heterogeneity of the menisci. Thus, the index increases slowly during this period. The menisci radii are almost the same when Se ~ 0, meaning that the heterogeneity of the menisci is then minimum. Figure 9 also demonstrates the index soaring at the lower limit of effective saturation. Therefore, it would be more reasonable to consider the index as a meniscus size distribution index in the Brooks and Corey capillary pressure model. Figure 10. The ratio of standard deviation of the meniscus throat/pore radius distribution (σm) to the given standard deviation of the radius size distribution of the network (σ0) changes with effective network saturation; ( 0 / ̅ = 5/50 (µ m), Z = 6). Figure 9 reveals the correlation of the conventionally defined meniscus size distribution index with the structural properties of the porous medium. The meniscus size distribution index is a function of the ratio 0 / ̅ and pore connectivity. Figure 9a shows that almost the same profiles are observed for networks with identical 0 / ̅ and the same connectivity. The index λ escalates with the Figure 10. The ratio of standard deviation of the meniscus throat/pore radius distribution (σ m ) to the given standard deviation of the radius size distribution of the network (σ 0 ) changes with effective network saturation; (σ 0 /r = 5/50 (µm), Z = 6).

Summary and Conclusions
In this paper, we explored the influence of the morphology and topology of pore network on the macroscopic parameters (wetting phase residual saturation or irreducible saturation S r , entry pressure P e , and the index parameter λ) of the Brooks and Corey capillary pressure model. For this purpose, pore network simulations were carried out for the limiting case of the capillarity-controlled invasion percolation process. This limiting case is relevant in processes such as drainage or slow early-stage drying.
According to our pore network simulation results, the irreducible saturation S r is not a function of the network pore/throat size but depends on the ratio standard deviation/mean radius (σ 0 /r), the ratio pore radius/throat length (r p /L t ), and the pore connectivity (Z). The entry pressure P e depends on the largest meniscus radius in the network. It is close to the mean capillary pressure for a network with narrow size distribution, and constant for a constant value of σ 0 /r. The entry pressure decreases with expanding pore size distribution, as well as at pore connectivity above six. The pore size distribution index λ in the Brooks and Corey capillary pressure model is easier to correlate with the meniscus size distribution profile, which depends on the invasion process and the structural properties. Profiles of the distribution index λ are related to the ratio σ 0 /r. They are the same for networks with the same σ 0 /r and are shifted upwards with the decline of σ 0 /r. Pore connectivity also plays a role in the distribution index. Moreover, the pore network simulations of the capillary pressure have led to results that are compatible with measurements of capillary pressure-effective saturation relationship conducted with water for porous media [39,40].
The capillary pressure profile can be further linked to the transport properties of two-phase flow, i.e., the relative permeability by Darcy's law. Srikiatden and Roberts [23], as well as Joekar-Niasar and Hassanizadeh [41] pointed out that the nonlinear properties, i.e., capillary pressure-saturation and permeability-saturation, are critical for realistic modeling of two (or multi)-phase flow in porous media. Thus, evaluating the transition from the parameter of the capillary pressure-effective saturation profile to the profile of relative permeability over effective saturation could be a future research direction.
One should keep in mind that the present study only considered invasion percolation in the capillary controlled regime. This regime corresponds to situations where capillary pressures are very large compared to the average viscous pressure drops due to the fluid flow. However, the gravity force, liquid viscosity, as well as the secondary capillary structures (rings, films, etc.) can play significant roles in biphasic fluid displacement processes such as drying [40][41][42][43][44]. Additionally, correlating more microstructure aspects of the porous medium with the macroscopic parameters deserves to be considered for further investigation. Moreover, pore network simulations presented in this work deserve to be appraised by experimental results on a real object, e.g., a three-dimensional printed structure.