As an unconventional gas resource with massive reserves, shale gas has received wide attention since it entered the commercial development stage (Zou et al., 2015a; Zou et al., 2016). According to the US Energy Information Administration (EIA) forecast, shale gas reserves in China account for about 20% of the world, ranking first globally (Ma et al., 2017), and thus shale gas has a broad development prospect in China (Dong et al., 2015; Hao et al., 2017; Sun et al., 2020). In this decade, shale gas production in China has grown steadily every year, but it mainly comes from middle and shallow shale gas reservoirs with a depth of smaller than 3,500 m. Deep shale gas reservoirs with a depth of over 3,500 m have received particular attention and have been regarded as the new growth point for China’s shale gas production in the future (Ma et al., 2021). Deep shale gas reservoirs are difficult to exploit because of the deep shale’s low porosity and ultra-low permeability. Horizontal well drilling technology (Nie et al., 2013) and hydraulic fracturing technology (Rongved et al., 2019) are the key technologies for the commercial development of shale gas reservoirs. Fracturing technology has been continuously improved and developed in recent years (Wang et al., 2015; Li et al., 2018; Cheng et al., 2019; Wu et al., 2022; Yang et al., 2022). After hydraulic fracturing, the flowback stratagem of fracturing fluid is one of the critical factors for shale gas production (Wang et al., 2020a; Qu et al., 2021).

The flowback rate of fracturing fluid is generally low and varies significantly for different shale gas reservoirs (Lutz et al., 2013; Vengosh et al., 2014; Gao et al., 2017). Compared with conventional gas reservoirs, deep shale gas reservoirs under high-temperature and high-temperature conditions have rich micro-nano pores with various pore structures (Guo et al., 2014; Chen et al., 2019), which affects the gas-liquid flow law and the flowback characteristics of gas wells. Because the flowback characteristics of gas wells are the macroscopic manifestation of the microcosmic mechanism of the imbibition and flowback, it is crucial to clarify the microcosmic mechanism of imbibition and flowback of fracturing fluid under reservoir conditions for the efficient development of deep shale gas reservoirs (Zheng et al., 2021a; Qu et al., 2021; Shao et al., 2022; Xu et al., 2022).

Lattice Boltzmann method (LBM) is considered as a powerful tool for simulating many complex flows such as multiphase flows (Wang et al., 2023), multicomponent flows (Ren et al., 2019), suspended particle flows (Jiang et al., 2022), blood flows (Cherkaoui et al., 2022), liquid slip flows (Ren et al., 2023), turbulent flows (Kareem and Asker, 2022), nanofluid heat flows (Khan et al., 2022), and so on. Recently LBM has been widely used to simulate pore-scale flows of multiphase fluid in porous media because it is easy to deal with complex geometric boundaries and the interactions between multiphase fluid and fluid and between fluid and solid (Ansumali and Karlin, 2005; Chikatamarla and Karlin, 2009; Zhang et al., 2009). At present, LBM models used for multiphase flow simulations mainly include pseudo-potential model (Shan and Chen, 1993), color gradient model (Gunstensen et al., 1991), free energy model (Swift et al., 1995) and phase field model (Shao et al., 2013). Combined with digital core technology (Zhu et al., 2012), LBM is widely used to study the seepage mechanism of multiphase fluid in porous media (Cantisano et al., 2013). Furthermore, compared with LBM, the molecular dynamics method is also a powerful tool to study the microcosmic mechanisms (Huang et al., 2020; Huang et al., 2021a; Huang et al., 2021b). But molecular dynamics is subjected to expensive computational costs, and so it is usually applied to multiphase flows in a simple nanoscale pore structure rather than in porous media with complex pore structure. In recent years, with the wide application of hydraulic fracturing technology in the development of unconventional oil and gas reservoirs (Zou et al., 2015b; Zhang et al., 2018; Li et al., 2022), some scholars have begun to use LBM to study the microscopic imbibition and flowback mechanism of fracturing fluid. In 2015, Zhang et al., 2015 used LBM to simulate the process of gas flooding in saturated water cores based on digital shale cores with a porosity of 16.5% and confirmed that most of the water would remain in the pores and could not be discharged. In 2018, Zheng et al., 2018 used LBM to simulate the influence of wettability heterogeneity distribution on the propagation and distribution of fracturing fluid based on a digital shale core with a porosity of 26.98%. In 2020, Wang et al., 2020b used LBM to simulate the process of water imbibition and oil displacement under the action of capillary force based on the tight sandstone with porosity ranging from 10.3% to 18.6% and studied the influence of pore structure and wettability on spontaneous imbibition. In 2021, Zheng et al., 2021b used LBM to simulate the imbibition and displacement of heavy oil under reservoir conditions based on the digital core of tight sandstone with a porosity of 9.47% and analyzed the distribution law of the oil-water phase. Zhang et al., 2021 developed a pseudopotential-based lattice Boltzmann (LB) method to simulate gas/water two-phase flow at pore scale and simulate gas/water two-phase flow in dual-wettability nanoporous media.

