Dynamic Pore-Scale Network Modeling of Spontaneous Water Imbibition in Shale and Tight Reservoirs

: Spontaneous water imbibition plays an imperative role in the development of shale or tight oil reservoirs. Spontaneous water imbibition is helpful in the extraction of crude oil from the matrix, although it decreases the relative permeability of the hydrocarbon phase dramatically. The dynamic pore-scale network modeling of water imbibition in shale and tight reservoirs is presented in this work; pore network generation, local capillary pressure function, conductance calculation and boundary conditions for imbibition are all presented in detail in this paper. The pore network is generated based on the characteristics of Barnett shale formations, and the corresponding laboratory imbibition experiments are matched using this established dynamic pore network model. The e ﬀ ects of the wettability, throat aspect ratio, viscosity, shape factor, micro-fractures, etc. are all investigated in this work. Attempts are made to investigate the water imbibition mechanisms from a micro-scale perspective. According to the simulated results, wettability dominates the imbibition characteristics. Besides this, the viscous e ﬀ ects including viscosity, initial capillary pressure and micro-fractures increase the imbibition rate, while the ﬁnal recovery factor is more controlled by the capillarity e ﬀ ect including the cross-area shape factor, contact angle and the average pore-throat aspect ratio. With the existence of micro-fractures, the imbibition rate increases signiﬁcantly while the ﬁnal recovery remains similar. The micro-factures contribute to the conductance of the pore network instead of the pore volume, and the increase of conductance increases the imbibition rate but not the ﬁnal recovery, as the ﬁnal recovery is controlled by the capillarity of the matrix pores and throats. found within the ﬁtting process. Since micro-fractures are observed in the imbibition experiments, the characteristic length is reduced signiﬁcantly. In the ﬁtting results, the e ﬀ ective characteristic length of Barnett shale samples is approximately several millimeters, and the fraction of oil-wet pores ranges from 70% to 50%.


Introduction
The spontaneous water imbibition phenomenon is commonly encountered in the development of shale or tight reservoirs. Water imbibition is helpful for the replacement of oil from the matrix and increases the oil recovery, while it decreases the relative permeability of the hydrocarbon phase dramatically, which is detrimental to the production rate and recovery factor [1,2]. Generally, the capillarity is believed to be the controlling factor for the imbibition of water in shale and tight formations. Analogous to waterflooded carbonate reservoirs [3][4][5], taking full advantage of water imbibition to yield the highest oil recovery is a major concern for the development of shale and tight reservoirs. Moreover, surfactant additives are proposed to alter the formation wettability and improve the performance of water imbibition [6][7][8]. Large numbers of experimental studies of spontaneous imbibition have been reported in the literature for both conventional [9][10][11][12][13] and unconventional reservoirs [14][15][16]. Large numbers of analytical models based on a bundle of capillary tubes have been proposed following the early works by Lucas [17]. Viscosity ratios, the tortuosity of capillary tubes, the shapes of cross-sectional areas and the assumption of fractal porous media are involved in these models [18,19]. An attempt has been made to find a universal upscaling equation since the imbibition for the studied porous media before simulating two-phase flow processes within them. A structured network consisting of nodes and connecting bonds is used, where nodes are assigned as pore bodies and bonds as pore-throats. A truncated lognormal density distribution function (c.f. Equation (1)) commonly yields a good fit for the profiles of the pore-throat distribution of reservoir rocks [38]. Since the larger pores tend to be connected with larger throats, we cannot randomly assign the pore bodies connected to the pore-throats; here, we sort the pore bodies by the mean values of the connected throats' radii, and then the corresponding sorted pore body radii are assigned.
where, r refers to the radius of pores or throats, and r m , r max , r min refer to the mean, maximum and minimum radii of pores or throats, respectively; σ is the standard derivation of the radius distribution. Following the work by Mason and Morrow [40], in which they studied the capillarity behavior of drainage and imbibition in irregular perfect wetting triangular micro-tubes, the definition of the shape factor is used as shown in Equation (2).
where A and P are the cross-sectional area and the perimeter of a pore or throat. For a cross-sectional shape, if we define the corner half angles of the triangle as β 1 < β 2 < β 3 , two constraints need to be met, given a shape factor value: where β 2 is randomly chosen to range from the lower and upper boundaries in Equations (5) and (6).
The mean values of Equations (5) and (6) are used; thus, the smallest half corner angle can be calculated by Finally, β 3 is calculated by applying Equation (3). In this way, we can use a single shape factor for the whole pore network. Therefore, one single shape factor will help us to use a single dimensionless local capillary function for every pore or throat with variable wetting phase saturation, which will be presented later. The schematic of the pore network and the triangular cross-section of every pore element are presented in Figure 1.

