Hydraulic Fracture Conductivity Loss Mechanisms for Unconsolidated Sands Considering Fine Migrations and Proppant Embedments

1  Introduction

Significant oil and gas reserves have been identified in unconsolidated sandstone reservoirs, and as such, the efficient development of these reservoirs has attracted significant attention in recent years. Hydraulic fracturing is widely recognized as a critical technique for production stimulation in unconsolidated sandstone reservoirs, as high-conductivity hydraulic fractures provide sustainable flow pathways for hydrocarbons while also helping control sand production. Therefore, achieving sustained hydraulic fracture conductivity is crucial in the design of hydraulic fracturing treatments for the effective stimulation of unconsolidated sandstone reservoirs [1−5].

Typically, hydraulic fracture conductivity undergoes a decline over production time due to multiple factors, including in-situ formation geomechanical characteristics, closure stress, proppant placement quality, and even fracturing fluid-proppant incompatibility. The key mechanisms underlying fracture conductivity loss have been extensively studied in existing literature [6−10]. Beyond the mechanical properties of proppants themselves, the interaction between fracturing fluids and reservoir rocks also affects proppant embedment by altering rock surface properties—this phenomenon has been verified in shale reservoirs but remains unclarified in low-strength unconsolidated sands. Recent experimental studies highlight the importance of the apparent Young’s modulus of proppant packs in alleviating fracture aperture reduction under in-situ closure stress, demonstrating that smaller proppant sizes and higher concentrations enhance the stability of fracture packs and reduce deformation. Additionally, a few laboratory experiments have investigated the effect of shale-fluid interactions on proppant embedment, revealing the influence of fracturing fluids on fracture conductivity [11]. To provide further insights into fracture conductivity loss, numerical simulations have been proposed [12,13]. Some authors have used computational fluid dynamics to quantify fracture conductivity loss and indicated that proppant breakage, embedding, and particle migration can affect fracture conductivity [14]. Yan et al. [15] developed a simplified conceptual model to estimate fracture conductivity in key channel fracturing scenarios, demonstrating that the Brinkman equation effectively describes fluid flow within open channels bounded by proppant pillars in fractured formations. Instrumented indentation testing has emerged as a robust method for characterizing in-situ rock mechanical properties and predicting proppant embedment, particularly in soft reservoir formations such as the Montney siltstone—where conventional conductivity tests pose implementation challenges and often yield inconsistent results [16]. Guerra et al. [17] examined conductivity impairment caused by fluid-rock interactions in shale reservoirs and identified key mechanisms, including clay swelling, surface softening, excessive proppant embedment, and fines migration. Tang et al. [18] evaluated the effect of shear deformation on proppant embedment in tight gas sandstone reservoirs and found that shear stress and relative rock movement influence proppant behavior. Xu et al. [19] investigated the impact of fracture tortuosity and proppant sphericity on fracture conductivity. Furthermore, advanced analytical models developed using the Drucker-Prager yield criterion have been established to characterize the shale’s elastoplastic behavior under high in-situ stress, offering more accurate predictions of proppant embedment depth than the classical Hertzian contact theory [20]. However, most existing research has focused on proppant embedment in hard sedimentary reservoir formations (e.g., sandstones, carbonates, and even shales), and Most existing studies analyze the effects of fines migration or proppant embedment separately in unconsolidated sands, but overlook their synergistic effect—which is the core cause of rapid fracture conductivity decay in unconsolidated sands, leaving a clear research gap. Considering the inherently low mechanical strength of such unconsolidated reservoir formations, proppant embedment is a major but not the only concern, especially under low to medium in-situ closure stress conditions. In comparison, mobilizable formation fines can dislodge from the formation matrix and migrate within the fracture network, further impairing fracture conductivity [21,22]. These issues are particularly critical in weakly cemented unconsolidated sandstone reservoirs and hydrate-bearing sediments, where hydrate dissociation processes significantly diminish sediment strength and facilitate severe fines migration [23]. In unconsolidated sandstone reservoirs, shear failure often dominates, with tensile failure also occurring and posing significant sand production risks. The migration and deposition of formation fines within proppant packs can significantly reduce fracture conductivity, a phenomenon especially observed in coalbed methane reservoirs where fines aggregation and pore-throat plugging are prevalent.

To investigate the migration mechanisms of fines during hydraulic fracturing for unconsolidated sandstone formations, this study has conducted a series of experimental studies and numerical simulations. Numerous researchers have developed models to simulate fines migration within hydraulic fracture networks. Zou et al. [24] conducted experimental studies on coal fines, demonstrating that these fines tend to aggregate and become trapped within proppant packs. Mitchell and Leonardi [25] used LBM-DEM simulations to elucidate particle detachment mechanisms, while Ahamed et al. [26] highlighted the roles of proppant crushing and diagenesis in fines migration. Recent advances in coupled CFD-DEM and DEM-DFM methods for particle-scale modeling have facilitated high-fidelity simulations of particle migration and clogging behaviors in granular proppant media [27]. In particular, in hydrate-bearing sediments, the significant reduction in pore water salinity caused by hydrate dissociation can further destabilize clay minerals (e.g., illite), thereby exacerbating fines migration and pore-throat clogging [28]. Recent studies on shale reservoir formations have highlighted the critical role of time-dependent creep deformation in fracture closure and proppant embedment. For example, Benge et al. [29] demonstrated that clay-rich “ductile” zones in shale formations are highly prone to creep deformation, resulting in substantial long-term fracture conductivity loss. Additionally, fractional calculus-based creep models for shale have been developed to more effectively characterize the nonlinear viscoelastoplastic behavior of shale under sustained in-situ stress, thus improving predictions of long-term fracture conductivity decay [30]. These models highlight the importance of the formation’s viscous modulus and proppant properties in controlling long-term fracture performance—key insights that may also apply to unconsolidated sands, despite their distinct mechanical properties. Recent coupled CFD-DEM simulations focused on fines migration have demonstrated that cohesive fines tend to form agglomerated matrices under reservoir fluid flow, significantly increasing the clogging potential in rock media [31]. This process typically consists of four stages: fines mobilization, invasion into the proppant pack, fines retention, and subsequent flowback. The same simulations also revealed that fines migration and retention within the proppant pack significantly affect long-term hydraulic fracture conductivity. Recent experimental studies on weakly cemented unconsolidated reservoirs have revealed that both proppant embedment and fines migration markedly reduce fracture conductivity, with the latter being particularly severe under rapid depressurization conditions. However, a comprehensive model that integrates both proppant embedment and fines migration mechanisms for unconsolidated sand reservoirs is still lacking.