Although some scholars use LBM to simulate the imbibition (or flowback) process of fracturing fluid, they only simulate a single imbibition (or flowback) process and do not study the coupled imbibition-flowback process. The imbibition process affects the flowback behavior because the gas-liquid saturation distribution at the beginning of flowback is determined by the gas-liquid saturation distribution after imbibition. It is necessary to simulate the coupled imbibition-flowback process to thoroughly investigate the flowback mechanism of fracturing fluid. In addition, the porosity of the digital core used to simulate the imbibition (or flowback) process of fracturing fluid is relatively large, resulting in the connectivity of the pore structure being good, which is inconsistent with the actual pore structure of deep shale.

In this work, we reconstruct a digital core of deep shale, use LBM to simulate the coupled imbibition-flowback process of fracturing fluid, study the microcosmic mechanism of the imbibition and flowback, and analyze the influence law of critical parameters on the flowback of fracturing fluid under reservoir conditions.

The three-dimensional nineteen-discrete-velocity (D3Q19) pseudo-potential model with multiple relaxation times (MRT) (Shan and Chen, 1993; Zheng et al., 2018) is carried out to simulate the imbibition and flowback processes of fracturing fluid. The MRT model has much better numerical stability than the Bhatnagar-Gross-Krook (BGK) model in the simulation of the flows involving complex boundaries such as porous media flow. In the pseudo-potential model, k particle distribution functions are used to represent k component fluids. The evolution equation of the distribution function of the k component is: f α k x + e α Δ t , t + Δ t − f α k x , t = Ω c o l l k + Ω f o r c e s k (1) where f α k x , t is the equilibrium distribution function of the k component in e α direction, x is the spatial coordinate, e α is the lattice speed in α direction, Δ t is the time step, t is the time, Ω c o l l k is the collision operator of the k component, and Ω f o r c e s k is the external force operator of the k component.

The collision operator Ω c o l l k of the k component can be expressed as: Ω c o l l k = M − 1 Λ k M f α e q , k x , t − f α k x , t (2) where M is the conversion matrix; Λ k is a diagonal matrix about relaxation times, expressed as Λ k = s c k , s e k , s ε k , s c k , s q k , s c k , s q k , s c k , s q k , s v k , s π k , s v k , s v k , s v k , s v k , s v k , s m k , s m k , s m k T , f α e q , k x , t is the equilibrium distribution function, which is expressed as: f α e q , k x , t = w α ρ k x , t 1 + 3 e α ⋅ u e q c 2 + 9 e α ⋅ u e q 2 2 c 2 − 3 u e q 2 2 c 2 (3) where w α = 1 / 3 α = 0 , w α = 1 / 18 α = 1 − 6 , w α = 1 / 36 α = 7 − 18 , ρ k is the density of the k component, c = Δ x / Δ t is the lattice sound velocity with c = 1 in this study, and u e q represents the equilibrium velocity.

The conversion matrix M is expressed as M = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 − 30 − 11 − 11 − 11 − 11 − 11 − 11 8 8 8 8 8 8 8 8 8 8 8 8 12 − 4 − 4 − 4 − 4 − 4 − 4 1 1 1 1 1 1 1 1 1 1 1 1 0 1 − 1 0 0 0 0 1 − 1 1 − 1 1 − 1 1 − 1 0 0 0 0 0 − 4 4 0 0 0 0 1 − 1 1 − 1 1 − 1 1 − 1 0 0 0 0 0 0 0 1 − 1 0 0 1 1 − 1 − 1 0 0 0 0 1 − 1 1 − 1 0 0 0 − 4 4 0 0 1 1 − 1 − 1 0 0 0 0 1 − 1 1 − 1 0 0 0 0 0 1 − 1 0 0 0 0 1 1 − 1 − 1 1 1 − 1 − 1 0 0 0 0 0 − 4 4 0 0 0 0 1 1 − 1 − 1 1 1 − 1 − 1 0 2 2 − 1 − 1 − 1 − 1 1 1 1 1 1 1 1 1 − 2 − 2 − 2 − 2 0 − 4 − 4 2 2 2 2 1 1 1 1 1 1 1 1 − 2 − 2 − 2 − 2 0 0 0 1 1 − 1 1 1 1 1 1 − 1 − 1 − 1 − 1 0 0 0 0 0 0 0 − 2 − 2 2 2 1 1 1 1 − 1 − 1 − 1 − 1 0 0 0 0 0 0 0 0 0 0 0 1 − 1 − 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 − 1 − 1 1 0 0 0 0 0 0 0 0 0 0 0 1 1 − 1 − 1 0 0 0 0 0 0 0 0 0 0 0 1 − 1 1 − 1 − 1 1 − 1 1 0 0 0 0 0 0 0 0 0 0 0 − 1 − 1 1 1 0 0 0 0 1 − 1 1 − 1 0 0 0 0 0 0 0 0 0 0 0 1 1 − 1 − 1 − 1 − 1 1 1 (4)