Control Equations and Conductance Calculation
The proposed dynamic pore-scale network model is based on the following assumptions: (1) The volume of throats is negligible compared with that of pore bodies, while the throats control the conductances between pore bodies; (2) the pressure field is continuous, and an upwind scheme is used for the numerical treatment, while the saturation could be discontinuous, especially for the connecting throats; (3) fluids are assumed to be incompressible, and the pore network is rigid.
Applying the mass conservation for every pore body, the governing equations are obtained as where , , and refer to the volume, phase saturation and pressure for every pore body, respectively, and and are the absolute and relative conductance between the and pores.
The geometry of pores or throats controls the conductance of multiphase flow. For single-phase flow in an irregular triangle, the absolute conductance is obtained by [28] = 0.6 (9) where is the length of a pore or throat. For the two-phase flow scenario of the triangular cross-section, the non-wetting phase occupies the center, while the wetting phase resides in corners, which leads to extra resistances for wetting phase flow [41]. Following the study by Sheng and Thompson [29], the relative conductances for two-phase flow in a throat during the imbibition process are further derived in Equation (10a,b):

Control Equations and Conductance Calculation
The proposed dynamic pore-scale network model is based on the following assumptions: (1) The volume of throats is negligible compared with that of pore bodies, while the throats control the conductances between pore bodies; (2) the pressure field is continuous, and an upwind scheme is used for the numerical treatment, while the saturation could be discontinuous, especially for the connecting throats; (3) fluids are assumed to be incompressible, and the pore network is rigid.
Applying the mass conservation for every pore body, the governing equations are obtained as where V, S, and p refer to the volume, phase saturation and pressure for every pore body, respectively, and K ij and K r ij are the absolute and relative conductance between the i th and j th pores. The geometry of pores or throats controls the conductance of multiphase flow. For single-phase flow in an irregular triangle, the absolute conductance is obtained by [28] where L is the length of a pore or throat. For the two-phase flow scenario of the triangular cross-section, the non-wetting phase occupies the center, while the wetting phase resides in corners, which leads to extra resistances for wetting phase flow [41]. Following the study by Sheng and Thompson [29], the relative conductances for two-phase flow in a throat during the imbibition process are further derived in Equation (10a,b): where C(G, θ) is the function of the shape factor and contact angle, accounting for the extra resistance in corners of the wetting phase, and S w and S n refer to the saturations of wetting and non-wetting phases, respectively; F(G, θ) will be introduced later in Equation (13). Following the analytical derivation of Øren et al. [21] and Valvatne and Blunt [22], the expression of C(G, θ) is rearranged in Equation (11a-c). Note that C(G, θ) does not depend on the local capillary pressure and saturation, which is only related to the cross-sectional shape factor and contact angle. Given a shape factor and a contact angle, a constant value of C(G, θ) is calculated and will be used in Equation (10b).

Local Capillary Pressure Function
For every pore element (pore body or throat), the fluid distribution is controlled by the local capillary pressure function, which is dependent on the element radius, shape of the cross area and contact angle. Based on the previous studies [27,42,43], the local capillary pressure function is reformulated as shown in Equation (12): However, Equation (12) is only partially applicable, as either snap-off or piston-like-filling happens before the wetting phase saturation reaches unity. The snap-off occurs when the upstream pore element is not filled by the wetting phase, while piston-like-filling occurs when the upstream element is filled by the wetting phase. The corresponding criteria [21,22] are shown in Equations (13) and (14) for snap-off and piston-like-filling, respectively.
If using the dimensionless capillary pressure as p c σ/r , the local capillary function is shown in Figure 2. Before snap-off or piston-like-filling occurs, the local capillary pressure function follows Equation (12). After the pore filling occurs, the local capillary function becomes a horizontal line (Equation (13) or (14)) until the wetting phase saturation reaches unity. A smooth technique is used in this work to make the local capillary function monotonic. The general capillary fitting correlation by Andersen et al. [44] is used. Thus, given either saturation or capillary pressure, the other situation can be reached.