To fully address the key mechanisms underlying hydraulic fracture conductivity loss in unconsolidated sandstone reservoirs, this study aims to fill the identified gap by developing an integrated experimental and theoretical framework for assessing fracture conductivity decline—one that explicitly integrates both mechanisms. In addition, the mathematical model proposed in this study has been validated against the experimental data collected in this study, and the effects of fines concentration and reservoir fluid flow velocity on fracture conductivity have been quantitatively analyzed. To systematically reveal the synergistic effect of these two mechanisms, this study employs a combined experimental and numerical strategy: the experimental research aims to macroscopically quantify the decline law of fracture conductivity under varying closure pressures, fines concentrations, and proppant parameters; whereas the subsequent CFD-DEM numerical simulation focuses on elucidating, at the meso-scale, how fluid viscosity and flow rate govern the migration and clogging behavior of fines, thereby explaining the phenomena observed in the experiments and exploring the optimization space for fluid parameters beyond the experimental scope. These two parts are complementary, jointly constructing a comprehensive understanding from macroscopic phenomena to microscopic mechanisms. These findings offer valuable insights for optimizing hydraulic fracturing treatment designs for unconsolidated sandstone reservoirs.

The overall research framework, integrating experimental and numerical approaches to elucidate these synergistic mechanisms, is illustrated schematically in Fig. 1.

Figure 1: Schematic of the research framework

2  Materials and Experiment Procedures

2.1 Sample Preparation and Materials

Artificial unconsolidated sandstone slabs were prepared for this study to simulate the target reservoir rocks. The base composition of the sandstone material—by mass percentage—was 60% quartz sand, 24% illite, 8% plagioclase, and 6% montmorillonite. To simulate variations in rock strength observed in real reservoirs, different amounts of cement were added during the sample preparation process for this purpose. Higher cement content led to stronger inter-particle bonding, and consequently, higher compressive strength of the sandstone slabs. The proppant used in the experiments included three sets of mesh sizes (and their corresponding particle size ranges) to systematically investigate size-dependent embedment behavior: 20/40, 30/50, and 40/70 mesh. The carrier fluid used in the experiments was 2 wt% KCl brine. This design strategy leverages simulation to complement experimental limitations: the 1 mPa·s case anchors to the experimental fluid’s viscosity, while 10–20 mPa·s covers viscosified fracturing fluids commonly used for proppant transport in the field. By correlating experimental results (with field-relevant low-viscosity fluid) to simulation outcomes (with varied viscosities), the study effectively characterizes how fluid viscosity regulates fines migration and conductivity loss—without requiring additional physical experiments. Fines used for the migration studies were sourced from actual unconsolidated sandstone formations to ensure their representativeness.

2.2 Experimental Setup

The schematic diagram of the customized experimental apparatus is presented in Fig. 2. The core component was a conductivity cell designed to hold the artificial sandstone slabs and proppant pack, which complies with relevant API standards for fracture conductivity testing. A distinctive feature of the experimental setup was an intermediate container installed at the inlet of the conductivity cell. This container—equipped with a bottom stirrer—enabled precise mixing of fines with fracturing fluid at predetermined concentrations prior to their injection into the fracture. A hydraulic loading jack was employed to apply and maintain precise closure stress on the sample. Fluid injection was controlled by a precision pump with a flow rate accuracy of 0.01 mL/min, and the system was equipped with precise back-pressure control functionality. Pressure differentials across the fracture were monitored via high-accuracy pressure transducers throughout the experiments. The baseline behaviors of pure proppant packs without fines and the embedment characteristics of pure unconsolidated formation rock are well-documented in the literature. For instance, the work of Zhang et al. [32] provides fundamental data on the conductivity of pure proppant packs under various stresses, which serves as a reference for isolating the effect of fines migration in our study. Chen et al. [33] provides a well-validated reference for understanding the baseline behavior of a pure proppant pack. The experimental parameters in this study, including proppant sizes and rock strength, were designed to be consistent with these established benchmarks, allowing for a meaningful comparative analysis where the synergistic effects are evaluated against these known baseline behaviors.

Figure 2: Schematic of conductivity measurement setup and unconsolidated sand slab

2.3 Testing Procedure

The experimental procedure commenced with the assembly of the artificial sandstone slabs and the proppant pack (typically at a concentration of 5 kg/m3 of proppant) inside the conductivity cell. The desired closure stress was initially applied and maintained via the hydraulic loading jack. The carrier fluid (2 wt% KCl brine) was then injected into the system at a constant flow rate (e.g., 200 mL/min). For experiments focused on fines migration, fines were mixed into the carrier fluid inside the intermediate container at predetermined mass concentrations (e.g., 2%, 5%, 10%) via the stirrer. Fluid flow was temporarily stopped once a stable pressure differential between the fracture inlet and outlet was established. The fracture conductivity was determined at each closure stress level for a minimum of one hour at the target formation temperature, which was maintained, to ensure data stability. Subsequently, the closure stress was increased stepwise from 3 to 27 MPa with an increment of 6 MPa, and the measurement process was repeated at each stress level to obtain the fracture conductivity decline curve under progressively increasing effective stress. A separate set of experiments employed a sand-packed model in which fines were initially placed at the chamber bottom instead of being pre-mixed, enabling the study of fines mobilization and invasion into the proppant pack when the fluid velocity exceeded the critical lift velocity of fines, thereby enabling the isolation of the clogging mechanism from proppant embedment.

2.4 Data Processing

The hydraulic fracture conductivity (Kf) was calculated using Darcy’s law for linear flow through a fracture:

Kf=QμLwΔP(1)

where: Q is the volumetric flow rate (m3/s), μ is the fluid viscosity (Pa·s), L is the length of the fracture (m), w is the fracture width (m), ΔP is the pressure differential across the fracture (Pa).

