Impact of Proppant Embedding on Long-Term Fracture Conductivity and Shale Gas Production Decline

1 Introduction

In recent years, the petroleum industry has conducted extensive oil and gas exploration work in shale-bearing areas, and several large shale oil and gas fields have been discovered around the world. Due to the compactness of shale reservoirs, hydraulic fracturing technology must be used to modify the reservoirs. Artificial fractures generated under closing stress will gradually close. Therefore, proppants must be placed in the artificial fractures to support the fractures and maintain a certain flow capacity. However, the fractures will gradually close due to the closing stress, causing the proppants to embed into the fracture walls. This phenomenon causes changes in the flow capacity of the cracks, thereby affecting production. Li et al. [1] analyzed the laws of proppant embedding from both macro and micro perspectives through theoretical analysis and experimental research on single proppant embedding. Lu et al. [2] proposed a numerical model of the elastic-plastic deformation of proppants embedded in rough fracture surfaces and combined it with mechanical experiments. The results showed that the degree of proppant embedding varied in different areas of the rough fracture surface. Lu et al. [3]. conducted shale self-permeability experiments, established a mathematical model of hydration depth, and constructed a mathematical model of coarse crack proppant embedding that considers the entire elastic-plastic-plastic stage. The observation results show that as the closing stress increases, the degree of proppant embedding in the oil-based shale is significantly weaker than that in the distilled water-based shale. The embedding of proppants affects the flow capacity of fractures, which in turn affects overall production. Therefore, it is also very important to consider production assessment methods that take proppant embedding into account. Chen et al. [4] proposed a novel DeepONet embedded physical information neural network (DE-PINN), which demonstrated better performance than traditional data-driven machine learning models in predicting productivity and cumulative yield. Xue et al. [5] established three data-driven models to reevaluate the oil content of long-producing wells. The results showed that the radio frequency model performed best, with an accuracy rate of 95.07%, while the neural network model had the poorest predictive ability, with an accuracy rate of only 79.8%. Tian et al. [6] used Gaussian membership functions to define fuzzy levels and the entropy weight method to determine the weights of the principal components, thereby establishing a fuzzy comprehensive evaluation model for reservoir quality. Field applications showed that the average prediction accuracy rate for the cumulative oil production per kilometer of fracture length within 12 months for 10 injection wells was 91.8%. In summary, within the field of proppant embedding research, the integration of theoretical models and experimental studies has established a relatively comprehensive research framework. This spans from the microscopic behaviour of individual particles to the macroscopic design of multi-particle sand placements, from coal seams to shales, and from dry conditions to hydrated conditions. According to prior reports, factors influencing proppant embedding typically include proppant particle size, proppant modulus, surface roughness, rock modulus, and Poisson’s ratio. Notably, proppant embedding impacts long-term flow capacity; consequently, investigating its effects on sustained flow capacity and production decline is essential.

Enhancing the comprehension of the viscoelastic mechanical properties of reservoir rocks is becoming increasingly crucial due to the extensive development of unconventional oil and gas resources [7,8,9,10]. The current experimental investigation into the dynamic effects of viscoelasticity in rock materials has largely established that the power law distribution, characterized by a long tail and low attenuation, significantly outperforms the previously utilized exponential distribution for fitting creep or relaxation curves [11,12,13,14]. Fractional-order calculus, which preserves the convolutional properties of Volterra genetic integrals, is more suitable for describing time-dependent, path-dependent, and nonlocal properties [15,16,17], and the creep decay pattern obtained from the fractional-order eigenequations is precisely in the power-law form. Consequently, the utilization of fractional-order calculus for rock viscoelasticity has progressively increased in recent years [18,19,20]. Many scholars have made important contributions in the fields of creep test data correction and fractional-order numerical methods [21,22,23].

Di Paola et al. [24] posited that in relaxation trials, the variation of stressors aligns more closely with a power law distribution than with an exponential distribution. They experimentally verified the accuracy of the fractional order model and pointed out the advantages of the fractional order model over the classical component model. Ma et al. [25] proposed a rock salt creep damage model by introducing a non-steady-state modular component into the Burgers viscoelastic model. Wang et al. [26] established a new nonlinear viscoelastic plastic creep model for deep soft rock based on the Nishihara model, replacing the classical linear element with a nonlinear rheological element containing an exponential equation. This model can describe the three stages of rock nonlinearity. Similarly, Ou et al. [27] developed a novel nonlinear viscoelastic plastic soil rheological model based on the Nishihara model by incorporating a nonlinear Hookean element. This model effectively captures the rheological characteristics of both artificial and natural soft clays. Meng et al. [28] provided a fractional-order description of the time-dependent mechanical behavior of strain-softened materials. They derived leading order fitting formulas that can compute experimental curves of strain-softened materials at different degrees of strain, and the variable order problem was considered to be a superposition form of numerical integration.

