Introduction

Unconventional oil and gas are the main sources of global oil and gas production, compared with conventional reservoirs. However, unconventional formations have lower porosity, poorer permeability, lower pressure coefficient, lower oil saturation, complex structures, and lower natural productivity1,2,3,4. Reservoir stimulation is a significant technology for achieving unconventional hydrocarbon economic production and ensuring high-quality development5,6,7. The core theory of reservoir stimulation is to break the tight reservoir matrix through hydraulic fracturing and achieve three-dimensional stimulation of the reservoir in terms of length, width, and height8,9,10,11. The formed complex fracture network can expand the contact area between the fracture surface and the reservoir matrix. Therefore, the distance for crude oil and natural gas in the matrix to flow into the fracture from any direction can be decreased, and the required driving pressure difference can be minimized. greatly improving the overall permeability of the reservoir The ability to effectively support the hydraulic fractures and maintain the appropriate fracture conductivity directly determines the daily and cumulative production of post-fractured wells. While, the transport and distribution of proppants within the fractures are the fundamental issues.

For the study of proppant transport within fractures, researchers have conducted a series of experiments and numerical simulations using methods such as the Discrete Element Method (DEM)12,13,14, Euler-Euler Method15,16,17, and Euler–Lagrange Method18,19. The focus of research has evolved from single fractures to fracture networks, from smooth to rough fracture surfaces, and from small-scale laboratory experiments to field-scale fracture simulations. The study considers several influencing factors, including the performance of the sand-carrying fluid (such as proppant particle size20,21,22, density23,24, sand ratio25,26, viscosity, and rheology27,28), operational parameters (such as fluid displacement, wellbore structure, and perforation parameters)29,30,31, and reservoir and fracture characteristics (such as mechanical properties and fracture morphology and size)14,32,33,34,35,36. Extensive research has shown that perforation ___location can directly influence the shape of proppant packing. Increasing the displacement and viscosity of fracturing fluid, while reducing proppant density and size, can significantly extend the proppant transport distance. Proppants tend to be distributed in clusters rather than forming sand embankments within rough fractures. Additionally, the horizontal transport distance and vertical settling height of the sand bank exhibit different trends in relation to fracture tortuosity32,37. As for the study on proppant transport in field-scale fractures, Mao, S. et al.38 proposed an efficient 3D multiphase particle-in-cell (MP-PIC) method only to investigate the proppant transport among multiple field-scale fractures (fracture near the heel, middle fracture, fracture near the toe). Jiang, T. et al.39 established a planar three-dimensional (PL3D) model of multi-cluster fracturing in horizontal wells after the proppant migration simulation in the horizontal wellbore. Further, Chen, M. et al.40 developed a proppant transport simulator coupled with multi-planar 3D (multi-PL3D) fracture propagation to simulate the proppant transport among multi-cluster of fractures, considering wellbore friction, multi-fracture stress interaction, fluid leak-off, and multi-scale propagation. All the above studies only simulate the transport of proppant within multiple clusters of fractures individually and do not achieve the prediction of fracture conductivity.

In the study of the conductivity of sand-filled fractures, indoor measurements are primarily conducted using fracture conductivity test devices, while numerical simulations are performed using the Discrete Element Method41,42,43, lattice Boltzmann method44,45,46, and empirical formulas related to fracture porosity47,48. Various studies have primarily analyzed the effects of factors such as proppant particle size distribution, compaction, embedding, fragmentation, reservoir particle transport, fracture geometry and surface roughness, reservoir mechanical properties, and stress loading methods on fracture conductivity37,49,50. Research findings indicate that the compaction, breakage, and embedding of proppants under in-situ stress after fracturing are the main causes of reduced fracture conductivity. Utilizing smaller proppants can decrease fracture conductivity but also reduce stress sensitivity. Additionally, increased surface roughness of the fracture is associated with lower fracture conductivity37,46. As for the study on the large-scale fracture conductivity, Luo, Z. et al.51 conducted acid-etched fracture conductivity experiments using a novel large-scale apparatus, but the fracture was only 1 m long. Liu, F et al.52 performed experiments on acid-etched large-scale rock plate (1000 mm × 100 mm × 15 mm) to investigate the acid etched fracture length, pattern, and conductivity. All the above experiments focus on the conductivity of acid-etched fractures, with no proppant filling the fractures, and the size is not comparable to field fractures.