Key parameters used in the experiments and subsequent modeling are summarized in Table 1.

3  Experimental Results

This section presents the experimental results of this study focused on investigating the decline of hydraulic fracture conductivity under the combined effects of fines migration and proppant embedment. The findings are systematically presented to demonstrate the impact of key variables—including fines concentration, closure pressure, and proppant properties—on fracture conductivity.

3.1 Impact of Fines Migration on Fracture Conductivity

The aforementioned experimental results indicate that particle migration and proppant embedment are two critical mechanisms responsible for the abrupt decline in fracture conductivity of unconsolidated sandstone. The experiments have established the macroscopic influences of parameters including particle concentration and closure pressure; however, a mesoscopic-scale mechanistic elucidation is lacking regarding how fluid properties (e.g., viscosity) and operating parameters (e.g., injection rate) regulate the dynamic transport behavior of particles within fractures and the resultant plugging patterns. As shown in Fig. 3, experiments observed nearly complete loss of conductivity at a closure pressure of 27 MPa and a particle concentration of 10%; however, does this critical concentration vary with the fluid’s carrying capacity? To address such questions and gain in-depth insight into the fluid dynamic mechanisms underlying the plugging phenomena observed in experiments, numerical simulation tools capable of capturing particle-fluid interactions are required. Thus, this study will subsequently employ CFD-DEM coupled simulation to focus on investigating how two parameters—fluid viscosity and injection rate, which were fixed in experiments yet critical in field applications—finely regulate the migration scope and plugging configurations of particles. This will complement the findings of experimental studies and provide a more comprehensive foundation for optimizing fracturing treatment parameters.

Figure 3: The effect of fines concentration and closure pressure on fracture conductivity

3.2 Proppant Embedment Behavior

To ensure the representativeness and engineering applicability of the proppant elastic modulus values selected in this study, the moduli of 1000, 1500, 2000, 2500, and 3000 MPa were determined based on two key considerations: first, commercially available proppant specifications (e.g., ceramic proppants and resin-coated sand commonly used in unconsolidated sandstone fracturing operations, whose elastic moduli typically range from 800 to 3500 MPa according to industry standards such as API RP 60 and supplier technical data sheets); second, targeted screening of industrial products—proppant samples with the desired modulus values were sourced from leading proppant manufacturers, and their actual elastic moduli were verified using a nanoindentation tester following the ISO 14577 standard. The nanoindentation tests were performed on three randomly selected particles from each proppant batch, with 5 indentation points per particle, ensuring the measured modulus values deviated by no more than ±5% from the target values. This selection strategy balances the needs of experimental design (clear differentiation of modulus effects) and field practicality (using industrially viable proppants), avoiding arbitrary modulus assignments.

The embedment behavior of proppants into the artificial unconsolidated sandstone slabs was systematically assessed in this study. All embedment tests were conducted under the following controlled conditions to ensure consistency with target reservoir environments: the test temperature was maintained at 26°C, and the artificial sandstone slabs used had a compressive strength of 15 MPa. During testing, the closure pressure was stepwise increased from 0 to 35 MPa at an increment of 5 MPa, and each pressure level was held for 30 min to allow for stable proppant embedment (avoiding transient deformation effects). The proppant elastic modulus was pre-calibrated using a nanoindentation tester in accordance with the ISO 14577 standard, with a test deviation controlled within ±5%; proppant particle sizes (20/40, 30/50, 40/70 mesh) were verified using standard sieves, ensuring a particle size distribution deviation of ≤±3%. After each pressure stabilization period, the proppant embedment depth was measured using a laser scanning confocal microscope with a resolution of 0.1 μm—5 random observation points were selected for each sample, and the results were averaged to reduce random measurement errors; additionally, 3 parallel samples were tested for each proppant type, with a repeatability error controlled within 8%. Figs. 4 and 5 demonstrate the effects of closure pressure, proppant elastic modulus, and proppant size on proppant embedment depth. As shown in Fig. 4, the embedment depth consistently increased with increasing applied closure pressure for proppants of any elastic modulus or size, as the applied stress is the key primary driving force for embedment. Notably, under the same closure pressure, proppants with a higher elastic modulus (i.e., stiffer proppants) displayed a greater embedment depth. For instance, at 35 MPa, the embedment depth for proppants with a modulus of 3000 MPa was approximately 0.04 mm, in contrast to only 0.02 mm for proppants with a modulus of 1000 MPa. This is ascribed to the lower elastic deformation capacity of stiffer proppants, which results in more effective stress transfer to the surrounding rock matrix. Fig. 5 presents a comparative analysis of the embedment behavior of proppants with different mesh sizes (40/70, 30/50, 20/40) under a range of closure stresses. The results demonstrate that under identical closure pressure, larger-sized proppants (20/40 mesh) display a reduced embedment depth. This behavior is due to their reduced contact area with the rock surface, which facilitates more effective distribution of the applied stress, thereby alleviating the localized stress concentrations that drive embedment.

Figure 4: The effect of closure pressure on embedment depth

Figure 5: The effect of proppant size on embedment depth

3.3 Synergistic Effect of Fines Concentration and Filtration

The synergistic effect of fines mass concentration and fluid filtration process on hydraulic fracture conductivity over time is systematically presented in Fig. 6 of this study. The experiments utilized a 20/40 mesh proppant pack under a constant closure pressure of 15 MPa, with a uniform proppant concentration of 5 kg/m3. The carrier fluid was 2 wt% KCl brine injected at a steady rate of 200 mL/min. Fines were pre-mixed into the fluid at specified mass concentrations (1%, 2%, 3%, 4%) using the intermediate stirrer container. The fracture conductivity was monitored until stabilization at each condition, ensuring the measured values reflected a steady-state impairment. The results demonstrate that hydraulic fracture conductivity underwent a pronounced decline as the applied filtration coefficient increased, a trend further amplified at higher fines mass concentrations. For example, at an applied filtration coefficient of 35 m−1, the hydraulic fracture conductivity was sustained at approximately 150 μm2·cm at a 1% fines mass concentration, but it dropped to about 80 μm2·cm with the fines mass concentration increased to 4%. This highlights the critical role of fines retention and accumulation within the proppant pack in accelerating hydraulic fracture conductivity impairment during reservoir fluid flow.