A comparative analysis of the aforementioned studies on reservoir rock creep behavior indicates that some scholars have made progress in this field [29,30,31]. However, current research primarily focuses on optimizing classical models and validating theoretical frameworks through experimental data fitting, lacking a systematic integration of productivity models. Building on this foundation, this study derives empirical formulas for reservoir rock creep through experimental fitting and directly integrates them into productivity models to quantitatively analyze the dynamic impact of creep effects on gas well productivity. This approach avoids the limitations of complex theoretical derivations, clearly elucidates the correlation between the viscoelastic mechanical behavior of rocks and productivity decline, and provides actionable theoretical support for practical development.

This study establishes a conductivity attenuation model based on experimental fitting and simulates gas well productivity using an embedded discrete fracture model. Subsequently, the simulation results were compared with real-world data to validate the model’s accuracy. Finally, the effects of proppant embedment depth, rock creep, and various proppant parameters on fracture flow capacity and gas production were analyzed.

2 Experimental Fitting of Flow Conductivity Attenuation Model

In order to analyze the effects of proppant embedment and creep on fracture flow conductivity, this study investigated the decay patterns of fracture flow conductivity under different proppant filling conditions through long-term flow experiments. Empirical formulas are fitted based on experimental results to assign values to fracture unit grids in the gas well productivity model.

2.1 Origin of Experimental Rock Samples

The rock samples for the fracture flow conductivity experiments were collected from shale outcrops in the western Chongqing and Luzhou regions, which exhibit typical geological characteristics of marine shale reservoirs. Field geological surveys were conducted to select samples with intact bedding planes and well-developed natural microfracture networks.

2.2 Experimental Apparatus

The long-term flow conductivity experiments utilized the FCS-842 fracture conductivity testing system developed by TEMCO Inc. (Tulsa, OK, USA). This system is designed to determine fracturing conductivity or proppant-packed layer permeability under simulated downhole conditions in the laboratory. It evaluates the long-term fracture conductivity relationships between proppants, fracturing fluid systems, hydrocarbon types, or flow rates. Key parameters include a maximum closure pressure of 138 MPa, a maximum temperature of 177°C, and a fluid flow rate of 0–50 mL/min. The system meets API standards and satisfies the requirements for oil and gas field development (Fig. 1).

Figure 1: Fracture flow conductivity testing apparatus.

2.3 Rock Sample Preparation

Experimental rock plates were prepared perpendicular to the bedding direction to maximize fluid contact with shale bedding. The preparation steps are as follows:

Irregular shale outcrops were cut into square plates measuring 180–182 mm in length, 40–42 mm in width, and 15–17 mm in height using rock processing equipment. The end faces of the plates were polished and further refined to dimensions of 177.4 ± 0.4 mm in length, 37.6 ± 0.5 mm in width, and 12 ± 0.5 mm in thickness. The plate ends were curved to match the testing apparatus (Fig. 2), with thickness variations controlled within 0.02 mm (Fig. 3).

Figure 2: Dimensions of rock samples.

Figure 3: Processed rock plates.

2.4 Experimental Principles

The experimental apparatus operates based on Darcy’s Law: K=100QμLWdWfΔP(1) where, K is the liquid permeability of the proppant-packed layer, μm2; Q is the flow rate in the fracture, cm3/s; μ is the fluid viscosity at the test temperature, mPa·s; L is the distance between pressure ports, cm; w d is the width of the proppant-packed layer, cm; w f is the thickness of the proppant-packed layer, cm; Δ P is the pressure difference between the upstream and downstream ports, kPa.

Using an API-standard modified flow cell with L = 17.8 cm and w d = 3.8 cm, Eq. (1) is simplified to calculate the equivalent fracture permeability:

The conductivity of the proppant-packed layer is expressed as: FRCD=KWf=468QμΔP(3) where, F R C D is fracture conductivity, μm2·cm.

2.5 Experimental Procedure

Closure pressure was incrementally applied starting at 5 MPa, with steps of 15 MPa up to 50 MPa (4 pressure points total). Each pressure level was maintained for 50 ± 2 h. To simulate actual reservoir conditions, the temperature was set to 60°C, and the fluid flow rate was 2 mL/min (Table 1):