In summary, the previous research on fracture conductivity has relied on small-scale sand-filled core experiments, measuring just a few centimeters, or constructing digital models of particle accumulation a few millimeters in size to study the conductivity of different types of sand-filled fractures. However, actual fractures are normally over hundreds of meters long and have uneven sand distribution within them. The conductivity obtained from these small-scale sand-filled fractures cannot satisfy the requirements. Only by segmenting the fracture based on sand concentration, assigning corresponding conductivity to different regions, and then conducting a flow simulation of the whole fracture, can a more accurate large-scale fracture conductivity be obtained.

Therefore, this study first simulates the uneven distribution of proppants across multiple fracture clusters by the mixture model and segments the entire fracture based on sand concentration. Subsequently, assigning the mechanical parameters and permeability to the corresponding regions of the sand banks, fluid flow and fracture surface deformation are simulated throughout the entire fracture. Finally, the conductivity of field-scale fractures can be obtained.

Proppant transport in multi-cluster of fractures

The stimulation of unconventional reservoirs typically involves staged fracturing in the horizontal section following the drilling of a horizontal well53,54. This study’s simulation focuses on a single fracturing segment, where the physical model comprises a horizontal wellbore section and multiple clusters of field-scale hydraulic fractures. Due to the symmetry of the hydraulic fractures, proppant transport is calculated within half of each fracture cluster. The fractures are numbered based on their distance from the slurry inlet, ranging from nearest to farthest, specifically labeled as 1# fracture, 2# fracture, and 3# fracture, as illustrated in Fig. 1. The dimensions of the hydraulic fractures and the horizontal wellbore are detailed in Table 1.

Fig. 1
figure 1

Geometric configuration of multi-cluster fractures and horizontal well.

Full size image
Table 1 Geometric parameters of multi-cluster fractures and horizontal well.
Full size table

Mixture model

The multiphase flow simulation using the Eulerian-Eulerian approach is a crucial model for simulating the distribution of proppant particles in hydraulic fractures. This approach primarily includes the Volume of Fluid (VOF) model, the Eulerian model, and the mixture model. The mixture model, which is a simplified version of the Eulerian model, is particularly suitable for simulating proppant transport. Due to its relatively low computational complexity and high accuracy, this study employs the mixture model to simulate proppant distribution across multiple fracture clusters.

In the mixture model, the particle–fluid mixture is treated as a single flow continuum with macroscopic properties such as density and viscosity. The mixture consists of a dispersed phase and a continuous phase, where the continuous phase is the liquid, and the dispersed phase comprises the proppant particles. The momentum equation used in the mixture model to calculate proppant distribution within fractures is:

$$\begin{gathered} \rho \frac{\partial j}{{\partial t}} + \rho \left( {j \cdot \nabla } \right)j + \rho_{{\text{c}}} \varepsilon \left( {j_{{{\text{slip}}}} \cdot \nabla } \right)j = - \nabla p - \nabla \cdot \tau_{{{\text{Gm}}}} + \rho g{ - }\nabla \cdot \left[ {\rho_{{\text{c}}} \left( {1 + \varphi_{{\text{c}}} \varepsilon } \right)u_{{{\text{slip}}}} j_{{{\text{slip}}}}^{T} } \right] \hfill \\ - \rho_{{\text{c}}} \varepsilon \left[ {\left( {j \cdot \nabla } \right)j_{{{\text{slip}}}} + \left( {\nabla \cdot \left( {D_{{{\text{md}}}} \nabla \varphi_{{\text{d}}} } \right)} \right)j} \right] \hfill \\ \end{gathered}$$
(1)

where, ρ is the mixture density, kg/m3; j is the volume average velocity of mixture, m/s; ρc is the density of continuous phase, kg/m3; ε is the reduction of density difference between two phases in the mixture; jslip is the sliding flux, m/s; p is mixture pressure, MPa; τGm is the sum of viscosity and turbulent stress in the mixture; uslip is sliding velocity, m/s; Dmd is turbulent diffusion coefficient, m2/s; g is gravitational acceleration, m/s2; ϕc and ϕd are volume fractions of continuous and dispersed phases in the mixture respectively.

The density of mixture can be calculated by:

$$\rho = \varphi_{{\text{c}}} \rho_{{\text{c}}} + \varphi_{{\text{d}}} \rho_{{\text{d}}}$$
(2)

where, ρd is the density of dispersed phase, kg/m3.

The volume average velocity of the mixture is:

$$j = \varphi_{{\text{d}}} u_{{\text{d}}} + \varphi_{{\text{c}}} u_{{\text{c}}}$$
(3)

where, uc and ud are velocity vectors of the continuous phase and the dispersed phase, respectively.

The reduction in density difference between two phases in the mixture is:

$$\varepsilon = \frac{{\rho_{{\text{d}}} - \rho_{{\text{c}}} }}{{\rho_{{\text{c}}} }}$$
(4)