Figure 6: The effect of filtration coefficient on fracture conductivity

4  Numerical Modeling Principles

The detachment of fines induced by sustained reservoir fluid flow can result in a decline in hydraulic fracture permeability. Many researchers have investigated the mechanisms of fracture permeability damage caused by fines in detail, yet the grain-scale mechanisms still remain not yet fully elucidated. As shown in Fig. 7, the typical process of fines migration during production involves several key stages from inside the hydraulic fracture to the wellbore. While this study primarily focuses on fines originating from the formation rock (as simulated in our experiments), it is important to note that the fine particles depicted in this process could also stem from the crushing of proppants under sufficient closure stress, as investigated in other reservoir types [34]. The migration and clogging behavior of these crush-generated fines would similarly impair fracture conductivity.

Figure 7: The stages of fine migration from hydraulic fracture to wellbore

4.1 Fluid Phase Governing Equations

The fluid has independent continuous physical field parameters, including velocity field, concentration field, and pressure field. The influence of proppants relative to the fluid phase is primarily characterized by the fluid volume fraction αf and the interphase momentum exchange source term Mfs between the fluid and proppants.

The fluid continuity equation,

∂∂t(αfρf)+∇⋅(αfρfuf)=0(2)

wherein, αf denotes the fluid volume fraction; ρf denotes the fluid density, kg·m−3; and uf denotes the fluid velocity, m·s−1.

The fluid momentum equation,

∂∂t(αfρfuf)+∇⋅(αfρfufuf)=−αf∇p+∇⋅(αfτf)−αfρfg+Mfs(3)

wherein, p denotes the fluid pressure, Pa; τf denotes the fluid viscous stress tensor, N·m−2; t denotes the time, s; g denotes the gravitational acceleration, m·s−2; and Mfs denotes the interphase momentum exchange source term between the fluid and proppants, which includes the collision force between solid particles and the interphase momentum transfer term.

τf can be expressed as follows:

τf=μf[(∇uf)+(∇uf)T](4)

wherein, μf denotes the liquid-phase dynamic viscosity, Pa·s.

For the fluid phase, the smallest mesh generated in the model corresponds to the fluid computational control volume element, in which αf and Mfs are physical quantities associated with proppants.

The fluid volume fraction is determined using the point-in-cell method: points are sampled within the computational grid cell; if a sampling point lies within the surface boundary of solid phase particles, it is excluded from the count. This process is applied to all points in the computational grid cell, and the proportion of fluid within the computational grid cell is determined using the following formula, so as to derive the fluid volume fraction.

αf=nfN(5)

wherein, nf denotes the number of fluid points within the computational grid cell; and N denotes the total number of points within the computational grid cell.

With the continuous injection of proppant-carrying fluid, proppants deposit at the bottom and form a proppant bank. This proppant bank reduces the available space for the free flow of proppant-carrying fluid, increases the flow velocity in the fracture, and the surface inhomogeneity of particles accumulated in the proppant bank also induces significant turbulence. The turbulence intensity is employed to characterize the intensity of turbulence fluctuations, with the formula given below,

I=uu¯=0.16(Re)−1/8(6)

Re=ρfULμf(7)

wherein, U denotes the cross-sectional average velocity of the fluid, m·s−1; and L denotes the characteristic length, m. The proppant transport process belongs to the category of fracture flow.

The characteristic length L is defined as:

L=wH2(w+H)(8)

wherein, w denotes the fracture width, m; and H denotes the fracture height, m.

During turbulent flow, intense energy and momentum exchange takes place between proppants and fracturing fluid in the hydraulic fracture. Turbulent numerical simulation for hydraulic fractures primarily consists of two types: transient simulation and steady-state simulation. Transient simulation directly solves the instantaneous governing equations for turbulence, yet it is only applicable to cases involving low Reynolds numbers and relatively simple boundary conditions. In contrast, steady-state simulation simplifies and approximates the characteristics of turbulence, thereby achieving faster computational speed. The transient simulation method was employed for calculations in this study, and the k-ε model is typically adopted as the turbulence model for particle-fluid two-phase flow systems. The Standard k-ε model is widely used, with moderate computational cost and computational accuracy sufficient to meet the requirements of general engineering calculations.

The turbulent kinetic energy and diffusion equation,

∂∂t(ρk)+∂∂xi(ρkui)=∂∂xj[(μ+μtσk)∂k∂xj]+Gk−ρε(9)

∂∂t(ρε)+∂∂xi(ρεui)=∂∂xj[(μ+μtσε)∂ε∂xj]+C1εεkGk−C2ερε2k(10)

wherein, k denotes the turbulent kinetic energy, m2·s−2; ε denotes the turbulent dissipation rate, m2·s−3; σk and σε denote respectively the turbulent Prandtl numbers corresponding to turbulent kinetic energy k and dissipation rate ε; C1ε and C2ε denote empirical constants; Gk and Gb denote respectively the production terms of turbulent kinetic energy induced by average velocity and buoyancy, kg·m−2·s2; μ denotes the dynamic viscosity of the fluid, Pa·s; and μt denotes the viscosity increment due to turbulence, Pa·s.

4.2 Particle Phase Governing Equations

The proppant phase in the proppant-carrying fluid is a discontinuous phase. The Discrete Element Method (DEM) is a numerical method applicable to addressing problems in discontinuous medium mechanics, first proposed by Cundall and Strack [35]. In DEM, the motion of each particle is governed by Newton’s second law of motion; detailed updates to displacement and angular velocity are made at each time step, and the resultant force and resultant moment of particle collisions are computed using the soft-sphere model. The motion equation for a single particle is given by:

mpdupdt=FC+FA+FV(11)

Ipcdωpdt=Tpc(12)

wherein, mp denotes the particle mass, kg; up denotes the particle linear velocity, m·s−1; Ip denotes the particle moment of inertia, kg·m2; ωp denotes the particle angular velocity, rad·s−1; and Tp denotes the contact torque generated by the contact between the particle and other particles, N·m.