Solution Scheme and Boundary Conditions
The solution scheme follows the IMPES algorithm. Firstly, Equations (8a) and (8b) are combined and discretized to yield the implicit pressure equation as where the local capillary pressure is treated in an explicit form (using the value from current time step ); i.e., , = , − ( ). After solving the pressure equations, the next step is to update the saturation explicitly using Equation (8a) a in discretized form: Setting the maximum time step size for every time step interval is necessary to ensure that the solution scheme is stable. Here, we only allow a single pore body to be filled by the wetting phase for every step. The maximum time step size estimation from the current step cannot ensure the stability of the calculation of the next step + 1 due to the severe nonlinearity of this problem. This instability issue can be prevented by rechecking the maximum water saturation changing at the + 1 step. If this is less than unity, the calculation proceeds to the next step + 1; otherwise, the water saturation at the + 1 step is recalculated using a smaller time step size. The details can be found in the work by Aziz and Settari [45]. After obtaining the saturation field, the capillary pressure is updated; then, the calculation proceeds until the preset maximum time. Moreover, the throats are treated explicitly by applying the Equations (8) and (15), which significantly reduces the nonlinearity of the conductivity of non-wetting phase.
The boundary conditions for spontaneous imbibition are set as follows: where the initially capillary pressure of every pore element is set as the maximum pressure .

Solution Scheme and Boundary Conditions
The solution scheme follows the IMPES algorithm. Firstly, Equation (8a,b) are combined and discretized to yield the implicit pressure equation as where the local capillary pressure is treated in an explicit form (using the value from current time step l); i.e., p n,l+1 = p w,l+1 − p c S l w . After solving the pressure equations, the next step is to update the saturation explicitly using Equation (8a) a in discretized form: Setting the maximum time step size for every time step interval is necessary to ensure that the solution scheme is stable. Here, we only allow a single pore body to be filled by the wetting phase for every step. The maximum time step size estimation from the current step l cannot ensure the stability of the calculation of the next step l + 1 due to the severe nonlinearity of this problem. This instability issue can be prevented by rechecking the maximum water saturation changing at the l + 1 step. If this is less than unity, the calculation proceeds to the next step l + 1; otherwise, the water saturation at the l + 1 step is recalculated using a smaller time step size. The details can be found in the work by Aziz and Settari [45]. After obtaining the saturation field, the capillary pressure is updated; then, the calculation proceeds until the preset maximum time. Moreover, the throats are treated explicitly by applying the Equations (8) and (15), which significantly reduces the nonlinearity of the conductivity of non-wetting phase.
The boundary conditions for spontaneous imbibition are set as follows: where the initially capillary pressure of every pore element is set as the maximum pressure p max c .

Dynamic Pore-Scale Network Modeling of Water Imbibition in Shale and Tight Formations
Firstly, the pore network of Barnett shale is generated using the truncated lognormal distribution of Equation (1), and the single-phase flow simulation is validated. Moghaddam and Jamiolahmady [46] presented the experimental study of the apparent gas permeability of Barnett shale samples; the corresponding pore-throat distributions were obtained from high-pressure mercury injection capillary pressure (HPMICP) experiments. We use Equation (1) to fit the experimental throat radius distribution data of HPMICP, and the mean, minimum and maximum of the throat radii are 6.2, 1 and 50 nanometers, respectively. Since the pore throat aspect ratio is non-trivial to be obtained for shale, the value is set as 1.7 following the works of Valvatne and Blunt [22]. For gas flow in shale and tight formations, the smaller the pore size, the more important the effect of a gas-solid collision; therefore, the corresponding non-Darcy effect needs to be considered. After considering the gas non-Darcy flow (the details are provided in our previous work [47]) for every pore element, the calculated apparent gas permeability and experimentally measured data show a good match, as shown in Figure 3, which demonstrates that the generated network is representative for Barnett Shale, which forms the foundation for the following two-phase flow simulation.

