Abstract
The imbibition experiment is an effective approach for measuring petrophysical properties of porous media, with many such experiments performed over the past decade. Quite some empirical, analytical, and numerical models have been developed to simulate spontaneous imbibition of the wetting phase fluid into porous media, but limitations still exist. In previous studies, the imbibition process has been considered to give a piston-like displacement or the porous medium modeled as multiply-sized pores linked with bonds; both approaches fail to yield comprehensive results due to their neglect of the presence of irregular fractures or nonuniform flow paths through the matrix. By building a numerical model for simulating laboratory-scale experimental data, we performed imbibition tests on several fractured Barnett Shale samples having fractures either parallel () or transverse () to the bedding plane and used MATLAB to build a new numerical model by combining the imbibition process in fractures and the matrix using concepts from percolation theory. The experimental data show that the rocks with -direction fractures have a more steady increase of imbibition rates than the case of -direction one. As the shale matrix with low pore connectivity hampers the upward water movement, the imbibition rate of shales with -direction fractures will decrease suddenly after the bottom layer in contact with water is saturated during the initial period. This wetting phase movement (WPM) model can simulate 3D porous media with 2D fractures. The rate of imbibition by fractured porous media is associated with physical parameters such as porosity and fracture distribution (e.g., the number and angle of fractures). Using Monte Carlo methods, we examined fracture parameters and predicted elapsed time and cumulative water imbibition, for the Barnett Shale samples. The results show that the rate of imbibed water mass is sensitive to the number of fractures directly connected to water source, and the connectivity between two neighboring grid cells is a key parameter for the wetting-front progression. The findings of this study can help to better understand the imbibition process with multiple influencing processes and factors in fractured-matrix rocks. Although the experiments, data simulation, and prediction results are based only on Barnett Shale samples, the model is readily applicable to imbibition tests of other fractured rocks to show the spatial and temporal behavior during a dynamic imbibition process that are not easily captured experimentally.
1. Introduction
Many reservoirs have been utilized in various fields of energy and environmental geosciences, such as groundwater exploitation, safe storage of high-level nuclear waste in geological repository, geothermal mining, and petroleum exploration in both conventional and unconventional rock formations. As these reservoirs form, regularly and irregularly sized fractures at different scales and orientations could be developed, providing preferential flow pathways for transmitting the petroleum and water [1, 2]. Unconventional reservoirs show the significance in petroleum production in recent decades after stimulated fracturing of tight shales [3, 4]. Because of the extremely low permeability and connectivity of shale reservoirs, some field-scale stimulation techniques are necessary to expand the natural, and initiate artificial, fractures, in order to enhance the petroleum production from the shale [5, 6]. As naturally fractured formation with the tight matrix block can cause the low mobility of stored fluids [7], a better understanding of how the fluid flows in fractured low-permeability media has been a challenging problem in the porous media studies. Spontaneous imbibition is a critical process during the hydraulic stimulation and wells shut-in for production preparation, which cause over 50% of the injected fluids to flow through the matrix and fractures by capillary force [8]. This displacement process affects many reservoir properties, for example, the relative permeability and consequent petroleum movement into the producing well [4]. Thus, it is vital to investigate the performance of shales during imbibition process and simulate the fluid flow in shale formations [9].
It is important to clarify the invasion percolation-like motion evolution of the wetting fluid during imbibition into fractured porous media, where the viscous, gravitational, and capillary forces are the main drivers for fluid flow. Dominated by interfacial phenomena, the capillary-driven flow is extremely important in nm-sized pore system, influenced by the flow pathways of connected pore space [10]. The first imbibition model in the petroleum field was developed to predict oil recovery [11]. For a consolidated sand, relative permeability curves showed that, at any given saturation, the permeability to oil during imbibition is smaller than during drainage stages [11, 12]. The theoretical relationship between drainage and imbibition was applied to show relative permeability characteristics but failed to consider heterogeneity of fractures and the corresponding relationship with the matrix.
An analytical solution, including the well-known Lucas-Washburn equation [13, 14], was proposed for capillary rise in tubes when studying imbibition phenomena. Eden’s growth model [15] was proposed to describe the specific types of two-dimensional cluster growth, in which clusters grow stochastically based on accumulation of materials on the boundary instead of the inside. The growth process along with the surface, especially in some natural processes, has inspired many researchers to focus on the growth site, and which particles are the new ones to be added into the whole system. The study of diffusion limited aggregation reflected the stochastic nature of particle walk [16]. A third growth model, invasion percolation, was motivated by fluid displacement under capillary forces in porous media [17]. Besides two-phase flow, three-phase flow in porous media was investigated; a model was built for nonaqueous phase liquid flow instead of allowing all phases to flow, and it only can be applied to two-dimensional rock with uniform pore-throat and -body sizes [18].
In response to the shale revolution since the first decade of the 21st century, spontaneous imbibition into low-permeability shale has been attracting much attention (e.g., [19–22]). Suo et al. [23] presented a two-dimensional model to simulate the imbibition process in heterogeneous porous media, by tracking the moving wetting front with a new interface integral method. Based on Hagen-Poiseuille flow, spontaneous imbibition in tortuous capillaries with variably-shaped apertures was simulated by Cai et al. [19], and their model also considered different geometrical factors as the limitations while those variabilities may cause the changes of imbibed water. Based on the model of Cai et al. [19], the model of Wang et al. [24] for porous media in tight reservoirs calculated the imbibition length, which is hard to measure experimentally; however, as the model took more geometrical parameters into account than other imbibition models, even a small changeable reason may lead to differences in results. However, overparametrization and complex governing equations in a model create a computational burden, especially at large scales of investigation. The rock properties may vary widely even over short distances, such that a model having uniform properties will be unable to fit field data properly, and an overparametrized model will lack flexibility and predictability. An analytical model focused on the imbibition process for confined nanofractures in shales, as Zeng et al. [9] modified the Navier-Stokes equation and considered the whole sample with same and uniform physical properties instead of discriminating the difference between the matrix and microfractures. In microscale studies, the multicomponent flow in heterogeneous porous media was simulated by the modified lattice Boltzmann method [25]. Extending the continuum numerical framework to complex domains is necessary, especially among low-permeability porous media [23].
Several papers reported that a disconnection of the pore network usually occurs only in media with low porosity [26–28], but some occasional reasons may cause the successive probability of water to pass through. Thus, fluid flow and chemical transport in porous media with low connectivity can be described using percolation theory, which contains lattice and bonds stochastically [29, 30]. Furthermore, neutron and X-ray tomography experiments were used to study the fracture development of Callovo Oxfordian claystone and provided a feasible approach to visualizing the development of water-induced fractures during imbibition [31], but recent simulators still lack the visualization of imbibition process in three-dimensional (3D) scale [9, 32–34].
As presented in Figure 1, this work will start with samples and methods for the Barnett Shale and imbibition experiments as well as basic theory and assumptions of modeling concepts and implementation, continue to mathematical modeling for two- and three-dimensions and simulations with fractured Barnett Shale with either - or -direction fractures, proceed to model predictions with cubic and cylindrical sample shapes of the Barnett Shale, and end with summary and conclusions. The wetting phase movement (WPM) model coded with MATLAB addresses the problem by separating the calculation of water imbibition in the matrix and capillary force in fractures. The model also provides a visualization of the imbibition process in 3D scale and successfully predicts the imbibition experiments of Barnett Shale with multisized parallelogram fractures.
2. Samples and Methods
Located in the Bend Arch of the Fort Worth Basin, the Barnett Shale is known as a low-permeability shale gas reservoir. This field covers with plenty of natural gas and oil as a commercial exploration place [35]. Barnett Shale samples were taken from an outcrop located in the Fort Worth Basin, a regional syncline of North Central Texas. The description of Barnett Shale core samples and its properties have been provided by Loucks and Ruppel [36] and Gao and Hu [10]. The samples used for imbibition tests were extracted from those layers parallel () or transverse () to the bedding plane, which are termed as - or -fractured samples.
Firstly, we measured the physical properties of Barnett Shale, such as porosity, permeability, fracture numbers and aperture widths, and air-liquid contact angle. Secondly, we took several pairs of - or -fractured samples to carry out the imbibition tests. The procedures and corresponding data processing of the imbibition test were presented in previous papers [21, 35, 37–39]. The cylindrical samples are around 15 mm in diameter and 20 mm in height. Every imbibition test was not ended until the additional mass of imbibed water can be ignored. Thirdly, after testing several pairs of fractured shale samples, a numerical model, combining the spontaneous water movement in fractures and in the matrix by capillary force, was built as below.
3. Basic Theory and Assumptions
The imbibition process, in which a wetting phase fluid displaces a nonwetting phase fluid, is driven by capillary forces [8, 10, 40]. Figure 2(a) illustrates the distribution of pores and the matrix. For better simulating the water imbibition into rock samples, reconstruction and simplification of the connectivity relationship between pores and fractures are significant to imbibition simulation. Figure 2(b) reflects the diversity of pore sizes and pore throats in 3D fashion for a shale sample, but the simulated results cannot be guaranteed if the computational burden is too heavy for calculation processes. Thus, finding the balance point between reality and efficiency is significant for simulation builders.
At the start of the imbibition test, the rock sample is dried under 60°C for two days, and no water has yet imbibed. With time, water penetrates into the samples though the matrix and fractures, as shown by the growing dark regions. The hydraulic conductivity of fractures is much (~5-6 orders of magnitudes) larger than that of the tight matrix; so, the rate of water imbibition in fractures is much faster. Thus, the pore size diversity of the rock can become a problem for the model. The major complex problem is how to simulate the imbibition process which may be affected by fracture aperture and distribution. This research will investigate the influence of fractures and pore properties on the imbibition process, especially in terms of rate of imbibed water over time.
In this study, the model combines the capillary force and spontaneous imbibition process in both the matrix and fractures. The virtual view of rock sample construction with grid cells is shown in Figure 3. At the start of the imbibition tests, rock samples had an initial water saturation of “zero” (or residual saturation). The rock sample was connected to the electronic balance, with only its bottom face contacting the water inside a closed chamber [35]. The sample was coated on the wall and loosely covered on the top face to minimize vapor absorption inside the chamber full of humid air with water reservoir. The experimental data used in this study were corrected for submergence depth decline and associated buoyancy gain by the method of Hu et al. [21].
The assumptions for the WPM model include the following: (1)For the matrix of the Barnett Shale, the average porosity is around 13%, and permeability ranges from 7 to 50 nano-Darcy. The fracture permeability, obtained from the parallel plate model [42], is six orders of magnitude larger than the matrix. Because the cylindrical samples are small (around 2.0 cm in diameter and 1.5 cm in height), the WPM model considers the matrix to be homogeneous, especially under the condition of much disparate permeability between fractures and the matrix(2)The rock sample can be divided arbitrarily into equal-sized cubic elements (grid cells), with each grid cell assumed to contain the same volume of pores(3)The volume of a fracture depends upon its aperture, width, and length(4)The bonds need to exist to make pore bodies connect with each other. Bond volume can be considered to be part of the pore’s volume; so, it need not be considered separately(5)The aperture widths can be drawn at random from a distribution that is independently derived from thin section petrography(6)The main driving mechanism for wetting-phase imbibition is capillary force. Other driving mechanisms, such as water adsorption by clay minerals [43], electrical double layer expansion [44], and osmotic effect [45], can be ignored(7)During one model time step, the wetting front advances simultaneously into all neighboring connected grid cells. In other words, the wetting front advances one unit of chemical distance [46] in one time step(8)The permeability reflects the rate of imbibition steps in the model; a larger permeability means one step represents longer time(9)Air can be trapped in dangling ends of the pore space during imbibition. Consequently, an individual grid cell may not attain 100% saturation even if water has been flowing through it(10)The wetting front stops moving once it reaches the upper boundary of the sample
4. Mathematical Modeling
The wetting phase movement (WPM) model was derived from the percolation theory which has been applied to water absorption into soils, geological, and other porous media samples [47]. A rock was viewed as a composition of some same-sized grid cells. Every two neighboring grid cells have a connection probability to make fluid go through. Once the pore size or permeability is too small, the corresponding connection probability will be less than the critical value, and the fluid movement will stop. Previous models, such as Hunt et al. [47] and Cai et al. [19], focus on the complex relationships of pores and pore throats; once there is a slight change on the initial values of basic parameters, the simulation results could be very different; thus, the results may be unpredictable and unreliable. In the WPM model, the rock sample was divided into multiple grid cells having identical properties (for example, pore numbers and pore size) to ensure that the moving path has the same probability. The model could embed various-sized or directional fractures to mimic real samples. Each two neighboring grid cells are either connected or not, according to the stochastic probability value. The probability that two neighboring grid cells are connected, termed as connection probability, provides a reliable possibility to simulate the paths that water follows during the imbibition process. In this work, the WPM model was implemented by MATLAB software. First, the model was built on 2D and 3D systems with fractures. Then, the model simulates the experimental data for the Barnett Shale samples. The samples have two kinds of fracture orientation with respect to the upward imbibition direction, parallel, and transverse to the bedding plane; these - and -direction fractures are shown as lines in 2D models and parallelograms in 3D models (see Appendix A).
4.1. 2D Model for Quantifying Water Percolation and Capillary Forces during Imbibition Process
Figure 4 lays out the structure of the 2D simulation models implemented in this work. The parametrization step sets the material properties (water surface tension and contact angle) and decorates the lattice to set fracture locations and apertures and intergrid cell connections. Since the sample bottom is in contact with water, the grid cells of lowest layer have 100% water saturation, barring the small contribution of trapped air. The second step is to begin the percolation process, as the water imbibes from the grid cells of first layer up into the second layer. Each square-shaped grid cell has four directions (up, down, left, and right) for the water molecules to move with an equal probability, but grid cells in the bottom layer have only one direction, which is upward. If one grid cell is crossed by a fracture, the pores in this cell become water-saturated. Then, the water will quickly fill that fracture until it reaches a certain level which is calculated according to Jurin’s law [48]. The neighboring grid cells near the fracture will be filled with water.
In the first step, the sample size and the real fracture position(s) are noted, and all grid cells containing the fractures are identified. Since the bottom layer of rock sample is in contact with water, in the first step, the water saturation of the first grid layer is set to 1 (100% water saturation). Thus, all grid cells in the first layer have no chance for more imbibition, while only the grid cells of second layer can be saturated in the next step. The rate of imbibition is related to the permeability and capillarity of rock sample. Each step represents a period of real time; if the permeability is higher, each step represents longer time. As all models are based on the Barnett Shale, the permeability and corresponding imbibition rates are the same in the simulations, with one step representing one hour. There are two domain theories in this model: the theory of invasion percolation is used within the porous matrix and Jurin’s law in the fractures as shown in Eq. (1) [48]:
As the pore surface of natural rock could be rough, the model adopts a friction coefficient to adjust the capillary rise in the model. Thus, when calculating the real capillary rise height, a parameter, named as friction coefficient (FC), is used in the WPM model; the modified capillary rise is shown in Eq. (2): where is the liquid-air surface tension (force/unit length), θ is the contact angle of water with the solid, ρw is the water density (mass/volume), is the local acceleration due to gravity (length/time2), and is the radius of the pore throat. We used at 20°C, 3, and [33].
Although the fracture volume can be ignored in grid cells the fracture width is significant in the calculation of capillary rise in fractures, the first layer contacts the water source directly; so, the first layer is seen as the free water surface. The model can calculate the capillary rise according to the Jurin’s law. Thus, we can know how high the water will rise in the fracture and in which grid cell it will stop.
Table 1 shows the basic information of a 2D water movement model, which includes the location of fractures and the parameters used for calculating the capillary force. The fracture volume can be neglected when calculating the imbibed volume in the model because the aperture size is much smaller than a grid cell. But in the capillary rise, the aperture width is significant. Thus, the rate of water imbibition will match the observed rate in the experiments.
Figure 5 shows several example steps of the imbibition process in a 2D model, with the blue cells showing the advancement of wetting front. Once a grid cell containing a part of fracture is saturated, the capillary force will make water rapidly fill other parts of the fracture. Water in a matrix cell advances stochastically from any interface cells, but the water can only invade a neighboring cell if the connection between two cells is active. We have assumed that all matrix cells are identical; so, the weight of water imbibition at any time can be obtained simply from the number of saturated cells, volume of a cell, and water density.
(a)
(b)
(c)
4.2. 3D Simulation Model of Water Advance under Capillary Forces
Flow in the 3D simulation model, shown in Figure 6, is much more complex than that in the 2D model, starting with counting how many imbibed cells are crossed by fractures over time; the method of embedding 2D fractures is described in Appendix A. Imbibition in the 3D model will also need to include vapor transport and absorption.
The total volume of imbibed grid cells is the simulation result of mass of water imbibed over time on the current simulation step, which is shown in Eqs. (3) and (4). The ratio of accessible porosity to total porosity (RAPTP) of connected pore volume is a new definition to describe how much of the grid can be filled with imbibing fluid. If a grid cell has several pores and paths linked to neighboring cells, it does not mean that all the connected void space in this grid could be filled with fluid during the imbibition process. RAPTP is an effective parameter to consider the percentage of a grid cell being imbibed by water.
and are the radius and height of the rock sample; is the number of grid cells used in the simulation; is the volume of a single grid cell; is the average porosity of the sample; is the number of grid cells filled with the wetting-phase fluid; is the mass of cumulative imbibed water.
Figure 7 shows the process of imbibition with the red planes representing the fractures and yellow cells crossed by the fractures. The fracture shape is a parallelogram to avoid errors in settling the complex boundaries of fractures and the computational burden in counting the crossed cells of fractures. The model could generate stochastically virtual fractures instead of calculating the positions of fractures.
(a)
(b)
(c)
(d)
The dark blue cells on the bottom layer are in direct contact with a full saturation in step one (Figure 7(a)). The light blue cells on the surface of other five directions contact air directly. These air-contact cells can also absorb water vapor due to the wettability of rock sample, as a thin water layer on the surface. For conveniently observing the invaded cells by vapor absorption, the color of cells invaded by vapor is chosen as 0% transparency. Figures 7(b)–7(d) show the saturated grid distribution of simulation results at steps 2, 4, and 6, respectively.
As shown in Figure 7, all the grid cells belonging to the bottom layer are therefore filled with water and counted into “movement tracing space” (a storage space in the model to track the positions of all imbibed cells), which stores the location of invaded grid cells at each simulation step. Based on chemical distance, the wetting front at the next step advances from the cells which are being filled with water at the current step. If water fails to move to a neighboring cell because the connection probability between those two cells was less than the threshold probability, the model will consider that those two cells are not directly connected. Grid cells that are completely unconnected are impermeable and never attain saturation greater than zero.
If some grid cells are invaded by vapor absorption before imbibition, the imbibed water will fill the rest of void space. The water vapor from the air has an impact on the imbibed mass, with the extent depending upon the rock properties and relative humidity. For example, the rock may absorb more water than normal when the humidity is high or when the rock contains swelling clays [2]. From Figure 7, we can see that the water moves along the fractures. Then, the cells crossed by fractures are filled with water faster than the matrix cells near the fracture cells to be filled with water in the next step by the connection probability by imbibition.
5. Simulations with Fractured Barnett Shale
The WPM model was then used to describe and interpret the laboratory imbibition tests of Barnett Shale. The sample bottom was in contact with water for the imbibition tests of Barnett Shale samples, either with (transverse to the bedding plane) direction fractures or with (parallel to the bedding plane) direction fractures, and the imbibition results are shown in Figure 8. The initial contact of imbibing deionized water causes the steep and noisy phase which should be not considered in the total imbibition process [38]; the WPM model only uses the stable experimental data when the imbibition tests for -and -direction fractures are at 30 and 50 seconds, respectively. The sample mass obtained by the electronic balance over time was greater than actually imbibed mass obtained from the independent weighing. The submergence depth decrease of rock sample in a fluid reservoir from imbibition and evaporation will lead to the reduction of buoyant force, and the corrected data after buoyant change following the approach of Hu et al. [37] was used for the WPM model.
Figure 9 shows four images of fractures with their apertures noted. Most fracture apertures are in the range of 5 to 80 μm, and these values are used in the WPM model.
(a)
(b)
(c)
(d)
There are two main fracture layouts in this study on the Barnett Shale and and directions. Though the shale samples have heterogeneous properties due to many factors, for example, the directions and the distribution of fractures, the connection probabilities in the same layer can be considered to be similar. The probability is larger through vertical axis than horizontal axis in direction sample, while the probability has the opposite performance in direction sample. According to percolation theory [29], there are two kinds of probability in the imbibition process, one between every two neighboring cells, and the other is the description of connectivity of a sample by Monte Carlo methods. When the 1st probability is given, the 2nd probability can be obtained from statistical tests. In other words, the 1st probability is the one for each cell to have active connections (site percolation) or for each connection between neighboring cells to be active (bond percolation), and the 2nd probability is a value that represents the connectivity of the whole sample. Therefore, the 2nd probability is comprised of the value of the 1st probability.
5.1. Dual Matrix-Fracture System with -Direction Fractures
Figure 10(a) shows the sketch of rock surface and virtual fractures. To avoid the uncertainty bias during simulation, the average simulation result is obtained by running the Monte Carlo experiments 1000 times to obtain statistically reliable results, but the running time will be longer than a single run in the same simulation. The laboratory imbibition data after evaporation and buoyancy correction and the average simulated results are compared in Figure 8(b) which shows a close match by the WPM model; the error difference between data and simulation is less than 3% in the long run.
(a)
(b)
To understand the importance of influencing factors, the model was evaluated by various sensitivity analyses (Figure 11). Every simulation line is obtained from the average of 1000 Monte Carlo runs. Figure 11(a) shows the sensitivity analysis of the ratio of accessible porosity to total porosity (RAPTP). The porosity was determined by water immersion porosimetry after vacuum pulling and Archimedes’ method [49] which could measure the volume of all pores connected to the sample’s surface [21]. The accessible pore volume which could be filled with wetting-phase-fluid in a dynamic imbibition test is less than the measured one after full saturation. If there are some pores connected to other large pores, but the pore throat is too small to let water move in, the small pore volume will be considered as nonconnective pores instead of being a part of effective pore volume. It is shown that when the RAPTP is 0.89, the simulation matches the best, with 89% surface-accessible pore space participating the dynamic imbibition process. Though the rates of water imbibed over time with various are almost equal at each step, the minor impact of becomes larger with the increasing imbibition time. The connection probability has the same impact on the mass of water imbibed (shown in Figure 9(b)). The greater value of the probability of imbibition, the more cells are able to be filled with wetting-phase fluid, and the greater value of the mass of water imbibed over time.
There is another important influencing factor during imbibition, which is the water vapor absorption into the sample. The shale sample has been coated with epoxy, a kind of water-proof polymer, on the exterior side. The vapor still could enter the tight shale sample through the top side though it was loosely covered with aluminum foil. Then, the vapor movement is the same as imbibition movement. If a grid cell is permeable, it allows vapor to move through the pores. The probability of vapor movement accessibility also has an impact on the mass of water imbibed over time though the influence that is not as prominent as the imbibition probability (Figure 11(c)). There is an essential parameter to represent the ability of a sample to let water vapor attach on surface and vapor maximum value, shown in Eq. (5). If the vapor enters the interior of a rock, the wetting phase will become a thin layer to attach onto the pore surface. Though the volume of such a layer is small, it still takes a part of pore space. The maximum volume of vapor that could attach onto the pore surfaces depends on the rock wettability, and its influence becomes more prominent over time.
5.2. Dual Matrix-Fracture System with -Direction Fractures
Figure 12(a) shows the sketch of rock surface and virtual fractures, and Figure 12(b) shows the comparison of the laboratory data and simulation results of water imbibition over time. As shown in Figure 12(b), the mass of imbibed water over time has a sudden increase in the first hour. If the fractures are mostly in direction, as settled to be vertical, the water movement may face less restriction than the sample with -direction fractures. When the wetting-phase liquid enters the cells crossed by -direction fractures, the fracture (all the cells at the same level) will be filled with liquid over a short time. In the meanwhile, the upper matrix cells have a relatively lower probability than the fracture cells at the same height level to be filled with liquid because of the gravity. As a result, if the -direction fractures are exposed to the exterior at the bottom layer, those grid cells will be filled with water suddenly because they directly contact the surrounding water phase. During the imbibition process, once the water moves from matrix cells to proceed to fracture cells, the water will fill all cells crossed by fractures and cause a sudden increase of water imbibition.
(a)
(b)
Figure 13(a) shows the sensitivity analysis of ratio of RAPTP in Barnett Shale sample with T-direction fractures. The initial stage is unstable as the spontaneous imbibition is controlled by capillary force [10]. The capillary force is determined by the number of fractures inserted into the cells of bottom layer. To mimic several fractures located at the bottom layer, some fractures shown in Figure 12(a) may be established to contact water at the first step even though their positions are not in the bottom layer. Figure 13(b) shows that the influence of RAPTP; the higher ratio of saturation, the more water being imbibed into the rock. Figure 13(c) indicates that fractures are sensitive to a direct contact of water. The probability of imbibition in Barnett matrix cells is extremely low. Since the vapor could enter small pores easier than liquids, the probability of vapor absorption is much higher than the probability of imbibition. Figure 13(d) shows that the probability of vapor absorption has little effect on the mass of imbibed water at the beginning, but it tends to have an increasing influence over time.
This phenomenon demonstrates that vapor absorption is not the main factor on mass of water imbibed rate at first but exists during the whole experiment. The maximum value of vapor saturation sensitivity analysis shows the same tendency of simulation results which are shown in Figure 13(e). Though it is hard to know the distribution of fractures in the rock sample without X-ray or neutron tomography, the WPM model could well match the mass of water imbibed over time with the laboratory data.
6. Model Prediction
6.1. Predictions for Cubic Barnett Shale Samples
The WPM model is able to predict the imbibition process in the same rock with different geometries. For the simulation accuracy, all the parameters used in the previous benchmark simulation tests (Figures 12 and 13) of the Barnett Shale are applied in the prediction (Figure 14). Though the bulk volumes of the samples are different, the size of the virtual simulation grid cells is set to the same value. For example, if the height of the cylindrical sample is 15 mm while the height of the cubic sample is 10 mm, the simulation grid cells are set as the same value (0.1 mm). Figure 14 shows four cubic samples, which have -direction fractures, -direction fractures, -direction matrix, and -direction matrix. Due to the smaller scale of cubic samples, there are only two visible fractures in -fractured and -fractured rocks (Figures 14(a) and 14(b)), and there is no visible fracture in the matrix-only samples (Figures 14(c) and 14(d)).
Figure 14 shows the comparisons of laboratory data and simulation results. All parameters are almost the same as the parameters used in the simulations of cylindrical samples, which are presented in the sensitivity analyses (Figures 11 and 13).
The prediction results have the same trends with the cylindrical runs. At the beginning of these experiments, the rate of water imbibition is related to the number and position of fractures. After the rapid imbibition stage, the rate of water imbibition tends to rely on the matrix properties, likely the pore size and pore throat in controlling the capillary pressure and permeability of shale matrix.
6.2. Prediction of Cylindrical Barnett Shale
To further assess the accuracy and viability of the WPM model, we tested six pairs of cylindrical Barnett rock samples cored from the same large block, while the dimensions of these samples are shown in Table 2.
Figure 15 shows the comparisons of laboratory data and simulation results for cylindrical Barnett samples. In these figures, sample size is a significant factor to influence the ultimate mass of imbibed water in 24 hours. The sizes of Barnett Shale samples are listed in Table 2. For a given sample size and sampling location, samples with -fractures have a better ability to imbibe water. Even though the -fractured samples are prone to take more water compared to -fractured samples, the imbibition rate in -fractured samples will slow down over time while the rate of -fractured samples tends to be steady.
7. Summary and Conclusions
In this work, we introduce a new model for wetting phase imbibition in the fractured porous media. The significant advantages over previous model are embedding simplified multiple fractures into matrix, calculating the water mass imbibed and visualizing the water movement in every step. The significance between the previous models and current one is that the percolation theory was introduced to simulate the stochastic movement by Monte Carlo approach. Therefore, this new model can test the simulation cases for thousands of times, in order to obtain the most possible results in the simulation. The size, position, and distribution of fractures are based on the laboratory observations and experiments which make the simulation results more reliable. The creative method of emplacing virtual fractures also helps the tracing of water movement paths statistically; thus, the model can have a more reliable prediction of the cumulative imbibed mass. However, as the Barnett Shale samples are reasonably considered in this work as homogeneous in the matrix, other scenarios, such as heterogeneous and irregularly shaped rocks with microfractures and large difference in pore size distribution, may not be predicted by this current model. Overall, from integrated experimental and modeling approaches, this work concludes with the following main advantages and contributions of the new model: (i)The model adopted the aperture width frequency from thin sections, which means that the capillary rise depends on each fracture itself. The WPM model is developed to be able to change the ratio of accessible porosity to total porosity and to simulate the water movement and mass of water imbibed in every step according to the heterogeneous fractures and homogeneous matrix. The new model reflects the performance of percolation process in rocks. The governing equations, e.g., capillary rise, are easily applied into the model to relieve the computational burden. This advantage satisfies the needed of studying complex fractures in rock to ensure the simulation accuracy of imbibed mass(ii)The distribution of fractures has a significant impact on the imbibition process. The multiple crossed fractures could provide effective channels to enlarge the conductivity in the matrix. The water moves easier along with fractures; so, the mass of imbibed water over time is larger in - than -direction rocks. Horizontal fractures prevent water from moving forward, and even though the matrix is homogeneous, the upper grid cells are unable to contact water percolated from the bottom grid cells because of the heterogeneous fractures(iii)The sensitivity analysis shows that the percolation probability determines the rate of water imbibition. The vapor saturation and probability in matrix grid cells affect the value of mass of water being imbibed in every step, especially at the beginning of the imbibition process. The values of vapor saturation and corresponding probability have more influence on the rock samples with -direction fractures than -directions(iv)The reason for the different rates of imbibed water may be due to the fracture distribution. If most fractures are in contact with the bottom layer, like-fractures, where the grid cells have 100% water saturation, the matrix cells have a higher probability to be saturated. With -fractures, the fractures in the bottom layer are saturated, but the water has less probability to move upwards than other directions or to encounter the next fracture. Thus, the steady rate of imbibition may be due to the sample lacking fractures which may cause the acceleration of water imbibition. Overall, the WPM provides capabilities of describing the experimental results and predicting the spatial and temporal distribution of water mass during an imbibition process, which are not easily shown experimentally, for geological samples
Appendix
A. Embedded 2D Fractures in 3D Models
Although natural fractures are irregular in shapes, virtual fractures in the WPM model are assumed to be parallelograms. It is unlikely to have four stochastic points which are located on the same one plane. However, if the coordinates of three stochastic points are known, it is easy to find the fourth point. Thus, the WPM model only needs the coordinates of three known points; the fourth point’s coordinate can be calculated. Assuming that the initial three points are noncollinear, there exists a unique plane which passes through all three points. The equation of that plane is (Eq. (A.1))
The distance between the two parallel lines can be represented as where is the point location of one line; , , , and are the normal vectors of a plane.
One of four boundary lines, which crosses two points (, , ) and (, , ), can be represented as
The point on the boundary line can be written as
The parameter t in Eq. (4) guarantees that the slope of the line has only one direction in , , and axis coordination. Eq. (5) is the transformation of Eq. (4), which guarantees that all the grid cells crossing the four boundary lines could be counted. When the calculation process counts the grid cells crossed by fractures, the model takes the minimum value of points as the value of , , or by the tree traversal method [50]. Since the points (, , ) and (, , ) are given, based on Eq. (5), for any certain value of , we can obtain the formula for the boundary line. Then, we can find all integer points on the line, which are the index of each grid to be either marked as imbibed or not. In the condition of Eq. (5), the model also avoids recognizing the grid cells outside the four boundary lines.
Since the fracture planes are crossed through a number of continuous grid cells, the eight points of cubic grid cells are not typically in the plane of the fracture. For precisely counting all the crossed grid cells, the WPM model adopts the tree traversal method [50]. As the direction and the location of the fractures may vary from sample to sample, making sure that the slope is not infinity is one of the significant procedures in calculation. If there is a fracture plane parallel to one of the coordinate planes, the calculation system will transfer the 3D fractures to 2D fractures for reducing the computing burden. When the length of four lines in one plane has a huge difference, the model will assign a suitable gradient to count every integer-number grid. Only the values of every point in two-dimensional view (such as , , ) are close neighbors, and the third dimensional value of points is continuous. The fracture plane is constrained by a 3D lattice, and every point on the plane has a unique coordinate to determine their locations. The model will judge whether the fracture plane is parallel to any coordination axis or not. If one plane is orthogonal or parallel to plane, the plane will be viewed as a 2D plane. As each individual grid cell is characterized by its (, , ) location, if the coordinates of and are known, then the coordinate can be found by Eq. (2). For example, if there are two points (1, 2, 1) and (1, 2, 4) belonging to the same plane, the intermediate integer points, such as (1, 2, 2) and (1, 2, 3), are counted as the same plane. In this way, the WPM simulator can find all the integer points belonging to one plane and then precisely generate every grid cell crossed by fractures.
This method counts all grid cells between the minimum and maximum values in the same dimension. There are many grid cells in the range of fractures, where their coordinates are between the minimum and maximum values, but they are not in the range of four lines. For example, the grid cells with dark blue outlines in Figure 16(a) and 16(b) are the cells which should be cut off (four points of one fracture are shown in Table 3). Thus, the model not only calculates the normal vectors of every fracture plane to guarantee the points are on the plane surface but also represents four boundary lines as numerical equations. During the traversal process, for example, when the value is established, the line of will cross four points with four boundary lines in axis plane. The space will arrange the four points from the least to greatest values. As the middle two points are the boundary points of the line, the model will only count grid cells precisely between these two end points of a boundary line during the traversal process.
(a)
(b)
Nomenclature
| , , , : | Normal vectors of a fracture plane |
| : | Distance between two points, L |
| : | Capillary rise height, L |
| : | Local acceleration due to gravity, LT-2 |
| : | Radius of pore throat, L |
| : | Radius of rock sample, L |
| : | Height of rock sample, L |
| : | Friction coefficient |
| : | Mass of cumulative imbibed water, M |
| : | Number of grid cells filled with the wetting-phase fluid |
| : | Number of grid cells in a rock sample during the simulation |
| : | Ratio of accessible porosity to total porosity |
| : | Slope of one fracture boundary line |
| : | Maximum ratio of wettability caused by vapor in the air |
| : | Volume of a single grid cell, L |
| , , : | Coordinates of the first point |
| , , : | Coordinates of the second point |
| , , : | Coordinates of any point |
| : | Average porosity of the sample, fraction |
| : | Liquid-air surface tension, MT-2 |
| θ: | Contact angle |
| ρw: | Water density, ML-3 |
Data Availability
The laboratory data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This project was completed with funding provided by the Nuclear Energy University Program, Office of Nuclear Energy, U.S. Department of Energy, award number DE-NE0008797.