The detachment of fines induced by reservoir fluid flow can cause a decline in hydraulic fracture permeability. Many researchers have studied the mechanisms of fracture permeability damage caused by fines in detail, though grain-scale mechanisms still remain not yet fully elucidated.

5  Numerical Modeling: Setup, Parameters and Results Analysis

5.1 Model Setup and Parameter Configuration

To investigate the synergistic effect of proppant embedment and fines migration, a high-fidelity CFD-DEM coupled model was developed. The model geometry was designed to rigorously replicate the experimental apparatus, representing a rectangular fracture channel with a length of 10 cm, height of 5 cm, and a consistent width (thickness) of 0.4 cm, bounded by argillaceous unconsolidated sandstone slabs (see Fig. 8).

Figure 8: Schematic of the geometry and observation plane of the CFD-DEM model

The computational domain was discretized using a structured mesh with a minimum cell size of 0.1 mm, determined through a mesh independence study to ensure solution accuracy without excessive computational cost. The boundary conditions were defined to mimic the experimental flow regime: a velocity inlet for specified flow rates, a pressure outlet, and no-slip walls for the fracture boundaries. The proppant pack was initialized as a random assembly of spherical particles with size distributions corresponding to the experimental 20/40, 30/50, and 40/70 mesh sizes. The fluid-proppant interaction was resolved via a two-way coupled Eulerian-Lagrangian approach, where the fluid phase solves the Navier-Stokes equations and the particle motion is tracked using the Discrete Element Method (DEM). All other key parameters, such as fluid properties and proppant densities, were maintained strictly consistent with the experimental conditions detailed in Section 2.4 and Table 1, thereby ensuring the validity and comparability of the simulation results.

5.2 Results and Analysis of Numerical Simulation

Building upon the experimental observations in Section 3, it is established that fines migration severely impairs fracture conductivity, and the degree of impairment is significantly influenced by fines concentration and closure pressure. However, the fluid viscosity (1 mPa·s) and injection rate (200 mL/min) used in the experiments were fixed. To uncover the intrinsic controlling mechanisms of fluid properties and flow conditions on the key process of fines migration, this numerical simulation section systematically varies the fluid viscosity and injection rate. The simulation results not only explain the clogging phenomena observed in the experiments (e.g., pore-throat plugging induced by high fines concentration) at the meso-scale but, more importantly, quantify how fluid parameters affect the migration front morphology and clogging location of fines. This provides theoretical support for mitigating conductivity decline by optimizing fluid parameters in field operations.

To elucidate the intrinsic mechanisms of fine particle blockage, advanced visualization techniques were employed, including particle tracking, velocity vector fields, and volume rendering of particle concentration. These methods enabled the identification of critical clogging zones and the dynamic evolution of pore-throat blockage. For instance, particle tracking revealed that fines initially migrate along high-velocity paths but deposit in low-velocity regions due to drag reduction and gravitational settling. The visualization results further demonstrated that clogging is initiated by particle bridging at pore throats, where the ratio of particle size to throat diameter exceeds 0.1, leading to rapid permeability reduction. Moreover, the spatial distribution of fines accumulation was correlated with local turbulence intensity, indicating that turbulent fluctuations enhance particle-particle collisions and agglomeration. This in-depth analysis provides a meso-scale understanding of how fluid parameters and particle properties collectively govern the clogging behavior, complementing the macroscopic experimental observations.

5.2.1 Effect of Fluid Viscosity on Fine Migration and Fracture Conductivity

The CFD-DEM coupled model developed in this study—specifically tailored for hydraulic fractures in unconsolidated sandstones—is designed to reveal the mechanism of hydraulic fracture conductivity loss under the synergistic effect of fine particle migration and proppant embedment. 1-s fine particle migration simulations were conducted under different fluid viscosities (1, 10, and 20 mPa·s). Results indicate that fluid viscosity exerts a distinct regular influence on the migration characteristics of fine particles and hydraulic fracture conductivity, with its underlying mechanism closely linked to fluid particle-carrying capacity, flow field properties, and fine particle-proppant interactions. At a fluid viscosity of 1 mPa·s, the drag force exerted by the low-viscosity fluid on fine particles is insufficient to overcome the self-gravity of fine particles and the weak adhesive forces between fine particles and proppants in the fracture inlet region. Consequently, fine particles exhibit a short overall migration distance, with a large quantity accumulating at the fracture inlet. Meanwhile, under low fluid viscosity conditions, fluid flow more closely approximates uniform laminar flow. Fine particles experience highly consistent forces within the flow field; thus, the leading edge of fine particles invading the hydraulic fracture forms a regular “simultaneous advancement” pattern (see Fig. 9a). As the fluid viscosity increases to 10 mPa·s, the fluid viscosity-induced force is markedly enhanced, and the drag force can drive fine particles to migrate toward the deep section of the hydraulic fracture. Additionally, the turbulence intensity of the flow field increases, resulting in unbalanced forces acting on fine particles within the non-uniform flow field. Consequently, the leading edge of fine particle invasion gradually transitions from the “simultaneous advancement” pattern to a gentle arc shape (see Fig. 9b). At this point, while fine particles begin to penetrate the pores of the proppant pack, large-scale clogging has not yet developed. As the fluid viscosity further increases to 20 mPa·s, the drag force of the high-viscosity fluid attains a maximum value, enabling it to drive fine particles to migrate to a position corresponding to 50% of the hydraulic fracture length. During this process, fine particles undergo significant sieving and deposition within the proppant pack, leading to an increased clogging rate of pore throats in the proppant pack. In accordance with Darcy’s Law, the hydraulic fracture permeability decreases accordingly, ultimately leading to a significant reduction in hydraulic fracture conductivity. Furthermore, under high fluid viscosity conditions, the velocity gradient effect of the flow field is more pronounced, the variation in the spatial distribution of fine particles increases, and the arc-shaped feature of their invasion leading edge becomes increasingly distinct (see Fig. 9c). This confirms that fluid viscosity—by regulating the migration and clogging behaviors of fine particles—acts as a key factor affecting the long-term conductive performance of hydraulic fractures in unconsolidated sandstones.