The external force operator Ω f o r c e s k of the k component can be expressed as: Ω f o r c e s k = Δ t 2 f α F , k x + e α Δ t , t + Δ t − f α F , k x , t (5) where f α F , k = 3 F k ⋅ e α − u e q ρ k c 2 f i e q , k (6)

Let f ¯ α k = f α k − Δ t 2 f α F , k , then Eq. 1 can be rewritten as: f ¯ α k x + e α Δ t , t + Δ t − f ¯ α k x , t = M − 1 Λ k M [ f α e q , k x , t − f ¯ α k x , t − Δ t 2 f α F , k ] + Δ t f α F , k (7)

Fluid density ρ k , fluid velocity u , and fluid viscosity μ k are defined as follows: ρ k = ∑ α f ¯ α k (8) ρ k u k = ∑ α f ¯ α k e α + Δ t 2 F k (9) u = ∑ k ρ k u k ∑ k ρ k (10) μ k = ρ k c 2 3 1 s v k − 0.5 Δ t (11)

In the pseudo-potential model, the total force F k of the k component mainly includes the gas-liquid force F c , k and the fluid-solid force F a d s , k , which are defined as: F k = F c , k + F a d s , k (12) F c , k x , t = − G c ρ k x ∑ i w i ρ k ¯ x + e α Δ t , t e α (13) F a d s , k x , t = − G a d s , k ρ k x , t ∑ i w i s x + e α Δ t e α (14) where G c is the gas-liquid interaction parameter, which is related to the interfacial tension; G a d s , k is the fluid-solid interaction parameter, which is related to wettability; s x + e α Δ t is a marking function, if x + e α Δ t is a solid point, then s x + e α Δ t = 1 , otherwise s x + e α Δ t = 0 .

Compared with the other reconstruction methods, Markov Chain-Monte Carlo (MCMC) method (Wu et al., 2006) has higher computational efficiency and can obtain the pore structure of porous media with good connectivity. Therefore, the MCMC method is used to reconstruct the shale digital core in this study. Deep shale cores are from Luzhou block in southern Sichuan. The scanning electron microscope (SEM) scanning image of the deep shale core is shown in Figure 1A. The SEM scanning image is denoised by the median filter algorithm (Dong and Blunt, 2009). The images before and after denoising are shown in Figure 1A and Figure 1B, respectively. It is obvious that the median filter algorithm can obtain a clearer grayscale image. Then the grayscale image is converted to a binarization image, shown in Figure 2, by the Otsu algorithm (Otsu, 1979). The conditional probabilities in three orthogonal directions are obtained by the traversing algorithm. Based on the conditional probabilities, a 3D digital core is reconstructed by the MCMC method. Figure 3 shows the reconstructed 3D digital cores of deep shale with and without isolated pores. The size of the digital core is 130 × 130 × 100 with the resolution of 34 nm per pixel. The effective porosity of the digital core is 5.09%. Figure 4 shows the comparison between autocorrelation functions of the SEM scanning image and digital core. It is seen that the autocorrelation function of the digital core is consistent with that of the SEM scanning image, which confirms the representativeness of the reconstructed digital core (Wu et al., 2006).

According to the reservoir condition of Luzhou block in southern Sichuan, the basic parameters of LBM simulation are set as follows: the reservoir pressure p = 80 MPa , the reservoir temperature T = 130 ° C , and the contact angle θ = 10 ° . Based on National Institute of Standards and Technology (NIST) chemical database, some physical parameters of gas and water under p = 80 MPa ; T = 130 ° C are obtained as follows: the methane density ρ g = 259 kg / m 3 , the water density ρ w = 972 kg / m 3 , the methane viscosity μ g = 33.9 μ Pa ⋅ s , the water viscosity μ w = 233 μ Pa ⋅ s , and the gas-water interfacial tension σ = 40 mN / m . In LBM simulation, 15 layers of pore lattices are added to the inlet and outlet of the 3D digital core to eliminate the end effect. The final size of the digital core used for simulation is 130 × 130 × 130 lattices. The pressure boundary condition (i.e., the non-equilibrium bounce-back boundary condition) is applied to the inlet and outlet, and the bounce-back condition is applied to the solid walls.