The slipping flux is:

$$j_{{{\text{slip}}}} = \varphi_{{\text{d}}} \varphi_{{\text{c}}} u_{{{\text{slip}}}}$$
(5)

where, uslip is relative velocity between the continuous and dispersed phases, m/s.

The sum of the viscosity and turbulent stress of the mixture is:

$$\tau_{{{\text{G}}m}} = \left( {\mu + \mu_{{\text{T}}} } \right)\left[ {\nabla j + \nabla j^{T} } \right] - \frac{2}{3}\left( {\mu + \mu_{{\text{T}}} } \right)\left( {\nabla \cdot j} \right){\text{I - }}\frac{{2}}{{3}}\rho k{\text{I}}$$
(6)

where, μ is the mixture viscosity, Pa·s; μT is the turbulent viscosity, Pa·s; k is available turbulent kinetic energy, m2/s2.

The turbulent diffusion coefficient of the mixture is:

$$D_{{{\text{md}}}} = \frac{{\mu_{{\text{T}}} }}{{\rho \sigma_{{\text{T}}} }}$$
(7)

where, σT is Schmidt number of turbulent particles.

The continuity equation of the mixture is:

$$\frac{\partial \rho }{{\partial t}} + \nabla \cdot \left( {\rho u} \right) = 0$$
(8)

The dispersed phase transport equation is:

$$\frac{\partial }{\partial t}\left( {\varphi_{{\text{d}}} \rho_{{\text{d}}} } \right) + \nabla \cdot \left( {\varphi_{{\text{d}}} \rho_{{\text{d}}} u_{{\text{d}}} } \right) = \nabla \cdot \left( {\rho_{{\text{d}}} D_{{{\text{md}}}} \nabla \varphi_{{\text{d}}} } \right)$$
(9)

The slipping velocity is:

$$\frac{3}{4}\frac{{C_{{\text{d}}} }}{{d_{{\text{d}}} }}\left| {u_{{{\text{slip}}}} } \right|u_{{{\text{slip}}}} = - \frac{{\left( {\rho { - }\rho_{{\text{d}}} } \right)}}{{\rho_{{\text{c}}} }}\left( { - \frac{\partial j}{{\partial t}} - \left( {j \cdot \nabla } \right)j + g} \right)$$
(10)

where, Cd is particle resistance coefficient, which can be calculated by:

$$C_{{\text{d}}} = \frac{24}{{{\text{Re}}_{{\text{p}}} }}\left( {1 + A{\text{Re}}_{{\text{p}}}^{B} } \right) + \frac{C}{{1 + D/{\text{Re}}_{{\text{p}}} }}$$
(11)

where, Rep is particle Reynolds number; Sp is proppant particle sphericity; A, B, C and D are empirical coefficients.

$$0 < S_{{\text{p}}} = \frac{{A_{{{\text{sphere}}}} }}{{A_{{{\text{particle}}}} }} \le 1$$
(12)
$$\begin{gathered} A = e^{{2.3288 - 6.4581S_{{\text{p}}} + 2.4486S_{{\text{p}}}^{2} }} \hfill \\ B = 0.0964 + 0.5565S_{{\text{p}}} \hfill \\ C = e^{{4.905 - 13.8944S_{{\text{p}}} + 18.4222S_{{\text{p}}}^{2} - 10.2599S_{{\text{p}}}^{3} }} \hfill \\ D = e^{{1.4681 + 12.2584S_{{\text{p}}} - 20.7322S_{{\text{p}}}^{2} + 15.8855_{{\text{p}}}^{3} }} \hfill \\ \end{gathered}$$
(13)

where, e is natural constant.

A similar model was employed in Yue’s simulation of proppant distribution within fracture networks55. To validate the accuracy of this model, the same parameters used in Liu’s experiment were used in the simulation56, as shown in Table 2. The comparison between experimental results and numerical results is presented in Fig. 2. It can be seen that the sand bank height at the entrance of fracture, sand bank equilibrium height and its position, as well as the sand bank height inside the fracture are basically consistent with the experimental results, which proves the accuracy of this model.

Table 2 Parameters in the experiment and validation simulation.
Full size table
Fig. 2
figure 2

Comparison between experimental results (Liu et al., 2022) and numerical results for model validation.

Full size image

However, during actual hydraulic fracturing operations, hydraulic fractures extend while sand-carrying fluid is injected into the fractures, and the placement of the proppant and the fracture propagation are dynamically correlated. This study used pre-set fixed-size fractures to simulate the transport of proppants without considering the dynamic changes of the fracture configurations, which is a limitation of this study.