Figure 9: Fine particle migration characteristics in hydraulic fractures of unconsolidated sandstone at different fluid viscosities: (a) Viscosity of 1 mPa·s; (b) Viscosity of 10 mPa·s; (c) Viscosity of 20 mPa·s

5.2.2 Effect of Fluid Displacement Rate on Fine Migration and Fracture Conductivity

Based on the CFD-DEM coupled model established in this study specifically for hydraulic fractures in unconsolidated sandstone, under the premise of maintaining consistent parameters—including proppant size (200~300 μm), proppant placement thickness (0.4 cm), initial fine particle concentration (4%), and closure pressure (15 MPa)—fine particle migration simulations of 1 s duration were conducted by varying only the fluid displacement. Due to differences in fluid drag force and flow regime, different displacements exert a distinct differential impact on the migration behavior of fine particles and hydraulic fracture conductivity: When the displacement is 100 mL/min, the fluid within the hydraulic fracture exists in a stable laminar flow regime. At this point, the fluid drag force can only overcome 60%~70% of the self-gravity of individual fine particles, failing to drive fine particles to migrate toward the deep section of the hydraulic fracture. This results in the deposition of a large quantity of fine particles at the hydraulic fracture inlet, forming a “shallow plugging zone”. The hydraulic fracture permeability decreases by only approximately 35%, and the low flow rate results in low construction efficiency, as illustrated in Fig. 10a; When the displacement is increased to 200 mL/min, the flow field transitions to a weak turbulent regime. The moderate drag force enables fine particles to migrate effectively to 50%~60% of the hydraulic fracture length. A portion of fine particles can penetrate into the middle part of the proppant pack via the “throat screening” effect of proppant pores. Furthermore, weak turbulence only results in the dispersed accumulation of the fine particle migration front in a “gentle arc” shape, with no continuous plugging zone formed within the proppant pack. Eventually, the hydraulic fracture permeability decreases by 55%~60%, which not only avoids excessive fine particle accumulation at the inlet but also prevents deep plugging. This aligns with the field practical principle of “balancing construction efficiency and hydraulic fracture damage” in the field fracturing of unconsolidated sandstone reservoirs, as illustrated in Fig. 10b; When the displacement reaches 300 mL/min, the flow enters a regime of strong turbulence. The excessive drag force not only drives fine particles to migrate to over 80% of the hydraulic fracture length, but the “vortex effect” of strong turbulence also causes fine particles to form a “continuous plugging zone” in the deep section of the hydraulic fracture. Meanwhile, the scouring action of the high flow rate on the proppant induces the disordered arrangement of proppant particles, which further intensifies the “bridging plugging” of fine particles. Eventually, the hydraulic fracture permeability decreases by more than 80%, and the high flow rate tends to exacerbate “sand production” in unconsolidated sandstone reservoirs. Consequently, this displacement is rarely employed in field fracturing operations for unconsolidated sandstone reservoirs, as illustrated in Fig. 10c.

Figure 10: Fine particle migration characteristics in hydraulic fractures of unconsolidated sandstone at different fluid displacements: (a) 100 mL/min; (b) 200 mL/min; (c) 300 mL/min

6  Discussions

This study systematically reveals the mechanisms of hydraulic fracture conductivity decline in unconsolidated sandstone reservoirs by integrating experimental measurements and numerical simulations. The experimental results demonstrate that fines migration and proppant embedment are two primary, mutually reinforcing causes of conductivity loss (Figs. 3–5). The numerical simulations further elucidate how fluid properties (viscosity) and operational conditions (displacement rate), as key regulatory means, exacerbate or alleviate the aforementioned damage by influencing the migration and clogging patterns of fines (Figs. 9 and 10). First, fines migration and proppant embedment are two primary, mutually reinforcing causes of hydraulic fracture conductivity decline in unconsolidated sands, and their synergistic effect significantly accelerates conductivity loss—especially under 27 MPa closure stress and 10% fines concentration, where conductivity is nearly eliminated. Second, fluid viscosity exerts a critical regulatory effect on fines migration in hydraulic fractures: low-viscosity fluids (1 mPa·s) lack sufficient drag force, causing fines to accumulate near the hydraulic fracture inlet (forming “shallow plugging zones”); high-viscosity fluids (20 mPa·s) drive fines to the fracture deep section, leading to widespread proppant pack pore-throat clogging. Third, proppant properties directly influence proppant embedment in unconsolidated sands—especially under low-to-medium closure stress (15–35 MPa): stiffer proppants (higher elastic modulus, e.g., 3000 MPa) exhibit greater embedment resistance than softer ones (1000 MPa); larger proppants (200–425 μm) reduce embedment via smaller rock matrix contact areas, which disperse applied stress.

For field fracturing operations in unconsolidated sandstone reservoirs, two key optimization strategies are recommended to mitigate hydraulic fracture conductivity loss: (1) Adjust fluid viscosity to 10–15 mPa·s to balance fines control and flow efficiency; (2) Select proppants with elastic modulus of 2000–3000 MPa and particle size of 200–250 μm to reduce embedment while avoiding excessive pore blocking. In addition, future research should focus more on the quantitative characterization of the synergistic effect of fines migration and proppant embedment—specifically, exploring the quantitative relationship between closure stress (5–40 MPa) and synergy intensity, and developing a predictive model for conductivity decline based on this relationship.

7  Conclusion

This study systematically investigates the synergistic mechanisms of fines migration and proppant embedment on hydraulic fracture conductivity loss in unconsolidated sandstone reservoirs through integrated experimental and CFD-DEM numerical approaches. The key findings are summarized as follows:

(1)   A coupled model incorporating both fines migration and proppant embedment was developed and validated against experimental data, effectively predicting the long-term conductivity decline in unconsolidated sandstone reservoirs. The model addresses the critical research gap concerning the synergistic effect of these two mechanisms, which has been largely overlooked in previous studies.

(2)   Both fines migration and proppant embedment significantly impair fracture conductivity. High fines concentrations (e.g., 10%) lead to severe pore-throat clogging, while increasing closure stress amplifies embedment, especially in low- to medium-stress regimes (15–35 MPa). Stiffer proppants (3000 MPa) exhibit twice the embedment depth of softer ones (1000 MPa), and larger proppant sizes effectively reduce embedment by distributing stress more efficiently.