In the simulation of the imbibition process, the pore space of the digital core is fully saturated with shale gas (i.e., methane). The fracturing fluid (i.e., water) derived by the capillary force infiltrates into the digital core and replaces the shale gas. The driving force in the imbibition process is the capillary force, which is related to the size of the pore throat. The stopping time of the simulation is 6 million time steps when the stable distribution of the two-phase fluid in the digital core can be obtained. Figure 5 shows the gas-liquid distribution in shale at different time steps and the water saturation versus time steps in the imbibition process. It is seen that the fracturing fluid advances almost uniformly throughout the pore space during the imbibition process. More liquid phase is distributed near the inlet than near the outlet. Water saturation increases with time, but the imbibition rate decreases until it finally becomes almost zero. This observation is because capillary force is the driving force in the imbibition process, and the greater the distance of fracturing fluid penetrating the core, the larger the viscous resistance of fluid migration to be overcome. Therefore, the imbibition rate gradually decreases with time, and eventually, the imbibition becomes equilibrium.

Figure 6 shows the slice diagrams of gas-liquid distribution at the equilibrium state ( t = 6.0 × 10 6 ) in the imbibition process. It is seen that the fracturing fluid is easy to adsorb on the pore wall, and the shale gas is located in the middle of the pore space. Figure 7 shows the gas-liquid distribution of slice 3 in Figure 6 at different time steps. It is observed that the fracturing fluid first invaded the pore along the pore wall and then gradually occupied the central position of the pore. The larger capillary force in tiny pores results in the fracturing fluid preferentially filling smaller pores and subsequently invading larger pores.

The flowback of fracturing fluid is simulated by LBM. the gas-liquid distribution in the digital core at the initial moment of the flowback process is the same as that at the end of the imbibition process, shown in Figure 6. The shale gas is injected continuously at the right outlet, and thus the shale gas displaces the fracturing liquid from right to left. The flowback pressure difference is set as 3.0 MPa. The boundary conditions and the basic parameters are the same as those applied in the case of the imbibition process. The stopping time of the simulation is 4 million time steps when the gas-liquid distribution in the digital core almost tends to be stable.

Figure 8 shows the gas-liquid distribution in shale at different time steps and the water saturation, gas saturation, and flowback rate versus time steps in the flowback process. It is seen that the viscous fingering is clearly observed during the flowback process, where shale gas flows through large pores to form a flow channel, and the fracturing fluid stays in tiny pores, resulting in a low flowback rate of the fracturing fluid. The flowback rate increases gradually with the flowback time and eventually tends to be almost constant.

Figure 9 shows slice diagrams of gas-liquid distribution at the equilibrium state ( t = 4.0 × 10 6 ) in the flowback process. The capillary force is resistance in the flowback process. Because of a smaller capillary force in a large pore, the gas phase preferentially breaks through in the middle of the macro-pores, and the fracturing fluid in small pores and blind ends is difficult to flowback.

Figure 10 shows the gas-liquid distribution of slice 2 in Figure 9 at different time steps. It is seen that with the increase of time, the proportion of the area occupied by the gas phase in the slice gradually increases, and the increased gas phase was mainly concentrated in the large pores of the slice.

In this section, based on the LBM simulations of the imbibition and flowback processes, the influence of the wettability, flowback pressure difference, pore structure, and fracturing fluid viscosity, on the imbibition and flowback behaviors are studied in detail. The simulation conditions and basic parameters are the same as those applied in imbibition and flowback cases in section 4 and 5, respectively.

In this paper, based on the SEM images of deep shale cores in Luzhou block, a 3D digital core was reconstructed by using the MCMC method. The LBM simulation of the imbibition and flowback processes in the reconstructed core was carried out, and the micro-mechanism of the imbibition and flowback of fracturing fluid was deeply studied. At the same time, the influence law of related parameters on the imbibition and flowback of the fracturing fluid was analyzed. The main results and understandings are as follows.

(1) In the imbibition process of the fracturing fluid, the distribution of the liquid phase near the inlet is more than that at the outlet, the imbibition rate gradually weakens with time, and finally the imbibition tends to be balanced. The fracturing fluid shows a uniform invasion, and the fracturing fluid is easily adsorbed on the pore wall, and the gas phase is located in the central position of the pore.