The parameters used in the numerical simulation are presented in Table 3 to examine the effects of fracture height and fracture width on proppant distribution. The proppant transport simulation is conducted using the two-phase flow module in a multiphysics coupling simulation software called COMSOL Multiphysics.

Table 3 Parameters used in the simulation of proppant distribution.
Full size table

Effect of fracture height on proppant distribution

Using the aforementioned mixture model, simulations were conducted to analyze proppant distribution within hydraulic fractures under varying fracture heights and widths. The simulation results demonstrated that the proppant distribution across fracture clusters in a horizontal well section exhibits several common characteristics: (1) Within each fracture cluster, from the first fracture to the second and subsequently to the third, there is a sequential decrease in overall proppant concentration; (2) In each fracture cluster, areas with a sand concentration between 0 and 3% constitute more than half of the total fracture area, while areas with a sand concentration exceeding 12% generally account for less than 10%. (3) The sand concentration at the bottom of each fracture cluster is relatively high, and the configuration of this sand bank is regular. Conversely, the sand concentration in the middle and upper sections of the fracture is lower, and the shape of the sand bank in these regions is more irregular, as depicted in Figs. 3 and 5.

Fig. 3
figure 3

Effects of fracture height on the proppant distribution.

Full size image

Beyond these general characteristics, the specific impact of fracture height on proppant distribution in multi-cluster fractures is notable. As fracture height increases, the proppant concentration within the fracture decreases, and the vertical extent of proppant placement diminishes, particularly when the fracture height exceeds 10 m. However, with increasing fracture height, the shape of the sand bank within fractures becomes more regular, as shown in Fig. 3.

A more detailed quantitative analysis of the effect of fracture height on proppant concentration within fractures is presented in Fig. 4. The sand concentration within the fractures is categorized into five levels: 0% < ϕd ≤ 3%, 3% < ϕd ≤ 6%, 6% < ϕd ≤ 9%, 9% < ϕd ≤ 12%, and 12% < ϕd ≤ 15%. The simulation results indicate that as fracture height increases, both the overall sand filling ratio and the sand concentration in each fracture cluster decrease significantly (see Fig. 4(a)). At a fracture height of 6 m, the area with higher sand concentration in the first fracture occupies a larger portion of the fracture, with areas exceeding 12% sand concentration accounting for 32.26% (see Fig. 4(a)). Additionally, some regions in the third fracture exhibit sand concentrations exceeding 3%, accounting for 5.95% (see Fig. 4(c)). As fracture height increases, the overall sand filling ratio across the three fractures decreases sharply from 94.8% at a fracture height of 6 m to 65.5% at a fracture height of 12 m (see Fig. 4(d)). Furthermore, at a fracture height of 12 m, the sand concentration in the second and third fractures is less than 3% (see Figs. 4(b) and 4(c)).

Fig. 4
figure 4

Proppant filling ratio of different proppant concentrations in each cluster of fractures under different fracture heights.

Full size image

Based on the sand concentration contour lines depicted in Fig. 3, the sand-filled regions within each fracture cluster at varying fracture heights are segmented into distinct parts, as illustrated in Fig. 8. In Fig. 8, the different colored areas represent the displacement of fracture surfaces in different sections under the closure pressure, indirectly illustrating the segments of fractures based on proppant concentration.

Effect of fracture width on proppant distribution

The proppant concentration within multiple fracture clusters of varying widths is shown in Fig. 5. The results indicate that as fracture width increases, the sand concentration along the fracture height decreases consistently, particularly in the first fracture. When the fracture width reaches 6 mm, the sand concentration exceeds 5% in most areas along the height of the first fracture. However, as the fracture width increases to 10 mm, regions with high sand concentration become confined to the lower one-third of the fracture height. Additionally, the maximum sand concentration within the second and third fractures also decreases with increasing fracture width.

Fig. 5
figure 5

Effects of fracture width on the proppant distribution.

Full size image

Figure 6 illustrates the effect of fracture width on proppant concentration within each cluster of fractures. While variations in fracture width have a relatively minor impact on the overall sand filling ratio of the fractures, they significantly influence the distribution of sand concentration within each cluster. With increased fracture width, the sand concentration in the first and second fractures, as well as across all three clusters, decreases notably (as shown in Figs. 6(a), 6(b), and 6(d)). Specifically, the proportion of sand concentration exceeding 6% in the first and second fractures decreases from 55.48% and 24.61% at a fracture width of 6 mm to 22.41% and 1.61% at a width of 10 mm, respectively. Across varying fracture widths, the sand concentration in the third fracture remains below 3%, but the sand filling ratio decreases as fracture width increases (as depicted in Fig. 6(c)).

Fig. 6
figure 6