(3)   Fluid parameters play a decisive role in regulating fines migration and associated conductivity damage. An optimal fluid viscosity of 10 mPa·s combined with a displacement rate of 200 mL/min promotes uniform fines distribution without causing deep clogging or excessive inlet accumulation, thereby balancing fracture damage and operational efficiency.

The proposed model and parameter optimization strategy provide practical guidance for fracturing design in unconsolidated and hydrate-bearing reservoirs, supporting long-term productivity assurance. Future work should focus on the dynamic effects of reservoir creep and clay content to further refine the predictive capability of the model.

Acknowledgement: Not applicable.

Funding Statement: This work had been financially supported by the National Natural Science Foundation of China (51704324, U1762213), Natural Science Foundation of Shandong Province (ZR2022ME068) and the Government of Perm Krai, research project No. CЭд-26-08-08-32 from 25.01.2024.

Author Contributions: Xian Shi: Conceptualization, Methodology, Writing—Original Draft. Botao Zhang: Investigation, Data Curation. Weidong Zhang: Validation, Formal Analysis. Zenghua Ma: Software. Bo Zhang: Project Administration. Ahmad Ramezanzadeh: Modeling, Simulation. Bin Li: Validation. Jian Mao: Data Curation. All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials: The data and models used in this study are available from the corresponding author upon reasonable request.

Ethics Approval: This article does not contain any studies with human participants or animals performed by any of the authors.

Conflicts of Interest: The authors declare no conflict of interest to report regarding the present study.

References

1. Bai T, Chen Z, Aminossadati SM, Li L, Liu J, Lu H. Dimensional analysis and prediction of coal fines generation under two-phase flow conditions. Fuel. 2017;194:460–79. doi:10.1016/j.fuel.2017.01.040. [Google Scholar] [CrossRef]

2. Bedrikovetsky P, Vaz A, Machado F, Zeinijahromi A, Borazjani S. Skin due to fines mobilization, migration, and straining during steady-state oil production. Petrol Sci Technol. 2012;30(15):1539–47. doi:10.1080/10916466.2011.653702. [Google Scholar] [CrossRef]

3. Guo Z, Hussain F, Cinar Y. Physical and analytical modelling of permeability damage in bituminous coal caused by fines migration during water production. J Nat Gas Sci Eng. 2016;35:331–46. doi:10.1016/j.jngse.2016.08.031. [Google Scholar] [CrossRef]

4. Yan C, Chen Y, Chen T, Cheng Y, Yan X. Experimental study of hydraulic fracturing for unconsolidated reservoirs. Rock Mech Rock Eng. 2022;55(6):3399–424. doi:10.1007/s00603-022-02827-6. [Google Scholar] [CrossRef]

5. Yang Y, Li X, Ju Y, Li G, Li X. Multifactor coupled evaluation model of fracability based on a combinatorial weighting method: a case study of the low-permeability sandstone reservoirs of Yanchang Formation, southwestern Ordos Basin, China. Rock Mech Rock Eng. 2025;58(1):905–25. doi:10.1007/s00603-024-04193-x. [Google Scholar] [CrossRef]

6. Lehman LV, Blauch ME, Robert LM. Desorption enhancement in fracture-stimulated coalbed methane wells. In: Proceedings of the SPE Eastern Regional Meeting; 1998 Nov 9–11; Pittsburgh, Pennsylvania. doi:10.2118/51063-ms. [Google Scholar] [CrossRef]

7. Zou YS, Mang XF, Wang L, Lin X. Experimental evaluation of conductivity of fracturing in medium and high-rank coalbed. J China Coal Soc. 2011;36(3):473–6. [Google Scholar]

8. Bedrikovetsky P, Siqueira FD, Furtado CA, Souza ALS. Modified particle detachment model for colloidal transport in porous media. Transp Porous Medium. 2011;86(2):353–83. doi:10.1007/s11242-010-9626-4. [Google Scholar] [CrossRef]

9. Fu J, Jiang S, Zhang L, Xiong P, Zhang K. An improved fractional damage creep model for shale and its application to the long-term fracture conductivity prediction. Comput Geotech. 2025;188:107601. doi:10.1016/j.compgeo.2025.107601. [Google Scholar] [CrossRef]

10. Sun Y, Zeng S, Dai C, Ding X, Ding F, Wang C, et al. Artificial fractures conductivity impairment in gas shale when fracturing with surface water. J Hydrol. 2024;629:130557. doi:10.1016/j.jhydrol.2023.130557. [Google Scholar] [CrossRef]

11. Liu YX, Zhou HY, Guo JC, Zhang Q, Chen C. Controlling factors of fracture aperture reduction based on experimental study using proppant pack apparent Young’s modulus. J Petrol Sci Eng. 2022;208:109506. doi:10.1016/j.petrol.2021.109506. [Google Scholar] [CrossRef]

12. Zhang J, Kamenov A, Zhu D, Hill AD. Development of new testing procedures to measure propped fracture conductivity considering water damage in clay-rich shale reservoirs: an example of the Barnett Shale. J Petrol Sci Eng. 2015;135:352–9. doi:10.1016/j.petrol.2015.09.025. [Google Scholar] [CrossRef]

13. Ahamed MAA, Perera MSA, Ranjith PG. Implementation of an elastoplastic constitutive model to study the proppant embedment in coal under different pore fluid saturation conditions: a numerical and experimental study. Fuel. 2022;317:123488. doi:10.1016/j.fuel.2022.123488. [Google Scholar] [CrossRef]

14. Mojarad RS, Settari A. Coupled numerical simulation of reservoir flow with formation plugging. In: Proceedings of the Canadian International Petroleum Conference; 2005 Jun 7–9; Calgary, Alberta. doi:10.2118/2005-061. [Google Scholar] [CrossRef]

15. Yan X, Huang Z, Yao J, Song W, Li Y, Gong L. Theoretical analysis of fracture conductivity created by the channel fracturing technique. J Nat Gas Sci Eng. 2016;31:320–30. doi:10.1016/j.jngse.2016.03.038. [Google Scholar] [CrossRef]