2.6 Experimental Results

The experimental results for flow conductivity are shown below (Fig. 4):

Figure 4: Experimental results of flow conductivity.

An exponential function was used to fit the long-term conductivity results under segmented closure stresses: y=Ae(−tB)+C(4) where y is the fracture conductivity, μm2·cm; t is stress application time, min; A represents the initial conductivity of the self-propped fracture, B controls the attenuation rate, and C is the long-term stable conductivity.

The fitting results are as follows (Fig. 5):

Figure 5: Fitted results of flow conductivity experiments.

As shown in Table 2, Table 3, Table 4 and Table 5 below, the larger the particle size, the higher the sand concentration and the higher the initial flow capacity. Under the same conditions, overall, expanded clay aggregate provides better support than quartz sand. The larger the particle size, the higher the sand concentration, and the lower the flow capacity decay rate.

3 Embedded Discrete Fracture Model

This study develops a gas reservoir fracturing capacity model that accounts for the decline in conductivity to examine the impact of proppant embedding and creep on fracture conductivity and shale gas production based on the following assumptions:

• (1)The fracturing fluid is distributed in the artificial main fracture and microfracture network, which is a two-phase gas-water flow, and the fracturing fluid can percolate into the matrix under the action of capillary force.
• (2)Relative independence between the gas-water phases, disregarding water evaporation and gas dissolution in water.
• (3)Capillary forces are present only in the matrix and microfracture network; the main cracks are more open and permeable, and the capillary forces are ignored.
• (4)The reservoir model is two-dimensional and fluid gravity is ignored.
• (5)Changes in temperature during fluid flow are not considered.

Then, the flux exchange between the matrix and the crack can be obtained as: QmF,k=TIkPmi,j−PFkTIk=2KmLkhfμd¯d¯=1A∫lxydA(5) where: P m i , j is the pressure at the center of the matrix grid block, MPa ; P F k is the pressure at the midpoint of the crack grid k, MPa ; K m is the permeability of the matrix grid block, D; L k is the grid step length of the crack k, m; h f is the height of the crack, m; μ is the viscosity of the fluid, mPa ⋅ s ; d ¯ is the average distance between the points inside the matrix grid block and the crack k, m; A is the area of the matrix grid block, m2; l x y is the distance from the cell to crack k in the matrix grid block, m.

The flow exchange between intersecting fractures is: Qij=TijPFj−PFiTij=TiOTjO∑k=1lTkOTkO=−2KFwfkhfkμLk,k=i,j(6) where: P F i is the pressure at the midpoint of the crack grid, MPa; P F j is the pressure at the midpoint of the crack grid, MPa; K F is the crack permeability, D; w f k is the width of the first crack, m.

Based on the basic theory of seepage mechanics, the gas seepage equation in the matrix can be expressed as: ∇βKmμgBg∇Pm+qgmVm+QmFVm=∂∂tϕmBg(7) where: β is the unit conversion factor, non-dimensional; μ g is the gas viscosity, mpa ⋅ s ; B g is the gas volume coefficient, m3/m3; q gm is the source-sink term in the matrix grid, m3/s; V m is the volume of the matrix grid block, m3; Q mF is the flurry flow rate between the matrix and the cracks, m3/s; ϕ m is the matrix porosity, dimensionless; t is the time, s.

The gas seepage equation in the fractures is expressed as: ∂∂ξβKFμgBg∂PF∂ξ+qFgVF+QmFVF=∂∂tϕFBg(8) where: ξ is the local coordinate system along the crack direction; q gF is the source-sink term in the crack grid, m3/s; V F is the crack grid block volume, m3; ϕ F is the crack porosity, dimensionless.

4 Basic Parameters and Model Validation

To validate the correctness of the model, the physical property parameters of the reservoir in the target block were selected for the simulation. The parameters for the model validation are presented in Table 6.

This model is a dual-pore medium model, where the permeability of the main fracture far exceeds that of the reservoir and microfracture network. Therefore, capillary forces in the main fracture are neglected, and the gas phase pressure is equal to the water phase pressure. Based on the continuity equation and the momentum equation, the differential equations for gas phase and water phase flow in the main fracture can be derived. The gas phase flow equation is shown below: ∂∂ξβKFKFrwμwBw∂PF∂ξ+qFw+QfFwVF=∂∂tϕFSFwBw(9) where: ξ is local coordinate system along the fracture direction; KF is permeability of the main fracture grid, mD; PF is pressure of the main fracture grid, MPa; VF volume of the main fracture grid, m3; KFrw is the relative permeability of the water phase in the main fracture grid, dimensionless; μw is viscosity of the water phase, mPa·s; Bw is volume coefficient of the water phase, m3/m3; qFw is Source-sink term for the water phase in the main fracture grid, m3/s; ϕ F is porosity of the main fracture grid, dimensionless; QfFw is water phase flow rate between the main fracture and microfracture networks, m3/s; SFw is water phase saturation in the main fracture grid, dimensionless.