Proppant filling ratio of different proppant concentrations in each cluster of fractures under different fracture widths.

Full size image

The sand-filled areas within each cluster of fractures at different fracture heights are delineated in Fig. 10.

Fracture closure simulation

Following hydraulic fracturing, as the slurry pressure within the fracture decreases, the proppant particles undergo three distinct processes: suspension, compaction, and embedding into the fracture surface.

Fracture closure-contact model

For the variation of fracture width during the suspension to compaction of proppant particles, the Barree model is used for calculation, as follows:

$$w_{{\text{f}}} = 2.04C_{{\text{p}}} w_{1} - 7.51C_{{\text{p}}} m_{3} \sigma_{{\text{n}}} d_{{\text{p}}}^{{m_{4} }}$$
(14)

where, wf is the aperture of hydraulic fracturing, m; Cp is the concentration of proppants, g·cm−2; m3 and m4 are regression coefficients; σn is the effective normal stress on the fracture surface, MPa; dp is the diameter of the proppant particles, μm. The first item on the right side of the equation represents the initial height of the proppant filling in the fracture, which is not affected by stress and mainly depends on the concentration of proppants. The second item represents the decrease in aperture caused by pore structure compaction, which mainly depends on the geostress conditions and the diameter of the proppant particles.

Under effective closure stress, the fracture surface begins to deform, and the proppant particles become embedded in it. The interaction between the fracture surface and the proppant particles is modeled using elastic mechanics theory. For isotropic linear elastic materials representing both the rock matrix and the proppants, the stress–strain relationship under external forces is described as follows:

$$\rho_{{{\text{s}}i}} \frac{{\partial^{2} u_{i} }}{{\partial t^{2} }} = \nabla \cdot \sigma_{i} + F_{{{\text{v}}i}}$$
(15)

where, ρs is the solid density, kg/m3; u represents the displacement, mm; σ is the stress tensor; Fv is the external force exerted on a unit volume of rock, N. When i = 1, it represents the rock matrix, and when i = 2, it represents the proppant.

According to the principles of elasticity, the relationship between rock stress and strain is expressed as:

$$\sigma_{i} = C_{i} \cdot \varepsilon_{i}$$
(16)
$${\varvec{\varepsilon}}_{i} = \frac{1}{2}\left[ {\left( {\nabla {\varvec{u}}_{i} } \right)^{T} + \nabla {\varvec{u}}_{i} } \right]$$
(17)

where, ε is the strain tensor.

The stiffness matrix in the three-dimensional stress–strain relationship is a 6th order tensor, expressed as follows:

$${\varvec{C}}_{i} = \frac{{E_{i} }}{{\left( {1 + \upsilon_{i} } \right)\left( {1 - 2\upsilon_{i} } \right)}}\left[ {\begin{array}{*{20}c} {1 - \upsilon_{i} } & {\upsilon_{i} } & {\upsilon_{i} } & 0 & 0 & 0 \\ {\upsilon_{i} } & {1 - \upsilon_{i} } & {\upsilon_{i} } & 0 & 0 & 0 \\ {\upsilon_{i} } & {\upsilon_{i} } & {1 - \upsilon_{i} } & 0 & 0 & 0 \\ 0 & 0 & 0 & {\frac{{1 - 2\upsilon_{i} }}{2}} & 0 & 0 \\ 0 & 0 & 0 & 0 & {\frac{{1 - 2\upsilon_{i} }}{2}} & 0 \\ 0 & 0 & 0 & 0 & 0 & {\frac{{1 - 2\upsilon_{i} }}{2}} \\ \end{array} } \right]$$
(18)

where, C is the stiffness matrix; E is the Young’s modulus, GPa; υ is Poisson’s ratio.

Determining the Young’s modulus and Poisson’s ratio of sand beds with varying concentrations of proppants is essential for accurately simulating the deformation and closure of proppant-filled regions within fractures. These experiments were conducted using the RTR-1500 high-temperature pseudo-triaxial rock mechanics tester. The results are presented in Tables 4 and 5.

Table 4 Young’s modulus of sand bank of different concentrations.
Full size table
Table 5 Poisson’s ratio of sand bank of different concentrations.
Full size table

For fracture segments without proppant filling, when the effective stress applied to the fracture surface becomes sufficiently large, the fracture surfaces come into contact. A penalty function is employed to model this surface contact, which is represented as follows:

$$F_{{\text{n}}} = K_{{\text{n}}} \cdot x_{{\text{p}}}$$
(19)

where, Fn is the contact force between the contact surfaces, MPa; Kn is the contact stiffness between the contact surfaces; xp is the penetration depth between the contact surfaces, m.