Dynamic Pore-Scale Network Modeling of Water Imbibition in Shale and Tight Formations
Firstly, the pore network of Barnett shale is generated using the truncated lognormal distribution of Equation (1), and the single-phase flow simulation is validated. Moghaddam and Jamiolahmady [46] presented the experimental study of the apparent gas permeability of Barnett shale samples; the corresponding pore-throat distributions were obtained from high-pressure mercury injection capillary pressure (HPMICP) experiments. We use Equation (1) to fit the experimental throat radius distribution data of HPMICP, and the mean, minimum and maximum of the throat radii are 6.2, 1 and 50 nanometers, respectively. Since the pore throat aspect ratio is non-trivial to be obtained for shale, the value is set as 1.7 following the works of Valvatne and Blunt [22]. For gas flow in shale and tight formations, the smaller the pore size, the more important the effect of a gas-solid collision; therefore, the corresponding non-Darcy effect needs to be considered. After considering the gas non-Darcy flow (the details are provided in our previous work [47]) for every pore element, the calculated apparent gas permeability and experimentally measured data show a good match, as shown in Figure  3, which demonstrates that the generated network is representative for Barnett Shale, which forms the foundation for the following two-phase flow simulation. Using the obtained pore network for Barnett shale, the initialization is conducted by setting every pore and throat element to have the local capillary pressure = 10 . Then, boundary conditions (Equation (17a) through (17d)) are applied in our dynamic pore network model, and the upper and lower boundaries are set as periodic. The water imbibition process can then be simulated directly within this proposed dynamic pore network model representing Barnett shale. Figure 4 shows the effect of the network size on the simulated spontaneous imbibition recovery with respect to time. We can see that a size of 20 × 60 is sufficient to capture the accuracy and efficiency for the following studies. The water saturation profiles of the water imbibition process are shown in Figure  5, where the average water saturation equals 0.2239. Note that the contact angle is set as 60° as the advancing contact angle is larger than the intrinsic contact angle due to the hysteresis effect [48]. All the mentioned conditions are set as the base case.  Using the obtained pore network for Barnett shale, the initialization is conducted by setting every pore and throat element to have the local capillary pressure p max c = 10 MPa. Then, boundary conditions (Equation (17a) through (17d)) are applied in our dynamic pore network model, and the upper and lower boundaries are set as periodic. The water imbibition process can then be simulated directly within this proposed dynamic pore network model representing Barnett shale. Figure 4 shows the effect of the network size on the simulated spontaneous imbibition recovery with respect to time. We can see that a size of 20 × 60 is sufficient to capture the accuracy and efficiency for the following studies. The water saturation profiles of the water imbibition process are shown in Figure 5, where the average water saturation equals 0.2239. Note that the contact angle is set as 60 • as the advancing contact angle is larger than the intrinsic contact angle due to the hysteresis effect [48]. All the mentioned conditions are set as the base case. shows the effect of the network size on the simulated spontaneous imbibition recovery with respect to time. We can see that a size of 20 × 60 is sufficient to capture the accuracy and efficiency for the following studies. The water saturation profiles of the water imbibition process are shown in Figure  5, where the average water saturation equals 0.2239. Note that the contact angle is set as 60° as the advancing contact angle is larger than the intrinsic contact angle due to the hysteresis effect [48]. All the mentioned conditions are set as the base case.  Generally, a mixed wettability is assumed for shale and tight formations. Here, we considered mixed wet conditions for Barnett shale. We assigned three types of mixed wet conditions-30% oil-wet, 50% oil-wet and 70% oil-wet-as shown in Figure 6. The corresponding spontaneous water imbibitioninduced oil recoveries are shown in Figure 7 for dimensionless time = , where is the unit conversion factor, which is equal to 0.018849 if is in minutes, in md, in decimal, in / , in cp, and L is a characteristic length in cm [49]; we follow the same units here. According to the figure, with the increasing proportion of oil-wet pores, the imbibition rate and recovery factor decrease dramatically. The results indicate that wettability dominates the water imbibition characteristics.  Generally, a mixed wettability is assumed for shale and tight formations. Here, we considered mixed wet conditions for Barnett shale. We assigned three types of mixed wet conditions-30% oil-wet, 50% oil-wet and 70% oil-wet-as shown in Figure 6. The corresponding spontaneous water imbibition-induced oil recoveries are shown in Figure 7 for dimensionless time t D = Ct k φ σ µ w 1 L 2 , where C is the unit conversion factor, which is equal to 0.018849 if t is in minutes, k in md, φ in decimal, σ in mN/m, µ w in cp, and L is a characteristic length in cm [49]; we follow the same units here. According to the figure, with the increasing proportion of oil-wet pores, the imbibition rate and recovery factor decrease dramatically. The results indicate that wettability dominates the water imbibition characteristics. Generally, a mixed wettability is assumed for shale and tight formations. Here, we considered mixed wet conditions for Barnett shale. We assigned three types of mixed wet conditions-30% oil-wet, 50% oil-wet and 70% oil-wet-as shown in Figure 6. The corresponding spontaneous water imbibitioninduced oil recoveries are shown in Figure 7 for dimensionless time = , where is the unit conversion factor, which is equal to 0.018849 if is in minutes, in md, in decimal, in / , in cp, and L is a characteristic length in cm [49]; we follow the same units here. According to the figure, with the increasing proportion of oil-wet pores, the imbibition rate and recovery factor decrease dramatically. The results indicate that wettability dominates the water imbibition characteristics.     The existence of micro-fractures is considered in this work to investigate their effect on water imbibition, as natural and hydraulic fractures are presented within shale and tight formations. Four lines of micro-fractures are assigned within the Barnett shale pore network, and the pore radius is set as 20 nanometers (c.f. Figure 8), which is 3.23 times the mean radius of the Barnett shale matrix; therefore, the micro-fracture permeability is approximately 10 times that of the matrix permeability. The corresponding imbibition recovery with respect to dimensionless time is shown in Figure 9. With the existence of micro-fractures, the imbibition rate increases significantly while the final recovery remains similar. The micro-factures contribute to the conductance of the pore network instead of the pore volume, and the increase of conductance increases the imbibition rate but not the final recovery, as the final recovery is controlled by the capillarity of the matrix pores and throats.
Energies 2020, 13, x FOR PEER REVIEW 9 of 14 The existence of micro-fractures is considered in this work to investigate their effect on water imbibition, as natural and hydraulic fractures are presented within shale and tight formations. Four lines of micro-fractures are assigned within the Barnett shale pore network, and the pore radius is set as 20 nanometers (c.f. Figure 8), which is 3.23 times the mean radius of the Barnett shale matrix; therefore, the micro-fracture permeability is approximately 10 times that of the matrix permeability. The corresponding imbibition recovery with respect to dimensionless time is shown in Figure 9. With the existence of micro-fractures, the imbibition rate increases significantly while the final recovery remains similar. The micro-factures contribute to the conductance of the pore network instead of the pore volume, and the increase of conductance increases the imbibition rate but not the final recovery, as the final recovery is controlled by the capillarity of the matrix pores and throats.  Based on the above analysis, the simulated imbibition recovery factor with respect to dimensionless time is compared with the laboratory experimental work of Barnett shale samples by Morsy and Sheng [16]. The corresponding results are shown in Figure 10. Both the effects of wettability and micro-fractures are considered to fit the experimental data. Mixed wettability for Barnett shale is found within the fitting process. Since micro-fractures are observed in the imbibition experiments, the characteristic length is reduced significantly. In the fitting results, the effective characteristic length of Barnett shale samples is approximately several millimeters, and the fraction of oil-wet pores ranges The existence of micro-fractures is considered in this work to investigate their effect on water imbibition, as natural and hydraulic fractures are presented within shale and tight formations. Four lines of micro-fractures are assigned within the Barnett shale pore network, and the pore radius is set as 20 nanometers (c.f. Figure 8), which is 3.23 times the mean radius of the Barnett shale matrix; therefore, the micro-fracture permeability is approximately 10 times that of the matrix permeability. The corresponding imbibition recovery with respect to dimensionless time is shown in Figure 9. With the existence of micro-fractures, the imbibition rate increases significantly while the final recovery remains similar. The micro-factures contribute to the conductance of the pore network instead of the pore volume, and the increase of conductance increases the imbibition rate but not the final recovery, as the final recovery is controlled by the capillarity of the matrix pores and throats.  Based on the above analysis, the simulated imbibition recovery factor with respect to dimensionless time is compared with the laboratory experimental work of Barnett shale samples by Morsy and Sheng [16]. The corresponding results are shown in Figure 10. Both the effects of wettability and micro-fractures are considered to fit the experimental data. Mixed wettability for Barnett shale is found within the fitting process. Since micro-fractures are observed in the imbibition experiments, the characteristic length is reduced significantly. In the fitting results, the effective characteristic length of Barnett shale samples is approximately several millimeters, and the fraction of oil-wet pores ranges Based on the above analysis, the simulated imbibition recovery factor with respect to dimensionless time is compared with the laboratory experimental work of Barnett shale samples by Morsy and Sheng [16]. The corresponding results are shown in Figure 10. Both the effects of wettability and micro-fractures are considered to fit the experimental data. Mixed wettability for Barnett shale is Energies 2020, 13, 4709 10 of 15 found within the fitting process. Since micro-fractures are observed in the imbibition experiments, the characteristic length is reduced significantly. In the fitting results, the effective characteristic length of Barnett shale samples is approximately several millimeters, and the fraction of oil-wet pores ranges from 70% to 50%.
Based on the above analysis, the simulated imbibition recovery factor with respect to dimensionless time is compared with the laboratory experimental work of Barnett shale samples by Morsy and Sheng [16]. The corresponding results are shown in Figure 10. Both the effects of wettability and micro-fractures are considered to fit the experimental data. Mixed wettability for Barnett shale is found within the fitting process. Since micro-fractures are observed in the imbibition experiments, the characteristic length is reduced significantly. In the fitting results, the effective characteristic length of Barnett shale samples is approximately several millimeters, and the fraction of oil-wet pores ranges from 70% to 50%.  Moreover, the effects of the initial maximum capillary pressure and the viscosity of the non-wetting phase on water imbibition are studied in this work, and the corresponding results are shown in Figures 11 and 12. According to the figures, when the initial capillary pressure increases and the non-wetting phase viscosity decreases, the corresponding imbibition rates increase. However, the final imbibition recovery factor tends to be similar, which is because the water imbibition process is a capillary-dominated flow and the recovery factor is controlled by capillarity, which is affected by pore shape, aspect ratio, contact angle, etc. Therefore, the viscous forces in the reservoir conditions increase the imbibition rates but not the final recovery factor. Moreover, the effects of the initial maximum capillary pressure and the viscosity of the nonwetting phase on water imbibition are studied in this work, and the corresponding results are shown in Figures 11 and 12. According to the figures, when the initial capillary pressure increases and the non-wetting phase viscosity decreases, the corresponding imbibition rates increase. However, the final imbibition recovery factor tends to be similar, which is because the water imbibition process is a capillary-dominated flow and the recovery factor is controlled by capillarity, which is affected by pore shape, aspect ratio, contact angle, etc. Therefore, the viscous forces in the reservoir conditions increase the imbibition rates but not the final recovery factor.    To further investigate the capillarity effect on the water imbibition of Barnett shale, the following sensitivity studies are conducted: the determination of the pore throat aspect ratio, contact angle and shape factor of the cross-area. The corresponding results for the water imbibition recovery factor with respect to dimensionless time are shown in Figures 13-15, respectively. A higher aspect ratio tends to increase the percentage of the non-wetting phase trapped by the snap-off effect, which increases the residual non-wetting phase saturation and hence the final recovery factor (c.f. Figure 13). When the contact angle and cross-area shape factor decrease, leading to an increase of capillary pressure, the imbibition rate increases at the beginning. However, the final recovery factor decreases slightly (c.f. Figures 13 and 14) because snap-off tends to occur more frequently with a lower contact angle and shape factor, as shown in Equation (13). To further investigate the capillarity effect on the water imbibition of Barnett shale, the following sensitivity studies are conducted: the determination of the pore throat aspect ratio, contact angle and shape factor of the cross-area. The corresponding results for the water imbibition recovery factor with respect to dimensionless time are shown in Figures 13-15, respectively. A higher aspect ratio tends to increase the percentage of the non-wetting phase trapped by the snap-off effect, which increases the residual non-wetting phase saturation and hence the final recovery factor (c.f. Figure 13). When the contact angle and cross-area shape factor decrease, leading to an increase of capillary pressure, the imbibition rate increases at the beginning. However, the final recovery factor decreases slightly (c.f. Figures 13 and 14) because snap-off tends to occur more frequently with a lower contact angle and shape factor, as shown in Equation (13).

Summary and Conclusions
In this paper, a dynamic pore network model is proposed, including the pore network generation, local capillary pressure function, conductance calculation and boundary conditions for water imbibition. The proposed model is applied to Barnett shale formations, and the corresponding

Summary and Conclusions
In this paper, a dynamic pore network model is proposed, including the pore network generation, local capillary pressure function, conductance calculation and boundary conditions for water imbibition. The proposed model is applied to Barnett shale formations, and the corresponding

Summary and Conclusions
In this paper, a dynamic pore network model is proposed, including the pore network generation, local capillary pressure function, conductance calculation and boundary conditions for water imbibition. The proposed model is applied to Barnett shale formations, and the corresponding laboratory imbibition experiments are matched using the dynamic established pore network model. Then, systematical sensitivity studies are conducted. According to the simulated results using our dynamic pore network model, wettability dominates the water imbibition characteristics. Besides this, the viscous effects including viscosity, initial capillary pressure and micro-fractures increase the imbibition rate, while the final recovery factor is controlled to a greater extent by the capillarity effect including the cross-area shape factor, contact angle and aspect ratio.