The gas phase equation is: ∂∂ξβKFKFrgμgBg∂PF∂ξ+qFg+QfFgVF=∂∂tϕFSFgBg(10) where: KFrg is the relative permeability of the gas phase in the main fracture, dimensionless; μg is viscosity of the gas phase, mPa·s; Bg is volume coefficient of the gas phase, m3/m3; qFg is source-sink term for the gas phase in the main fracture, m3/s; QfFg is gas phase flow rate between the main fracture and the microfracture network grid, m3/s; SFg is gas phase saturation in the main fracture, dimensionless.

The gas phase saturation and water phase saturation in the main fracture satisfy the following relationship.

For microcrack networks, the continuity equation and momentum equation can be used to obtain the differential equations for gas and water phase flow in microcracks. The liquid phase flow equation is shown below: ∇βKfKfrwμwBw∇Pfw−QfFwδfVf+αβKmKmrwμwBwPmw−Pfw=∂∂tϕfSfwBw(12) where: Kfrw is relative permeability of the water phase in the microfracture network, dimensionless; Pfw is water phase pressure in the microfracture network, MPa; Vf is unit volume of the microfracture network, m3; α is shape factor, 1/m2; Kmrw is relative permeability of the water phase in the matrix, dimensionless; Pmw is water phase pressure in the matrix, MPa; φf is porosity of the microfracture network, dimensionless; Sfw is water phase saturation in the microfracture network, dimensionless; δf is connectivity determination coefficient, taken as 0 when there are no main fractures penetrating the computational cell; taken as 1 when main fractures are present.

The liquid phase seepage equation is as follows: ∇βKfKfrgμgBg∇Pfg−QfFgδfVf+αβKmKmrgμgBgPmg−Pfg=∂∂tϕfSfgBg(13) where: Kfrg is relative permeability of the gas phase in the microcrack network, dimensionless; Kmrg is relative permeability of the gas phase in the matrix, dimensionless; Pfg is gas phase pressure in the microcrack network, MPa; Pmg is gas phase pressure in the matrix, MPa; Sfg is gas phase saturation in the microcrack network, dimensionless.

The gas phase saturation and water phase saturation in a microcrack network satisfy the following relationship.

In microcracks, capillary forces caused by interfacial tension play an important role in fluid distribution and cannot be ignored. The gas phase and water phase pressures in the microcrack network satisfy the following equations: PfcSfw=Pfg−Pfw(15) where: Pfc is the capillary force in microcracks, MPa.

The continuous equation and the equation of motion can be used to obtain the differential equations for gas phase and water phase seepage in the matrix. The liquid phase seepage equation is shown below: ∇βKmKmrwμwBw∇Pmw+αβKmKmrwμwBwPmw−Pfwδm=∂∂tϕmSmwBw(16) where: Smw is water phase saturation in the matrix, dimensionless; δm is connectivity judgment coefficient, taken as 0 when there is no microcrack grid in the calculation unit, and 1 otherwise.

The gas phase seepage equation in the matrix can be expressed as: ∇βKmKmrgμgBg∇Pmg+αβKmKmrgμgBgPmg−Pfgδm=∂∂tϕmSmgBg+ρsVLPmgPmg+PL(17) where: Smg is gas phase saturation in the matrix, dimensionless; ρs is matrix density, kg/m3; VL is Langmuir volume, m3/kg; PL is Langmuir pressure, MPa.

The gas phase saturation and water phase saturation in the matrix system satisfy the following relationship:

The gas phase pressure, liquid phase pressure, and capillary force in the matrix system satisfy the following equation: PmcSmw=Pmg−Pmw(19) where: Pmc is the capillary force in the matrix, MPa.

The cumulative gas production of this study’s model and comparison curve of the cumulative gas production of real gas wells were obtained by setting the flow pressure at the bottom of the well to 10 MPa (Fig. 6). The figure illustrates that the two cumulative gas production curves are closely aligned, yet the disparity between them will progressively widen over time, consistent with the established laws, thereby validating the accuracy of the model presented in this paper.

Figure 6: Cumulative gas production comparison curve.