Effect of fracture height on fracture closure

After segmenting regions with varying sand concentrations within each fracture cluster, a numerical simulation analysis was performed to investigate the deformation and closure of fractures under different fracture heights and widths. This analysis utilized the established elastic mechanics model, along with the Young’s modulus and Poisson’s ratio of sand beds with varying sand concentrations. The results are displayed in Figs. 7and 9. The simulation findings reveal several common characteristics in the deformation of each fracture cluster: (1) Due to variations in the mechanical properties of sand beds with different sand concentrations, significant disparities in fracture deformation are observed across different regions. High sand concentration at the bottom of the fractures provides substantial support, whereas the top of the fractures, particularly at the fracture toe, exhibits low sand concentration, leading to severe deformation or even closure under high closure pressure. (2) Fracture areas without proppant filling can remain open when the closure pressure is below 10 MPa. However, they will close as the closure pressure increases, and the irregular shapes of these areas can cause the enclosed sand-filled regions to lose functionality due to their closure. (3) As the overall sand concentration decreases with increasing distance from the sand-carrying fluid inlet in each fracture cluster, the fracture aperture correspondingly deteriorates.

Fig. 7
figure 7figure 7

Fracture deformation and closure with different fracture heights under closure pressure of 70 MPa.

Full size image
Fig. 8
figure 8

Remaining fracture width with different fracture heights.

Full size image
Fig. 9
figure 9

Fracture deformation and closure with different fracture widths under closure pressure of 70 MPa.

Full size image

The deformation and closure of fractures at closure pressures of 70 MPa, with fracture heights of 6 m, 8 m, 10 m, and 12 m, are illustrated in Fig. 7. Additionally, two vertical cross-sections, A-A’ and B-B’, located at 10 m and 60 m along each fracture, were analyzed to compare fracture surface displacement at these positions. It is observed that the fracture surface displacement at the bottom of the fractures at A-A’ and B-B’ is relatively small, allowing the fractures to maintain a larger width. In contrast, the displacement at the top of the fractures is relatively large, with the displacement curve exhibiting fluctuations, indicating alternating closed and open areas. Furthermore, as the fracture height increases, the fluctuations in the displacement curves at these two positions intensify, resulting in larger closure areas, particularly at the A-A’ cross-section.

The remaining fracture width is defined as the area-weighted average of the fracture width after deformation in each region under in-situ stress. Figure 8 illustrates the effect of fracture height on the remaining fracture width for each cluster of fractures. The results indicate that, across different fracture heights, the remaining fracture width of each cluster decreases with increasing in-situ stress. When the in-situ stress is below 20 MPa, the reduction in fracture width per unit of closure pressure is more pronounced. Once the closure pressure exceeds 20 MPa, the remaining fracture width and closure pressure exhibit a linear relationship. Among the three fracture clusters, 1# fractures and 2# fracture generally maintain a relatively larger remaining fracture width. As the fracture height increases, the difference in remaining fracture width between these two clusters gradually diminishes. At a fracture height of 12 m, the remaining fracture width of 1# fracture is generally smaller than that of 2# fracture. Focusing on fracture 1# as an example, the influence of fracture height on the remaining fracture width is depicted in Fig. 8(e). It can be concluded that as the fracture height increases, the remaining fracture width of 1# fracture continuously decreases, although the rate of reduction diminishes.

Effect of fracture width on fracture closure

Figure 9 illustrates the displacement and closure of fractures with varying widths under a closure pressure of 70 MPa. The simulation results for fractures with different widths exhibit patterns similar to those observed in fractures with varying heights. When the fracture aperture is 6 mm, the fracture tip remains largely open, although some areas near the fracture toe may close. As the fracture width increases to 8 mm, a closed region begins to form at the tip of 3# fracture. With a further increase in fracture aperture to 10 mm, closed regions appear at the tips of both 2# fractures and 3# fracture.

Figure 10 illustrates the impact of fracture aperture on the remaining fracture width for each cluster of fractures. The results demonstrate that, regardless of the fracture type, the remaining fracture width for each cluster decreases as closure pressure increases. When the closure pressure is below 20 MPa, the reduction in fracture width per unit of closure pressure is more pronounced. Beyond 20 MPa, the remaining fracture width and closure pressure exhibit a linear relationship. Among the three clusters of fractures, 1# fractures and 2# fracture generally maintain a relatively larger remaining fracture width. 3# fracture, with an initial width of 6 mm, remains largely unchanged when the closure pressure increases from 40 to 50 MPa, after which it decreases linearly. This results in a difference in the remaining fracture width between 1# fractures and 2# fracture. This behavior may be attributed to the sudden closure of a small area within the fracture. Using 1# fracture as an example, Fig. 10(d) presents the influence of fracture width on the remaining fracture aperture. It is evident that as the fracture aperture increases, the reduction in the remaining fracture width of 1# fracture becomes more significant. When the closure pressure increases from 0 to 80 MPa, initial fracture widths of 6 mm, 8 mm, and 10 mm decrease by 42.31%, 61.58%, and 63.46%, respectively. However, due to differences in initial fracture width, the fracture with an initial width of 8 mm has the smallest remaining width when the closure pressure reaches 80 MPa.