16. Zheng W, Yu H, Yang S, Cui X, Zhang H. Estimation of proppant embedment from laboratory instrumented indentation tests on Montney formation in Western Canada. Int J Rock Mech Min Sci. 2023;170:105495. doi:10.1016/j.ijrmms.2023.105495. [Google Scholar] [CrossRef]

17. Guerra J, Zhu D, Hill AD. Impairment of fracture conductivity in the eagle ford shale formation. SPE Prod Oper. 2018;33(4):637–53. doi:10.2118/184857-pa. [Google Scholar] [CrossRef]

18. Tang Y, Ranjith PG, Wu B. Experimental study of effects of shearing on proppant embedment behaviour of tight gas sandstone reservoirs. J Petrol Sci Eng. 2019;172:228–46. doi:10.1016/j.petrol.2018.07.066. [Google Scholar] [CrossRef]

19. Xu J, Ding Y, Yang L, Liu Z, Gao R, Yang H, et al. Conductivity analysis of tortuous fractures filled with non-spherical proppants. J Petrol Sci Eng. 2021;198:108235. doi:10.1016/j.petrol.2020.108235. [Google Scholar] [CrossRef]

20. He B, Xie L, Che Y, Zhang Y, Liu J. A simplified analytical elastoplastic embedment depth model of proppant based on the Drucker–Prager yield criterion. Rock Mech Rock Eng. 2023;56(9):6207–17. doi:10.1007/s00603-023-03394-0. [Google Scholar] [CrossRef]

21. Zhang F, Wang T, Liu F, Peng M, Bate B, Wang P. Hydro-mechanical coupled analysis of near-wellbore fines migration from unconsolidated reservoirs. Acta Geotech. 2022;17(8):3535–51. doi:10.1007/s11440-021-01396-2. [Google Scholar] [CrossRef]

22. Mittal A, Rai CS, Sondergeld CH. Proppant-conductivity testing under simulated reservoir conditions: impact of crushing, embedment, and diagenesis on long-term production in shales. SPE J. 2018;23(4):1304–15. doi:10.2118/191124-pa. [Google Scholar] [CrossRef]

23. Li B, Shen Y, Sun Y, Qi Y, Shan H, Zhang G. Assessing hydraulic fracturing feasibility in marine hydrate reservoirs: impact of closure pressure and hydrate decomposition on fracture conductivity. Ocean Eng. 2024;309:118569. doi:10.1016/j.oceaneng.2024.118569. [Google Scholar] [CrossRef]

24. Zou YS, Zhang SC, Zhang J. Experimental method to simulate coal fines migration and coal fines aggregation prevention in the hydraulic fracture. Transp Porous Medium. 2014;101(1):17–34. doi:10.1007/s11242-013-0228-9. [Google Scholar] [CrossRef]

25. Mitchell TR, Leonardi CR. Micromechanical investigation of fines liberation and transport during coal seam dewatering. J Nat Gas Sci Eng. 2016;35:1101–20. doi:10.1016/j.jngse.2016.09.038. [Google Scholar] [CrossRef]

26. Ahamed MAA, Perera MSA, Li DY, Ranjith PG, Matthai SK. Proppant damage mechanisms in coal seam reservoirs during the hydraulic fracturing process: a review. Fuel. 2019;253:615–29. doi:10.1016/j.fuel.2019.04.166. [Google Scholar] [CrossRef]

27. Dou X, Liu Z, Yang D, Zhao Y, Li Y, Gao D, et al. 3D CFD-DEM modeling of sand production and reservoir compaction in gas hydrate-bearing sediments with gravel packing well completion. Comput Geotech. 2025;177:106870. doi:10.1016/j.compgeo.2024.106870. [Google Scholar] [CrossRef]

28. Lei J, Chen Z, Zhao J, Wang Y, Guo W. Effect of decrease in pore water salinity induced by hydrate decomposition in clayey silt sediment on the destabilization-migration-clogging characteristics of illite. Energy. 2024;308:132887. doi:10.1016/j.energy.2024.132887. [Google Scholar] [CrossRef]

29. Benge M, Katende A, Rutqvist J, Radonjic M, Bunger A. Creep properties of shale and predicted impact on proppant embedment for the caney shale, Oklahoma. Rock Mech Rock Eng. 2023;56(8):5903–21. doi:10.1007/s00603-023-03362-8. [Google Scholar] [CrossRef]

30. Peng Y, Luo A, Li Y, Wu Y, Xu W, Sepehrnoori K. Fractional model for simulating Long-Term fracture conductivity decay of shale gas and its influences on the well production. Fuel. 2023;351:129052. doi:10.1016/j.fuel.2023.129052. [Google Scholar] [CrossRef]

31. Doan T, Indraratna B, Nguyen TT, Rujikiatkamjorn C. Coupled CFD-DEM modelling of clogging of granular columns by cohesive fines. Comput Geotech. 2025;177:106902. doi:10.1016/j.compgeo.2024.106902. [Google Scholar] [CrossRef]

32. Zhang J, Ouyang L, Zhu D, Hill AD. Experimental and numerical studies of reduced fracture conductivity due to proppant embedment in the shale reservoir. J Petrol Sci Eng. 2015;130:37–45. doi:10.1016/j.petrol.2015.04.004. [Google Scholar] [CrossRef]

33. Chen D, Ye Z, Pan Z, Zhou Y, Zhang J. A permeability model for the hydraulic fracture filled with proppant packs under combined effect of compaction and embedment. J Petrol Sci Eng. 2017;149:428–35. doi:10.1016/j.petrol.2016.10.045. [Google Scholar] [CrossRef]

34. Yin B, Zhang Y, Lou Y, Liu S. Damage mechanism of proppant and conductivity reduction post fracturing in unconventional reservoirs. Fuel. 2025;380:133086. doi:10.1016/j.fuel.2024.133086. [Google Scholar] [CrossRef]

35. Cundall PA, Strack ODL. Discussion: a discrete numerical model for granular assemblies. Géotechnique. 1980;30(3):331–6. doi:10.1680/geot.1980.30.3.331. [Google Scholar] [CrossRef]