5 Simulation Results and Discussion

5.1 Simulation Results under Different Flow Diversion Capacity Decay Modes

The variable conductivity model, which accounts for proppant embedding, and the conductivity attenuation model, which considers rock creep, are utilized in the established capacity model to simulate and analyze the capacity under various conductivity attenuation scenarios. The simulation results show that when cracks occur, they affect production. That is, the geometric shape of the cracks is related to production, which is consistent with the conclusions of Marembo Micheal et al. [32].

The aforementioned figure illustrates that the disparity between the capacity simulation results and those obtained without accounting for conductivity attenuation was relatively small when only proppant embedding was considered. However, the results of the capacity simulations produced large differences from those of the other two sets of simulations after considering rock creep. A subsequent study indicated that the influence of rock creep on gas well production was more significant than that of proppant embedding, thereby highlighting that rock creep predominantly affects the fracture conductivity. As shown in Fig. 7d, the simulation results of the production capacity model that took into account proppant embedding and rock creep showed a trend of first increasing and then decreasing in daily production compared to the other two sets of simulation results.

Figure 7: Comparison of pressure cloud and capacity simulation under different conductivity attenuation models: (a) Considering proppant embedment and creep; (b) Considering proppant embedment; (c) The attenuation of conductivity is not considered; (d) Comparison of capacity simulations under different modes of conductivity attenuation.

5.2 Study on the Influence of Fracture Flow Conductivity

To determine the effects of proppant parameters on fracture conductivity and gas well performance, simulations were conducted for reservoirs with proppant modulus values in the range of 5000 < E2 < 40,000 (Fig. 8).

Figure 8: Comparison of capacity simulation with different proppant sizes: (a) Considering proppant embedment and creep; (b) Considering proppant embedment; (c) The attenuation of conductivity is not considered; (d) Comparison of capacity simulations under different modes of conductivity attenuation.

Impact of Proppant Concentration on Productivity. Reducing proppant concentration from 5 kg/m2 to 2.5 kg/m2 resulted in an average productivity decline of 9.20%. The maximum decline (11.11%) occurred with quartz sand (70/140 mesh), while the minimum decline (7.48%) occurred with quartz sand (40/70 mesh). The ceramic proppant group showed an average decline of 9.10%, and the quartz sand group 9.29%, with minimal difference between the two.

Impact of Particle Size on Productivity. Reducing particle size from 40/70 mesh to 70/140 mesh led to an average productivity decline of 3.52%. The maximum decline (5.13%) occurred with ceramic proppants at 2.5 kg/m2, while the minimum decline (0.57%) occurred with quartz sand at 5 kg/m2. The ceramic proppant group exhibited an average decline of 4.52%, and the quartz sand group 2.52%.

Impact of Proppant Type on Productivity. Switching from ceramic proppants to quartz sand caused an average productivity decline of 2.59%. The maximum decline (4.13%) occurred at 5 kg/m2 with 40/70 mesh quartz sand, while the minimum decline (0.80%) occurred at 5 kg/m2 with 70/140 mesh quartz sand.

Conclusions. The order of impact on productivity is proppant concentration (average decline 9.20%) > particle size (3.52%) > proppant type (2.59%). Proppant concentration is the key factor. Particle size has a greater impact on ceramic proppants than on quartz sand. Under the same conditions, increasing particle size enhances productivity more effectively for ceramic proppants. For optimal productivity, ceramic proppants should be used with high concentrations and large particle sizes. Considering cost, quartz sand is recommended with high concentrations and small particle sizes to ensure long-term reservoir stimulation.

6 Summary

This study established a fracture conductivity decay model based on experimental fitting. By integrating a nested discrete fracture model, it investigated production capacity simulations under different fracture conductivity decay patterns while examining the influence of proppant parameters on fracture conductivity. The following conclusions were drawn:

• 1.The initial crack drainage capacity of sand with a concentration of 5 kg/m2 is higher than that of sand with a concentration of 2.5 kg/m2. However, as the closing stress increases, the crack drainage capacity decreases rapidly and stabilizes at a stress of 50 MPa. Expanded clay aggregate has a higher initial crack drainage capacity than quartz sand.
• 2.The influence of rock creep on fracture conductivity is more significant than that of proppant embedment on fracture conductivity.
• 3.For damage-controllable reservoirs with a viscous modulus ranging from 5000 to 40,000 MPa, the order of influence of a single-proppant system on productivity is: proppant concentration (average attenuation of 9.20%) > particle size (3.52%) > proppant type (2.59%).