Fig. 10
figure 10

Remaining fracture width with different initial fracture widths.

Full size image

Fracture conductivity simulation

When the closure pressure is low, certain sections of the fracture that are not filled with proppants remain open, allowing fluid to continue flowing through these regions. Consequently, when simulating fluid flow in fractures, it is essential to perform a coupled analysis between the proppant-filled and unfilled sections.

Free and porous media flow model

The motion of fluid in fractures parts without proppants can be represented by the Navier–Stokes equation as follows:

$$\frac{{\partial \rho_{{\text{l}}} v_{{\text{l}}} }}{\partial t} + \rho_{{\text{l}}} \left( {v_{{\text{l}}} \cdot \nabla } \right)v_{{\text{l}}} = \nabla \cdot \left[ { - P_{{{\text{fl}}}} + \mu_{{\text{l}}} \left( {\nabla v_{{\text{l}}} + \left( {\nabla v_{{\text{l}}} } \right)^{T} } \right)} \right]$$
(20)

where, ρl is the fluid density, kg/m3; vl is the fluid velocity, m/s; Pfl is the fluid pressure in fracture, MPa; μl is the fluid viscosity, mPa·s.

The fluid flowing in fracture parts filled with proppants can be expressed as:

$$d_{{\text{f}}} \frac{{\partial \left( {\varphi_{{\text{f}}} \rho_{{\text{l}}} } \right)}}{\partial t} + \nabla \cdot \left( {d_{{\text{f}}} \rho_{{\text{l}}} v_{{\text{l}}} } \right) = d_{{\text{f}}} Q_{{\text{m}}}$$
(21)
$$v_{{\text{l}}} = - \frac{{k_{{\text{f}}} }}{{\mu_{{\text{l}}} }}\nabla P_{{{\text{fl}}}}$$
(22)

where, df is the fracture width, mm; ϕf is the porosity of fracture parts filled with proppants; kf is the permeability of fracture parts filled with different concentrations of proppants, 10−3μm2. The fracture permeability at various proppant concentrations was measured using the SL20-2D fracture conductivity tester, as shown in Fig. 11(a). Fracture conductivity was measured separately for proppant concentration ranges of 0%–3%, 3%–6%, 6%–9%, 9%–12%, and 12%–15%, as illustrated in Fig. 11(b).

Fig. 11
figure 11

Sand bank permeability equipment and experimental results.

Full size image

Therefore, the fracture conductivity can be calculated by:

$$K_{{\text{f}}} = - \frac{{Q\mu_{{\text{l}}} }}{\Delta pA}l \cdot w_{{\text{f}}}$$
(23)

where, Kf is the fracture conductivity, D·cm; Q is the fluid production, m3/d; Δp is the pressure drop within fractures, MPa; A is area of fluid flowing inside fractures, m2; l is the fracture length, m.

Effect of fracture height on fracture conductivity

Using the previously established free and porous media flow model, a numerical simulation analysis was conducted to investigate the fluid flow characteristics within fractures under varying parameters. The results indicate that fluid flow within each fracture cluster shares several common characteristics: (1) Due to variations in proppant concentration across different fracture clusters, the permeability of the sand banks in each region differs, leading to significant non-uniformity in fluid flow within the fractures. (2) Under low closure pressure, the sparse distribution of proppants at the top of the fracture results in high permeability, directing fluid flow primarily through this region. However, as closure pressure increases, the top of the fracture deforms and closes rapidly, causing a sharp decline in permeability and a gradual shift of the primary flow path towards the bottom of the fracture. (3) Compared to 1# fractures and 2# fracture, 3# fracture exhibits relatively uniform internal fluid flow due to its low proppant concentration and the presence of a single characteristic area.

