Abstract
A systematic numerical method was presented to investigate the effect of aperture distribution on the relation of capillary pressure versus fluid saturation (P-S relation) for a single fracture. The fracture was conceptualized as a two-dimensional lattice-grid model and its aperture field was described by a probability distribution. Based on the invasion percolation theory, a program was developed to simulate the quasi-static displacement. The simulation was verified validly by comparisons of the experimental results. The effects of the statistical parameters were further quantified. The results show that the largest local aperture on the fracture boundary determines the AEV. The larger mean decreases the variation coefficient, which causes the more uniform aperture field, smoother air invasion front, and steeper capillary pressure-saturation curve (CPSC). The larger standard deviation increases not only the range but also the contrast degree of the apertures, thus providing a nondeterministic rule in the P-S relation. The larger correlation length causes a more homogeneous aperture field and a dual connectivity of the fracture. The increase of the difference and contrast degree between the small and large apertures results in dual-aperture fields. The dual-aperture field and dual connectivity of the fracture both contribute to the bimodal characteristic of the CPSC.
1. Introduction
Two-phase flow simulation in fractured media has been of interest over the last five decades in many engineering applications such as petroleum reservoirs [1], geothermal energy production [2], underground disposal of nuclear waste [3], environmental contamination by nonaqueous phase liquids [4], and geological storage of carbon dioxide [5]. The relation of capillary pressure versus fluid saturation is required for these simulations (hereafter referred to as the P-S relation) [6].
Fluid flow in fractured media is often dominated by the highly permeable rock fractures as the rock matrix mostly has a very low permeability. Many studies have analyzed the P-S relation by means of two-phase flow experiments on single fractures. Romm [7] firstly conducted two-phase flow tests on a smooth fracture. In a subsequent study, Persoff and Pruess [8] performed air-water flow experiments on natural rough fractures but without measuring phase saturations because of the measuring difficulties involved. Karpyn et al. [9] used microcomputed tomography to observe the internal structure of the artificially Berea sandstone fracture and in situ water saturations. It was found that phase distribution was mainly determined by the fracture geometry. In a recent study, Babadagli et al. [10] conducted constant-rate immiscible displacement to assess the roughness effect on multiphase flow in a single fracture. The experimental results demonstrated that fracture surface roughness had a critical effect on the hydraulic conductivity and residual displaced phase saturations. Most of the methods are destructive and hence make it difficult to study the effects of aperture distribution on the same fracture. The special fracture replica with adjustable mean aperture was designed to study the P-S relation for different single fractures in our previous experimental work [11].
Although performing experiments is a direct approach, it shows the difficulties in making fracture replicas, controlling capillary pressure, and measuring phase saturations [7–11]. Simultaneously, some researchers tried to build the P-S relation via numerical percolation model. Broadbent and Hammersley [12] were the earliest ones to model capillary flow in fractures using the standard percolation (SP) theory. However, Pruess and Tsang [13] did not observe the residual wetting phase saturation after using the SP model. Although Kwicklis and Healy [14] imposed an accessibility criterion to consider the continuous connecting path, the authors did not obtain the residual wetting phase saturation either. Wilkinson and Willemsen [15] introduced an invasion percolation (IP) model to incorporate the wetting fluid trapped in a fracture in the early time. Amundsen et al. [16] used IP model to compare simulations with two-phase flow experiments in the fractures oriented horizontally and found that the IP model could predict the slow, immiscible displacement processes in narrow fracture accurately. However, Pyrak-Nolte et al. [17] indicated that the IP model would lead to a rough invading phase and high residual saturation. To include the effect of local aperture geometry and in-interface curvature, Glass et al. [18] proposed a modified invasion percolation (MIP) model. Yang et al. [19] also pointed out the importance of in-plane curvature in accounting for capillary pressure in many cases. However, significant similarities could still be observed through the comparison of P-S relation obtained from the simulations based on the IP and MIP models, respectively [20]. Likewise, Li et al. [21] suggested that the MIP model produced little improvement over the IP model owing to the high value of the curvature radius of the invading front.
Natural fractures always have spatially variable apertures [22]. Weerakone et al. [23] found that the apertures of rough-walled fractures not only change randomly but also show spatial correlation. It is desirable to describe the aperture field of a single fracture using a probability distribution and take the spatial correlation of the apertures into account. In addition, the experimental and numerical results both imply that the P-S relation is influenced significantly by the aperture distribution. Neuweiler et al. [24] theoretically investigated the effect of correlation structure of the aperture field on characteristic properties of the immiscible displacement of water by air. Yang et al. [25] indicated the larger standard deviation of apertures would give rise to smaller air-entry pressure and larger residual water saturation by numerical simulation. Moreover, Ye et al. [26] used the fractional Brownian motion to model the rough fracture surfaces and employed the IP approach to model the two-phase displacement. The simulations verified that fracture properties were closely related to residual saturations depending on spatial correlation of aperture distribution. Our previous experimental studies also revealed the mean aperture of the single fracture has a strong influence on the P-S relation. Nevertheless, it is obvious that the systematic research on the effect of aperture distribution on the P-S relation for a single fracture in existing works is still insufficient. Meanwhile, it is crucial to have insight numerical analysis on the effect of aperture distribution in view of the difficulty of the experimental research. Therefore, the objectives of this study are to (i) generate the single fracture with spatially variable apertures; (ii) build the P-S relation based on the invasion percolation theory; and (iii) quantify the effect of aperture distribution on the P-S relation.
2. Generate the Aperture Field of the Fracture
2.1. Single Fracture Model
The fracture was 100 mm × 100 mm in shape and conceptualized as a two-dimensional lattice-grid model with 10000 discrete square elements which were arranged in 100 rows and 100 columns. Each element was 1 mm × 1 mm in and directions on the microscopic scale, respectively, and behaved like a small parallel-plate having a constant aperture. The illustration of the fracture model is shown in Figure 1.
Generating the aperture field of the single fracture needs the knowledge of the aperture distribution. Advanced experimental techniques have been used to characterize the void geometry of single fractures. The aperture distributions were often found to follow lognormal distributions approximately in most cases excluding some contact points [9, 27]. Therefore, the apertures of the fracture were treated as log-normally distributed in this research. The spatial fluctuations of the apertures were described by the correlation length that showed the correlation degree between apertures at different locations. By definition, the local aperture values are likely to be similar if the distances of two elements in the fracture plate are smaller than the correlation length.
2.2. Random Aperture Field
The aperture field of the fracture was generated by assigning spatially correlated and log-normally distributed random numbers to the local apertures of the elements. The generation process of the desired random numbers is shown in Figure 2, in which is the coordinate location of the element mark as , , , and are the mean, standard deviation, and correlation length of the apertures, respectively.
The assembly of , , and together determined the covariance matrix , where is the separation lag. The Cholesky technique was then employed to decompose into a lower triangular matrix . On the other hand, the multiplicative congruential generator (MCG) was used to yield two sets of uniformly distributed random numbers, that is, and that were located in the interval of . Then the polar algorithm was employed to transform the forgoing numbers into pairs of random numbers, that is, and that obeyed the standard normal distribution with a mean of zero and standard deviation of one. and sequentially formed a random vector . The normally distributed numbers were calculated from the , , and ultimately. The log-normally distributed values were obtained by the natural exponential function of .
The realization of the generation process is available in COVAR program [28]. The program requires a set of two-dimensional coordinate locations and three main input parameters: mean, standard deviation, and spatial correlation length. To compare the simulations with the experiments, the fracture aperture field with mm, mm, and mm was firstly used, which were identical to the previous experimental settings. Meanwhile, the base setting of the statistical parameters justified the approximation of ignoring the influence of the in-plane curvature according to the discussion by Glass et al. [18]. The generated aperture field with a resolution of elements is depicted in Figure 3.
3. Build the P-S Relation for the Fracture
3.1. Description of the Previous Experiment
A laboratory device was designed to measure the P-S relation. The apparatus includes a rough-walled fracture, two end caps, and a water content monitoring and pressure control system (see Figure 4).
The fracture consists of two transparent epoxy plates of equal dimensions (150 mm long and 100 mm wide) that are clamped together. The air-water contact angle of this epoxy plate is about 53° indicating the epoxy is a hydrophilic material, which means the fracture walls are similar to hydrophilic rocks. The upper plate is smooth and the bottom one is equipped with the surface roughness accurately controlled by machining using a high-resolution digital lathe. The fracture aperture size changes spatially and follows a lognormal distribution. The effects of the mean aperture sizes on the P-S relation were investigated by adding the thin steel plates of different thickness between the two plates.
The fracture was saturated first by introducing deaired water through the water content monitoring system (see Figure 5). After the fracture was saturated, the water level in the horizontal tube was adjusted down, 1 mm at each step to decrease the water pressure, 0.01 kPa per step. The special end caps were fabricated to introduce air into the fracture uniformly and measure water volume changes accurately. The air pressure in the fracture was maintained at atmospheric pressure. The capillary pressure applied to the fracture was then equal to the difference between the air pressure and the water pressure. The drainage water volume from the fracture was used to calculate the saturation degree of the fracture. Finally, the capillary pressures and degree of saturations at different equilibrium steps were plotted as a capillary pressure-saturation curve (CPSC). In addition, the two-phase flow was observed visually through the upper transparent epoxy plate. Several snapshots were taken to show the process of water flow out the fracture by a digital camera.
The water flow was found to be mainly controlled by the capillary force considering the small aperture size and the horizontal arrangement of the fracture. The invasion percolation theory was used to simulate the quasi-static displacement process to build the P-S relation for the single fracture under drainage condition.
3.2. Definition of Capillary Pressure
Air and water were chosen as the “invader” and “defender,” respectively, and the water was considered as incompressible. As shown in Figure 6, the single fracture was placed horizontally and fully water saturated first. The top and bottom edges of the fracture plane were considered to be no-flow boundaries. The simulation began by air invading elements on the left boundary of the fracture plane, and the water was displaced from the right outlet. It is noted that water can also flow out the left boundary as long as the air has no hindering effect.
The fracture is assumed to be square lattice of pores in our simulation. When modeling the quasi-static displacement within a horizontal, variable aperture fracture, the physical law that controls the air-water displacement is the Young-Laplace equation which provides the local capillary pressure () for a set of principal radii of curvature of the air-water interface ( and ) [29]:where is the surface tension. represents the principal radius of curvature perpendicular to the fracture plane which is dependent on the local aperture and the contact angle of the fluid associated with local convergence or divergence angle of the fracture wall, and it can be expressed as is the radius of in-plane curvature, which is determined as a function of the spatial correlation length and the included angle . Then, it can be obtained as
The effect of in-plane curvature can be neglected in comparison to that of when the spatial correlation length is much larger than the mean of apertures. For the sake of simplicity, the local convergence or divergence angle can be ignored, which reflects the topology of the fracture surface. Thus the sufficient capillary pressure required for air to invade a parallel-plate element is given by
3.3. Calculation of Water Saturation
The key issue of building the P-S relation at a particular capillary pressure is to determine the elements that are filled with water. In consideration of the small values of the element area and local aperture, the elements were assumed to be only filled with water or air. A global aperture threshold can be deduced from the inverse of the capillary pressure:
Equation (5) is the first criterion to satisfy for air to invade an element based on the SP theory which is also called the capillary pressure criterion. Each element with the local aperture larger than the aperture threshold is occupied by air, and the remaining elements are occupied with water. The second requirement adopted is the air accessibility criterion; that is, the air must have continuous paths to invade the elements. If the elements with larger apertures occupied by air from the judgment of the capillary pressure criterion disconnect to an external air flow path, it will still remain water-filled. The last one is water trapping criterion; that is, the local apertures can be air-occupied only if there is an escape route for the water being displaced. In an invasion percolation process with trapping, water clusters can no longer be invaded if they are surrounded by air-occupied elements and cut off from the water outlet of the fracture plane. After the elements filled with water are confirmed, the water saturation can be defined as follows:where is the total volume of the elements filled with water, is the total volume of entire elements, is the number of the elements filled with water, is the number of entire elements, is the aperture of the element denoted as , and is the side length of the square element.
The quasi-static displacement was simulated to calculate the water saturation automatically by the author edited code based on Fortran 90. As the setting of input statistical parameters for random generator can create a multirealization of aperture field, each realization was obtained by changing the seed number of the generator. In order to build the ensemble averaged P-S relation, the average of ten results was taken. The environmental temperature was set as 25°C and the surface tension was about 71.79 mN/m [30]. In consideration of the hysteresis effect of contact angle, the drainage contact angle of the hydrophilic epoxy plates was assumed to be 35°. The capillary pressure was increased with a small interval of 0.01 kPa at each step. All of these settings agreed with the previous experiment.
The capillary pressure and the corresponding water saturation at different steps were plotted as the capillary pressure-saturation curve (hereafter referred to as the CPSC) which is often employed to reflect the P-S relation for the fracture as shown in Figure 7. In addition, the development of air occupancy patterns was recorded simultaneously by the water saturation map constructed for each step of the simulation.
3.4. Results and Analysis
As illustrated in Figure 7, the CPSC is “S” in shape and can be divided into three stages: water saturated stage, water rapid-losing stage, and water residual stage. There are three key parameters determining the shape of the CPSC including the air-entry value (AEV), the slope (SL) of the central section of the curve, and the residual water saturation (). As the central section of the curve is like a straight line, SL can be simply given bywhere is the residual capillary pressure corresponding to the . The parameter SL indicates the average water saturation reduction in response to a change in capillary pressure.
A comparison of experimental results and numerical data is also shown in Figure 7. The AEV was 0.23 kPa, which was equivalent to the measured one. SL and were 4.4 and 0.3, which were slightly smaller and larger than the experimental values of 6.6 and 0.26, respectively. This difference is caused by the air-water flow simulation and does not take the film flow of the water along the fracture surface into account. Another reason is that the invasion percolation model yields a larger residual water saturation without the implementation of the in-plane curvature. Though the two curves are not an exact match, the numerical approach shows approximately similar behavior and properties of the experimental results. The agreement validates the numerical method based on the invasion percolation theory.
Figure 8 shows a sequence of the water saturation maps of one realization and the snapshot of the residual water of the previous experiment. The white squares indicate the elements filled with air, while the black squares indicate the elements filled with water. Figure 8(a) shows the fracture stayed in the water saturated stage before the average capillary pressure exceeded the AEV. Figure 8(b) shows the air firstly intruded into the element with the largest local aperture on the fracture boundary. The phenomenon can be called the boundary effect since the local aperture on the boundary determines the value of AEV. When the capillary pressure continued to increase, more air came into the fracture. The water in all elements was initially continuous; that is, it connected with each other and extended to the fracture plane boundary. Soon, however, the water turned into being discontinuous, and a cluster of residual water was trapped in the elements near the fracture plane boundary as shown in Figure 8(c). The change is a sign of CPSC entering the water rapid-losing stage. The water saturation rapidly decreased with increasing the capillary pressure as illustrated in Figure 8(d). Air invaded the elements along the route with the least capillary barrier quickly. The air invasion route was always continuous but the water in the elements with small apertures was kept and lost outflow path to the fracture boundary gradually. Figure 8(e) shows that more elements were occupied by air. The swept area of air was greater with the increase of capillary pressure. After the water lost all the continuous routes, there was no water flow out the fracture even at higher capillary pressure, and the trapped water was left as residual water clusters. Owing to the randomness of the aperture field, the residual clusters were fractals scattered in the fracture and had nonuniform sizes. Figure 8(f) indicates that the experimental pattern had similar features. It is verified again that the numerical approach proposed is capable of building the P-S relation for a single fracture.
(a)
(b)
(c)
(d)
(e)
(f)
The water saturation is often expressed as a function of the capillary pressure in a porous medium by the van Genuchten (VG) equation [31]:where , , and are the fitting parameters. The numerical data were fitted to the van Genuchten function. It can be seen from Figure 7 that the VG function described the predicted data well. The best values of the fitting parameter were 0.27 for , 3.27 for , and 12.95 for and the Adj. -Square is 0.998. It demonstrates that a single fracture can be modeled as a two-dimensional porous medium since the P-S relation for the single fracture can be expressed well as a VG equation. By computation, the parameter was found to be approximately equal to the reciprocal of the AEV. Meanwhile, the exponential coefficient allows for rapid water loss from the fracture once the AEV is reached and further affects the value of SL. Since the fitting parameters entirely determine the VG function, the AEV, SL, and can be the representative indicators to quantify the effect of aperture distribution on the P-S relation for the fracture.
4. Parametric Effect Study and Discussions
The capillary behavior of two-phase flow depends on the contact angle and the pore size distribution in porous media. Analogously, the contact angle and the aperture distribution have an important influence on the P-S relation for the single fracture with variable apertures. Since the drying process is mainly considered, the contact angle hysteresis effect is not in the scope of discussion. The small drainage contact angle for the hydrophilic rock in the drying process is assumed as 0° for the convenience of calculations. The capillary pressure is applied to the fracture at an accuracy of 0.05 kPa since the narrow fracture is harder to lose water. The effects of the statistical parameters, including mean, standard deviation, and correlation length, are discussed in the following study.
4.1. Effect of Mean
The means of apertures were set at 0.15 mm, 0.20 mm, 0.30 mm, 0.40 mm, and 0.50 mm. All simulations had the same standard deviation of 0.10 mm and correlation length of 0.36 dm. The values concur well with the measurement of nature dolomitic limestone fractures [32]. The fracture is narrow to ignore the influence of film flow and in-plane curvature. The simulated capillary pressure-saturation curves for the fractures are shown in Figure 9.
As shown in Figure 9, the AEVs were 0.50 kPa, 0.45 kPa, 0.35 kPa, 0.25 kPa, and 0.25 kPa for the five fractures. It demonstrates that the AEV for the fracture decreases with the larger mean. The values of SL changed from 0.62, 1.05, 2.46, and 2.90 to 5.93 and it is manifested as the curves becoming steeper when the mean aperture increases. The values of were 0.22, 0.21, 0.14, 0.13, and 0.11 for the five simulations. It indicates that for a fracture is smaller with a larger mean aperture. The general trend of the simulations was accordant with the experimental results. Particularly, the CPSC was found to have a characteristic of the bimodal curve when the mean aperture was the minimum value of 0.15 mm. The CPSC was divided into two segments by the vertical dashed line. The left segment was like a reduced CPSC when the capillary pressure was less than 0.8 kPa, while the right one was also like a normal CPSC when the capillary pressure exceeded 0.8 kPa.
From the knowledge of (4) and the boundary effect, the AEV is inversely proportional to the maximum aperture on the boundary of the fracture. The maximum aperture increases with increasing the mean and hence the value of AEV decreases. Moreover, the variation coefficient (i.e., the ratio of the standard deviation to the mean: ) turns smaller when the mean becomes larger, but keeping the standard deviation of the apertures unchanged. The smaller variation coefficient gives rise to a relatively more uniform aperture field, which causes a larger value of SL.
As shown in Figure 10, the swept area of air increased gradually in the same capillary pressure interval as the mean was increasing. The increasing swept area connected to each other and led to a massive invasion region. Meanwhile, the air invasion front became smoother in a more uniform aperture field and resulted in a steeper CPSC. The amount of water held in the elements with small apertures decreased when the mean was increasing. It is reasonable because air can only occupy the elements in which the aperture is larger than a critical aperture. Figure 10(a) supports a remarkable behavior that air quickly invaded the elements on the top left corner of the fracture after exceeding the first AEV of 0.45 kPa and afterward, the first residual capillary pressure of 0.6 kPa was also fast reached. The behavior of rapid air invasion reappeared when the capillary pressure exceeded the second AEV of 0.95 kPa. These behaviors imply that a small mean and a large variation coefficient tend to increase the difference between the small apertures and large apertures, which results in a dual-aperture field to form the bimodal curve.
(a)
(b)
(c)
4.2. Effect of Standard Deviation
To investigate the P-S relation for the fractures with different standard deviations, the standard deviations were varied from 0.01 mm, 0.04 mm, 0.10 mm, and 0.14 mm to 0.24 mm, while the mean and correlation length were fixed at 0.30 mm and 0.36 dm, respectively. The capillary pressure-saturation curves for five different cases are shown in Figure 11.
There is not a deterministic rule between the P-S relation and the standard deviation from the observation of Figure 11. Figure 12 shows the summary of the detailed values of AEV, SL, and . It is shown that increasing standard deviation definitely caused smaller AEV, while there was no definitive correlation between the SL and the standard deviation. When the standard deviation was smaller than 0.14 mm, the SL had a negative correlation with the standard deviation. Meanwhile, when the standard deviation exceeded 0.14 mm, SL had a positive correlation with the standard deviation. The opposite feature of correlation with standard deviation existed in and SL. When the standard deviation was smaller than 0.10 mm, positively correlated with the standard deviation. Meanwhile, when the standard deviation exceeded 0.10 mm, negatively correlated with the standard deviation. It is notable that the CPSC for the fracture with the minimum standard deviation was of the steepest curve. The corresponding P-S relation seemed like an on-off function with 19.2 for SL and 0.04 for . Interestingly, the fractures with a larger variation coefficient had the same characteristic of the bimodal curve.
The unexpected results can be inferred from the function of the probability density (pdf) and the apertures under different standard deviations as shown in Figure 13. It can be observed that the pdf curve turned to be flat and got both larger and smaller apertures when standard deviation increased. Although the lognormal distribution is biased towards smaller apertures, taking the boundary effect into account, it is easier for the air to invade the elements with larger apertures on the boundary, thus causing the smaller AEV. When the standard deviation was the minimum value of 0.01 mm, the apertures lied in a narrow range. An extreme case is that when the standard deviation is equal to 0, the fracture will be a smooth parallel-plate with a constant aperture.
Figure 14(a) proves that there were only two swept areas in the fluid distribution map because the fracture was smooth. On the contrary, increasing the standard deviation caused a larger variation coefficient, which was the reason for the aperture field becoming nonuniform. As shown in Figure 14(b), the air invasion front became rough and a higher amount of water was trapped in the fracture and hence the SL decreased and increased. When standard deviation increased further, the SL and had a departure from the original correlation. The air swept area became larger in the same capillary pressure interval. Figure 14(c) shows the air invaded the fracture elements to form a vast region, which was about half area of the fracture plane when the capillary pressure reached the maximum value. The phenomenon also reveals a special on-off relationship between the capillary pressure and water saturation. Because the log-normally distributed aperture field tends to obtain smaller apertures inherently, increasing standard deviation brings about a higher probability of having small apertures but lower probability of large apertures. That is, increasing standard deviation increases not only the range but also the contrast degree between the small and large apertures. When the standard deviation was increased to 0.24 mm, the corresponding aperture field was like a special parallel-plate which had overwhelming majority value of apertures less than the mean aperture but a fraction of large apertures. This resulted in the increase of SL and decrease of when the value of standard deviation was increased from 0.14 mm to 0.24 mm. In addition, the increasing contrast degree led to a dual-aperture field, which resulted in the bimodal curve.
(a)
(b)
(c)
4.3. Effect of Spatial Correlation Length
The analysis of the effect of correlation length on P-S relation was simulated for fractures under cases with different correlation lengths ranging from 0.08 dm, 0.22 dm, 0.36 dm, and 0.50 dm to 0.64 dm but keeping the mean and standard deviation constant, which were 0.20 mm and 0.10 mm, respectively. Figure 15 plots the simulated capillary pressure-saturation curves for the fractures.
As shown in Figure 15, the AEV and for a fracture both decreased slightly with the increase of correlation length. The predicted SL had no dependable correlation with the change of correlation length. It fluctuated approximately around the value of 1.27. Although there were only some subtle differences between the curves, the effect of correlation length can be clearly seen from the aperture fields with different correlation lengths shown in Figure 16.
(a)
(b)
(c)
The realization of the aperture field with dm shown in Figure 16(a) had the most random distribution. The field looked heterogeneous as the small apertures and larger apertures randomly scatter in the field. In contrast to the case of small correlation length, it had a higher probability for the local apertures to have similar apertures at neighboring points with a larger correlation length. Therefore, the aperture field was more homogeneous when the correlation length became larger as shown in Figure 16(b). This results in larger apertures on the boundary of the fracture plane and then leads to a smaller AEV. When the correlation length was increased to 0.64 dm, the field showed an isolated area with large apertures in the upper part of the field, while the area with small apertures was observed in the lower part of the field as shown in Figure 16(c).
By comparing Figures 17(a) with 17(b), it is found that the significant difference between the aperture fields controls the locations of the air invading pathways; thus it can be concluded that the larger the correlation length is, the easier the continuous air invading pathways are to form. It is also the reason for smaller values of for the fractures. It is noted that the air invading pathways are isolated from the others although they are continuous under the large correlation length. Figure 17(c) shows a dual connectivity of the fracture that air prefers to invade the elements with strongly correlated apertures and avoid invading the ones with weakly correlated apertures, and therewith the curves have the characteristic of a bimodal curve especially when the correlation length is large.
(a)
(b)
(c)
5. Conclusions
A set of numerical methods was presented to study the effect of aperture distribution on the P-S relation for the single fracture. The aperture field of the fracture was generated using a probability distribution method. The CPSC for the fracture was built based on the invasion percolation theory. The simulations were similar to the previous experimental results. The effects of mean aperture, standard deviation, and correlation length on the P-S relation were then investigated, respectively. The following conclusions are drawn from the study:
() The CPSC for the single fracture is fitted well using the VG equation, which confirms that the single fracture can be viewed as a porous medium. The AEV, SL, and can be the indicators to quantify the effect of aperture distribution on the P-S relation. The AEV is controlled by the boundary effect and found to be approximately equal to the reciprocal of the fitting parameter . The SL can reflect the slope of the central section of the CPSC like the fitting parameter . When the applied capillary pressure becomes high, the water phase will turn to be discontinuous. can be calculated from the volume of residual water clusters scattered in the fracture.
() The larger mean of apertures causes the smaller AEV, larger SL, and smaller . The smaller variation coefficient of apertures leads to a relatively more uniform aperture field, which results in the smoother air invasion front and steeper CPSC. The increasing standard deviation increases the range of apertures, and it causes the smaller AEV. The SL and have no definitive correlation with the standard deviation. When the standard deviation is relatively small, the SL decreases and increases with the increase of the standard deviation. The fracture with the minimum variation coefficient behaves like a smooth parallel-plate. The corresponding P-S relation seems like an on-off function. The increasing standard deviation also increases the contrast degree between the small and large apertures. Therefore, the SL increases and decreases with the further increase of the standard deviation.
() The values of AEV, SL, and only have subtle fluctuations under different correlation lengths, but the significant differences between the aperture fields are found. The larger correlation length gives rise to a more homogeneous aperture field, which forms a dual connectivity of the fracture. The dual connectivity for air to invade contributes to the bimodal characteristic of the CPSC when the correlation length is large. In addition, the decrease of mean and increase of standard deviation both lead to a dual-aperture field. The dual-aperture field also results in the same characteristic of the bimodal curve.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research was supported by the Fundamental Research Funds for the Central Universities (Project no. 2015B36114).