Figures 12, 13, 14, and 15 illustrate the fluid velocity in each fracture cluster of varying heights at closure pressures of 10 MPa and 70 MPa, respectively. At a closure pressure of 10 MPa, the vertical cross-section A-A’ is positioned at 10 m along the fracture length, while at 70 MPa, the vertical cross-section B-B’ is similarly positioned, allowing for comparison of fluid flow rate changes under different closure pressures. The results suggest that as fracture height increases, the impact of closure pressure on fluid velocity within different fracture regions diminishes. Notably, when the fracture height is 12 m and closure pressures are 10 MPa and 70 MPa, there is no significant change in flow velocity within 2# fractures and 3# fracture, as shown in Figs. 15(c2) and 15(c3). This is primarily because these two fracture clusters contain only a single sand concentration area, limiting the effect of compaction on fluid velocity changes.

Fig. 12
figure 12

Fluid flowing velocity in each cluster of fractures with height of 6 m.

Full size image
Fig. 13
figure 13

Fluid flowing velocity in each cluster of fractures with height of 8 m.

Full size image
Fig. 14
figure 14

Fluid flowing velocity in each cluster of fractures with height of 10 m.

Full size image
Fig. 15
figure 15

Fluid flowing velocity in each cluster of fractures with height of 12 m.

Full size image

The effect of fracture height on the overall conductivity of fractures is illustrated in Fig. 16. The results indicate that, at lower closure pressures, the overall conductivity of 3# fracture is greater than that of 1# fractures and 2# fracture. However, the conductivity of 3# fracture decreases rapidly. Among the three fracture clusters, 1# fracture exhibits the lowest overall conductivity at low closure pressures but demonstrates the greatest stability. At higher closure pressures, 1# fracture maintains the highest overall conductivity among the fracture clusters.

Fig. 16
figure 16

Effects of fracture height on the fracture conductivity.

Full size image

Effect of fracture width on fracture conductivity

Figures 17, 18, and 19 illustrate the fluid velocity in each fracture cluster with varying widths at closure pressures of 10 MPa and 70 MPa, respectively. The data show that as fracture width increases, there is a more pronounced change in fluid velocity across different regions of 1# fractures and 2# fracture. Specifically, with increasing closure pressure, fluid velocity in the bottom region of these fractures increases, while the flow rate in the top region decreases. In contrast, due to the presence of a single sand concentration in 3# fracture, the longitudinal fluid velocity decreases under higher closure pressures. Additionally, as the fracture aperture increases, the influence of fracture width on fluid velocity diminishes.

Fig. 17
figure 17

Fluid flowing velocity in fractures with width of 6 mm.

Full size image
Fig. 18
figure 18

Fluid flowing velocity in fractures with width of 8 mm.

Full size image
Fig. 19
figure 19

Fluid flowing velocity in each cluster of fractures with width of 10 mm.

Full size image

Figure 20 further illustrates the effect of fracture width on the overall conductivity of fractures. The results indicate that, at lower closure pressures, the overall conductivity of 3# fracture is higher than that of 1# fractures and 2# fracture. However, the conductivity of 3# fracture decreases rapidly. When the fracture apertures are 6 mm, 8 mm, and 10 mm, and the closure pressures exceed 20 MPa, 23 MPa, and 20 MPa, respectively, the overall conductivity of 3# fracture among the three fracture clusters decreases to the lowest level.

Fig. 20
figure 20

Effects of fracture aperture on the fracture conductivity.

Full size image

Conclusions

This study divided the fractures into several regions based on the sand concentration. Utilizing the mechanical properties of sand banks in each region and applying elastic mechanics, the displacement and closure of the fracture surfaces under closure pressure were analyzed. Finally, using a free and porous media flow model, the overall fracture conductivity was calculated. The following conclusions were drawn:

(1) When the sand ratio in slurry is 15%, the maximum sand concentration in 1# fracture exceeds 12%, while the maximum sand concentration in most 3# fractures is below 3%. Fracture volume of sand concentration greater than 12% comprised less than 10% of the total fracture volume, while fracture volume of sand concentration less than 3% accounts for more than 50% of the total fracture volume.

(2) At the top of the fractures, especially at the toe, sand concentrations closed to 0% lead to significant deformation and potential closure under the closure pressures. Nonetheless, reducing the width of the fracture can maintain a larger fracture width when the closure pressures exceeds 60 MPa.

(3) Proppant distribution at the top of the fracture is sparse, leading to fluid flow primarily occurring in this region under the closure pressures of 10 MPa. As closure pressure increases to 70 MPa, the top area of the fracture deforms and closes rapidly, causing the main fluid flow channel to shift toward the fracture bottom.

(4) When the closure pressure is below 20 MPa, 3# fracture, characterized by lower sand concentration, exhibits the highest conductivity but also the greatest stress sensitivity. Whereas, 1# fracture, with higher sand concentration, experiences less conductivity decay under closure pressure, thereby maintaining greater fracture conductivity at higher closure pressures.