Mechanisms of Proppant Transport in Rough Fractures of Offshore Unconventional Reservoirs: Shale and Tight Sandstone

<jats:p>After hydraulic fracturing, unconventional reservoirs frequently encounter challenges related to limited effective proppant support distance and suboptimal proppant placement. Due to the strong heterogeneity of offshore reservoirs, which causes varying fracture roughnesses depending on different lithologies, a systematic study of the relationship between roughness and proppant transport could optimize operational parameters. This study incorporates the box dimension method for fractal dimension analysis to quantify roughness in auto-correlated Gaussian distributed surfaces created by true triaxial tests. Combined with the numerical analysis of (computational fluid dynamics) CFD-DEM (discrete element method) for bidirectional coupling, the laws of proppant deposition and transport processes within fractures with different roughnesses are obtained through comparative verification simulations. The results show that for rougher fractures of shale, the proppants are transported farther, but at JRC_52, (joint roughness coefficient), where there may be plugging in curved areas, there is a risk of near-well blockages. Compared to the smooth model, fluctuations in JRC_28 (tight sandstone) drastically increase turbulent kinetic energy within the fracture, altering particle transport dynamics. Moreover, smaller proppants (d/w ≤ 0.3) exhibit better transport capacity due to gravity, but the conductivity of the proppant is limited when the particles are too small. A d/w of 0.4 is recommended to guarantee transport capacity and proppant efficiency near the well. Additionally, proppants injected sequentially from small to large in shale fractures offer optimal propping effects, and can take advantage of the better transport capacity of smaller proppants in rough fractures. The large proppant (d/w = 0.8) is primarily deposited by gravity and forms a sloping sand bed, which subsequently ensures the aperture of the fractures. This research provides a fresh perspective on the influence of fracture roughness on proppant transport in offshore unconventional reservoirs and offers valuable considerations for the order of proppant injection.</jats:p>

section_content: Offshore reservoirs are characterized by rapid sedimentary facies changes, complex connectivity, and strong heterogeneity. Hydraulic fracturing offshore is limited by operational conditions, resulting in lower displacement and fracturing fluid viscosity, leading to poor sand transport capability and a high risk of blockage. After fracturing, the morphology and conductivity of fractures depend on the proppant deposition. Unconventional reservoirs, characterized by low porosity and permeability, often require hydraulic fracturing for production enhancement but frequently suffer from suboptimal fracture conductivity [1, 2]. During proppant transport in fractures, various factors affect the sedimentation rate of the proppant, such as the proppant concentration, particle size, shape, fluid rheology, pumping rate, and fracture width. Accordingly, scholars have also performed much research through both laboratory experiments and numerical simulations. However, due to the large differences in scale and the material of experimental instruments and target parameters, a unified standard for experimental instruments and experimental testing procedures has not yet been established. Researchers including Sahai et al. [3] [4] [5] have built physical simulators with multiple characteristics to visualize the proppant transport and deposition processes with different sensitivities. In the experimental process, hydraulic fractures were typically characterized by presetting the hydraulic fracture into two completely smooth parallel plates, which were utilized to study the influence of proppant deposition and accumulation on fracture closure. Extensive experimental simulations have also been conducted to analyze concentric cylinder devices, narrow-slit flow devices in plexiglas and complex fracture models that contain subslits in terms of fracture number, fracture morphology and fracture density [6]. Compared with the actual fracture, the material, morphology and dimensions of the test device involved in the study have certain limitations, which cannot systematically guide the fracturing process design. In contrast to the limitations of experimental design and conditions, numerical simulations are more feasible in proppant transport studies due to their low cost, diverse model forms and different boundary conditions [7, 8]. 

Currently, numerical simulations of liquid-solid biphasic flows are performed by the Euler-Euler and Euler-Lagrange methods [9]. In contrast, Zhang [10] systematically investigated the transport and deposition of multisize proppant in fractures in both vertical and horizontal wells. The study quantitatively characterized the effects of multisize particles on proppant placement relative to uniform-size particles. Based on previous studies, methods using CFD-DEM (computational fluid dynamics-discrete element method) have advantages in modeling particle loading flows more realistically. Meanwhile, the unresolved CFD-DEM has more potential in the discipline of proppant transport modeling than other simulation methods [11, 12]. Lu et al. [13] obtained the causes and mechanisms of uneven proppant distribution and low proppant efficiency in a single simplified fracture by CFD-DEM simulations of variable fracturing fluid viscosity and proppant density. The study demonstrated that the increase in fluid viscosity contributed to an increase in the proppant migration distance and the length of the effective support near the wellbore. 

Fractures are commonly found in hydraulic fracturing operations at different scales of roughness. Previous efforts have exploited the smooth flat plate fracture model, resulting in the effect of oil flow resistance and collision loss of fluid particles on rough walls being ignored, which in turn altered the proppant transport patterns. Fractal theory was first proposed by the American mathematician Benoît B in the 1970s, which was applied to fracture simulation. Barton [14] proposed the (joint roughness coefficient) JRC-JCS (joint compressive strength) model which first used the JRC to describe the degree of structural surface roughness. Raimbay [15] [16] [17] et al. found through horizontal rough fracture proppant transport experiments that fracture wall roughness controls not only the fluid flow path but also the stability of proppant distribution. Huang et al. [18, 19] conducted a series of experiments on sand transport in vertical rough fractures, which showed that fracture wall roughness significantly enhanced the vertical distribution of proppant in fractures. In a recent study, Jun Li et al. [20] modeled a double-sided, equal-width rough fracture wall employing 3D printing, which macroscopically compared the particle deposition of rough walls with a smooth surface. Using CFD-DEM, Zheng et al. [21] derived a fracture opening distribution function and found that local surface roughness, steady over long ranges, inversely affects the permeability of natural fracture fluid in mismatched composite profiles. Tomac I et al. [22] [23] [24] [25] enriched the theoretical models of proppant migration by quantifying flow mechanisms within rough fractures and at fracture intersections through experiments and simulations. By innovatively proposing a quantification process for fracture surface roughness, this study systematically researched the numerical simulation of proppant migration within rough fractures in unconventional reservoirs, including typical shales and sandstones, significantly saving manpower and resources. Some challenges emerge in quantitatively describing the flow mechanisms near complex structural surfaces. Studies on the effects of irregular rough fractures on fluid and particle flow, velocity distribution, and inlet-outlet pressure differentials are scarce, especially regarding the dynamic transport of particles in offshore unconventional reservoirs. Based on a literature review and the core scanning results from hydraulic fracturing of offshore unconventional reservoirs, a precise methodology for quantifying fracture roughness is proposed. Combined with the influence of roughness on the characteristics of the flow region, the paper is designed to reveal the transport law of particles coupled with the flow field in rough fractures of rocks in different lithologies, optimize the placement process and efficiently implement efficient filling of proppants in fractures. The coupled CFD-DEM technique is applied to investigate the particle retention location, stacking structure and coverage area in the fracture to reveal the deposition behavior of the uniform and hybrid particle sizes in the fracture and to provide a reference for the proppant size selection and injection order. Consequently, the study of the particle transport and retention behavior of the solid phase within rough fractures offers important guidance for the optimization of the fracturing process and parameters.
section_name: Introduction

section_content: When a particle-liquid mixture flows, the two phases interfere with each other to some extent. For continuous phases, finite element software (CFD-Fluent 19. 2) is applied to solve the hydrodynamic control equations and simulate hydrodynamic problems. As the first software used to monitor particle motion patterns via the discrete element method, discrete element software (EDEM 2018) enables the tracking of individual particle motion and the simulation of interparticle interactions. Our research primarily focuses on the numerical solution of the particle-fluid mixture flow problem between the walls of rough fractures. 

Figure 1a shows the process of CFD-DEM coupling and a schematic diagram of the solid-liquid motion during the coupling process (Figure 1b ). This reflects the CFD-DEM calculation process as a two-way coupling, where the pressure and momentum obtained are transferred from the fluid flow to the particle motion process by calculating the pressure and velocity fields of the flow channel at the same moment. Tension and Saffman lift forces are added to calculate the updated particle positions and velocities. The forces of the particle reaction flow field are calculated to update the fluid pressure and velocity fields. This iterative cycle is repeated. 

for fracture surface roughness, this study systematically researched the numerical simulation of proppant migration within rough fractures in unconventional reservoirs, including typical shales and sandstones, significantly saving manpower and resources. Some challenges emerge in quantitatively describing the flow mechanisms near complex structural surfaces. Studies on the effects of irregular rough fractures on fluid and particle flow, velocity distribution, and inlet-outlet pressure differentials are scarce, especially regarding the dynamic transport of particles in offshore unconventional reservoirs. Based on a literature review and the core scanning results from hydraulic fracturing of offshore unconventional reservoirs, a precise methodology for quantifying fracture roughness is proposed. Combined with the influence of roughness on the characteristics of the flow region, the paper is designed to reveal the transport law of particles coupled with the flow field in rough fractures of rocks in different lithologies, optimize the placement process and efficiently implement efficient filling of proppants in fractures. The coupled CFD-DEM technique is applied to investigate the particle retention location, stacking structure and coverage area in the fracture to reveal the deposition behavior of the uniform and hybrid particle sizes in the fracture and to provide a reference for the proppant size selection and injection order. Consequently, the study of the particle transport and retention behavior of the solid phase within rough fractures offers important guidance for the optimization of the fracturing process and parameters.
section_name: Mathematical Models

section_content: When a particle-liquid mixture flows, the two phases interfere with each other to some extent. For continuous phases, finite element software (CFD-Fluent 19. 2) is applied to solve the hydrodynamic control equations and simulate hydrodynamic problems. As the first software used to monitor particle motion patterns via the discrete element method, discrete element software (EDEM 2018) enables the tracking of individual particle motion and the simulation of interparticle interactions. Our research primarily focuses on the numerical solution of the particle-fluid mixture flow problem between the walls of rough fractures. 

Figure 1a shows the process of CFD-DEM coupling and a schematic diagram of the solid-liquid motion during the coupling process (Figure 1b ). This reflects the CFD-DEM calculation process as a two-way coupling, where the pressure and momentum obtained are transferred from the fluid flow to the particle motion process by calculating the pressure and velocity fields of the flow channel at the same moment. Tension and Saffman lift forces are added to calculate the updated particle positions and velocities. The forces of the particle reaction flow field are calculated to update the fluid pressure and velocity fields. This iterative cycle is repeated. The continuous phase of the incompressible fluid obeys the Navier-Stokes control equations. The equations of continuity of motion and conservation of momentum are as follows [27, 28] :

where ρ g is the continuous phase fluid density, kg/m 3 ; u g is the fluid velocity, m/s; p is the fluid pressure, Pa; τ g is the fluid viscous stress tensor within the flow region; g is the acceleration of gravity, g = 9. 8 m/s 2 ; F g-p is the force between the particle phase and the fluid phase, N; ε g is the volume fraction of the fluid; and t is time, s. 

Considering the collisions between particles and particles and between particles and walls, the governing equation of the discrete term can be calculated by combining Newton's second law as follows:

where I i is the rotational inertia of particle i, kg•m 2 ; m i is the mass of particle i, kg; v i is the translational velocity of particle i, m/s; w i is the rotational velocity of particle i, r/min; k i is the number of particles in contact with particle i; I i is the torque, kg•m; f p-g and f contact,ij are the fluid-solid interaction and contact force between particle and fluid, respectively, N. 

Based on previous experience with the treatment of pressure drop correlation equations, the particle-fluid traction model can be formulated in terms of the interphase momentum transfer coefficient and the slip velocity. The equations are as follows [29] :

where d p is the diameter of the particle, m; µ g is the viscosity of the fluid, Pa•s; u p is the velocity of the particle flow field, m/s; and C D is the coefficient of traction. The traction coefficient is calculated as follows:

The Reynolds number is calculated as follows:

The pressure gradient of the fluid is as follows:

where V p is the volume of the particle, m 3, and dp dx is the pressure gradient along the x direction.
section_name: Mathematical Models

section_name: Establishment of the Rough Fractures Model

section_content: To ensure the reliability of the computational model, physical fracture propagation experiments were conducted using large fresh outcrop core samples (500 mm × 500 mm × 500 mm) from the WC formation in the L block of the Pearl River Mouth Basin offshore, China. In unconventional reservoirs, sandstone-mudstone interbedding is prevalent. For sandstone, traditional tri-cluster and segmented perforation methods were used for hydraulic fracturing and well stimulation. The specific fracturing details are listed in Table 1. Combining the tight sandstone rock sample (#1) with a brittleness index (BI) of 0. 40 after fracturing, as depicted in Figure 2, it is evident that a vertically penetrating fracture oriented perpendicular to the minimum horizontal principal stress direction formed (length: 500 mm, width: approximately 100 mm). As illustrated in Figure 3a, the curve can be divided into three segments. The initial stage (0~t 1 ) is characterized by pressure accumulation. Upon activating the servo pump pressure system and gradually injecting fracturing fluid, the pumping pressure in the sample increased linearly through the preset perforations, reaching the rupture pressure of the sample at 18. 9 MPa after t 1, resulting in fracture initiation. The subsequent stage (t 1 ~t2 ) represents the main fracture formation phase, during which the creation of the primary fracture enlarges the internal space within the rock mass, leading to an instantaneous pressure drop. The final stage (t 2 ~t3 ) signifies the extension of the main fracture, displaying a smooth curve with a typical symmetrical bimodal expansion pattern, and the fracture plane exhibits distinct rough surface features. Fracture propagation tests were conducted to obtain fracture surfaces in conjunction with different sandstones and shales from unconventional reservoirs. Due to the strong heterogeneity of offshore reservoirs, incorporating the brittleness index offers a comprehensive assessment of elastic modulus and Poisson's ratio for characterizing reservoir properties. Figure 3b illustrates the post-fracturing morphology of deeper sample #3 (BI_0. 62) under identical stress boundary conditions. The results indicate that rock sample #1, with a lower brittleness index, has a simpler and smoother surface morphology. In contrast, rock sample #3, with a higher BI, exhibits a more complex morphology and rougher surface after fracturing, producing more fractures. Certainly, designing the same fracturing process for the different reservoirs does not produce the desired results during oil and gas production. Due to the strong heterogeneity of offshore reservoirs, incorporating the brittleness index offers a comprehensive assessment of elastic modulus and Poisson's ratio for characterizing reservoir properties. Figure 3b illustrates the post-fracturing morphology of deeper sample #3 (BI_0. 62) under identical stress boundary conditions. The results indicate that rock sample #1, with a lower brittleness index, has a simpler and smoother surface morphology. In contrast, rock sample #3, with a higher BI, exhibits a more complex morphology and rougher surface after fracturing, producing more fractures. Certainly, designing the same fracturing process for the different reservoirs does not produce the desired results during oil and gas production. 

Due to the strong heterogeneity of offshore reservoirs, incorporating the brittleness index offers a comprehensive assessment of elastic modulus and Poisson's ratio for characterizing reservoir properties. Figure 3b illustrates the post-fracturing morphology of deeper sample #3 (BI_0. 62) under identical stress boundary conditions. The results indicate that rock sample #1, with a lower brittleness index, has a simpler and smoother surface morphology. In contrast, rock sample #3, with a higher BI, exhibits a more complex morphology and rougher surface after fracturing, producing more fractures. Certainly, designing the same fracturing process for the different reservoirs does not produce the desired results during oil and gas production.
section_name: Physical Experiments of Fracture Propagation

section_content: Rough fractures have uneven surfaces, with walls composed of upper and lower rough sections. The joint roughness coefficient (JRC) represents fracture surface roughness, a parameter considered by scholars worldwide to analyze aperture, shear stiffness, and hydraulic behavior [30]. Fractal theory quantifies irregular shapes, aiding in the description of fracture surface concavity and aperture, and its correlation with JRC.
section_name: Numerical Modelling of Rough Fractures

section_content: Rapid lithological changes in offshore unconventional reservoirs, along with frequent interbedding of sandstone and mudstone, can shorten proppant transport distances. Based on the rock characteristics of Block L, shale and tight sandstone were selected for proppant transport analysis. Generally, tight sandstone, while a type of sandstone, has higher strength and density compared to conventional sandstone and shale.
section_name: Point Cloud Data Processing

section_content: Rapid lithological changes in offshore unconventional reservoirs, along with frequent interbedding of sandstone and mudstone, can shorten proppant transport distances. Based on the rock characteristics of Block L, shale and tight sandstone were selected for proppant transport analysis. Generally, tight sandstone, while a type of sandstone, has higher strength and density compared to conventional sandstone and shale. However, as our study focuses on the impact of fracture surfaces on proppant transport, these rock mechanical properties have minimal influence and are therefore not considered in the analysis. 

To accurately obtain the fracture surface roughness, ten rock core samples (diameter: 100 mm) from different depths of the vertical well in the target reservoirs were scanned postfracturing to capture surface roughness variations. By comparing small-sized samples of sandstone and shale after fracturing (Figure 4a ), significant differences in surface roughness between different lithologies in the oilfield can be observed. Specifically, shale exhibited more pronounced surface undulations per unit area and a greater degree of roughness. Consequently, these characteristics significantly affected the transport and placement of proppant within the fractures. The use of threshold segmentation and binarization on the height maps derived from point cloud data generated a two-dimensional image distinguished by black and white regions. The combination of height and texture feature parameters was employed to quantitatively characterize macroscopic roughness. Industrial analysis and visualization software is used for thresholding the processed gray images. The process enables a combination of two types of regions (target and background regions) of different grayscale to be selected as a more reasonable threshold to determine which region each pixel point in the image should belong to. 

image distinguished by black and white regions. The combination of height and texture feature parameters was employed to quantitatively characterize macroscopic roughness. Industrial analysis and visualization software is used for thresholding the processed gray images. The process enables a combination of two types of regions (target and background regions) of different grayscale to be selected as a more reasonable threshold to determine which region each pixel point in the image should belong to. Combined with the fracture surface characteristics of rough fractured rock obtained by 3D laser scanner, the original data are transformed into point cloud data of convex rough walls in the Cartesian coordinate plane by computing the height of the bulge. To closely approximate the actual fracture shape, the rough surface dimensions are transformed by a Gaussian autocorrelated function with Equation (10). Meanwhile, fractal interpolation was used to simulate the actual roughness of this rock fracture surface to reduce the amount of data. 

where h (, )

C u v is the autocorrelated function; Combined with the fracture surface characteristics of rough fractured rock obtained by 3D laser scanner, the original data are transformed into point cloud data of convex rough walls in the Cartesian coordinate plane by computing the height of the bulge. To closely approximate the actual fracture shape, the rough surface dimensions are transformed by a Gaussian autocorrelated function with Equation (10). Meanwhile, fractal interpolation was used to simulate the actual roughness of this rock fracture surface to reduce the amount of data. 

where C h (u, v) is the autocorrelated function; σ h 2 is the square of the overall roughness; u and v are the lags along the x and y directions, respectively; l c1 and l c2 are the correlation length coefficients along the x and y directions, respectively. 

The macroscopic undulation was quantitatively characterized using a combination of height feature parameters and texture feature parameters to obtain a histogram of the distribution of apertures shown in Figure 5. The figure shows a Gaussian distribution of apertures, with the maximum aperture near 5 mm and the mean close to 0. The validity of our algorithm is verified, which is consistent with the results of the preliminary research above. of our algorithm is verified, which is consistent with the results of the preliminary research above.
section_name: Point Cloud Data Processing

section_content: By studying the fractal theory model, Xie established an empirical formula between JRC and the fractal dimension D, as shown in Figure 6a [31]. Subsequently, the fractal approximation was applied to analyze 10 typical nodal contours to obtain their corresponding fractal dimension and JRC. The fractal approximation of JRC was used to estimate the surface JRC values using 10 contour line center JRC values as a control group. The maximum error was approximately 1. 4% and the minimum error was approximately 0. 34%, demonstrating the high accuracy of the validated values. Taking JRC_28 as an example, we employed the box-counting method with a fixed box size increment. These results demonstrate that a box width increment (boxwidth_incr) of 200 provides better fitting performance, yielding a fractal dimension (D) of 1. 17 based on the calcu-
section_name: Quantification of Fracture Roughness

section_content: By studying the fractal theory model, Xie established an empirical formula between JRC and the fractal dimension D, as shown in Figure 6a [31]. Subsequently, the fractal approximation was applied to analyze 10 typical nodal contours to obtain their corresponding fractal dimension and JRC. The fractal approximation of JRC was used to estimate the surface JRC values using 10 contour line center JRC values as a control group. The maximum error was approximately 1. 4% and the minimum error was approximately 0. 34%, demonstrating the high accuracy of the validated values. Taking JRC_28 as an example, we employed the box-counting method with a fixed box size increment. These results demonstrate that a box width increment (boxwidth_incr) of 200 provides better fitting performance, yielding a fractal dimension (D) of 1. 17 based on the calculated slope in Figure 6b.
section_name: Quantification of Fracture Roughness

section_content: By studying the fractal theory model, Xie established an empirical formula between JRC and the fractal dimension D, as shown in Figure 6a [31]. Subsequently, the fractal approximation was applied to analyze 10 typical nodal contours to obtain their corresponding fractal dimension and JRC. The fractal approximation of JRC was used to estimate the surface JRC values using 10 contour line center JRC values as a control group. The maximum error was approximately 1. 4% and the minimum error was approximately 0. 34%, demonstrating the high accuracy of the validated values. Taking JRC_28 as an example, we employed the box-counting method with a fixed box size increment. These results demonstrate that a box width increment (boxwidth_incr) of 200 provides better fitting performance, yielding a fractal dimension (D) of 1. 17 based on the calculated slope in Figure 6b. The distribution of JRC after fracturing exhibited distinct patterns for tight sandstones and shales in offshore unconventional reservoirs, with shale exhibiting greater roughness, averaging 48, while the tight sandstone showed lower values, with JRC at 37, as shown in Figure 7. The distribution of JRC after fracturing exhibited distinct patterns for tight sandstones and shales in offshore unconventional reservoirs, with shale exhibiting greater roughness, averaging 48, while the tight sandstone showed lower values, with JRC at 37, as shown in Figure 7.
section_name: Quantification of Fracture Roughness

section_content: To investigate the arrangement pattern of proppant within hydraulic fractures of different rock samples in offshore unconventional reservoirs, an autocorrelated Gaussian surface point cloud data and 3D surfaces under different roughness conditions were established based on the fractal dimension calculation process shown in Figure 4. The fractal dimension of the surface of the five binarized images was calculated by the box dimension method, yielding five joint roughness (JRC) values of 52, 48, 28, 20 and 0, as shown in Figure 8. The figure shows that surfaces with larger JRC values exhibit greater aperture fluctuations per unit area. Different joint roughness coefficients are defined for offshore shale (JRC_48,52) and tight sandstone (JRC_20,28).
section_name: Establishment of Flow Channels with Different Roughnesses

section_content: To investigate the arrangement pattern of proppant within hydraulic fractures of different rock samples in offshore unconventional reservoirs, an autocorrelated Gaussian surface point cloud data and 3D surfaces under different roughness conditions were established based on the fractal dimension calculation process shown in Figure 4. The fractal dimension of the surface of the five binarized images was calculated by the box dimension method, yielding five joint roughness (JRC) values of 52, 48, 28, 20 and 0, as shown in Figure 8. The figure shows that surfaces with larger JRC values exhibit greater aperture fluctuations per unit area. Different joint roughness coefficients are defined for offshore shale (JRC_48,52) and tight sandstone (JRC_20,28). 

To investigate the arrangement pattern of proppant within hydraulic fractures of different rock samples in offshore unconventional reservoirs, an autocorrelated Gaussian surface point cloud data and 3D surfaces under different roughness conditions were established based on the fractal dimension calculation process shown in Figure 4. The fractal dimension of the surface of the five binarized images was calculated by the box dimension method, yielding five joint roughness (JRC) values of 52, 48, 28, 20 and 0, as shown in Figure 8. The figure shows that surfaces with larger JRC values exhibit greater aperture fluctuations per unit area. Different joint roughness coefficients are defined for offshore shale (JRC_48,52) and tight sandstone (JRC_20,28). As previously mentioned, by combining with the point cloud data of rough fractures at different coordinates established by scanning, we established the rough surfaces of the rock fractures via interpolation. Meanwhile, the dimensions of the surfaces were designed to closely match the actual fracture shape. Linear interpolation is commonly used for data smoothing. The y-value of an unknown data point can be estimated using the slope m, m As previously mentioned, by combining with the point cloud data of rough fractures at different coordinates established by scanning, we established the rough surfaces of the rock fractures via interpolation. Meanwhile, the dimensions of the surfaces were designed to closely match the actual fracture shape. Linear interpolation is commonly used for data smoothing. The y-value of an unknown data point can be estimated using the slope m, m = (y 2y 1 )/(x 2x 1 ), with the formula y = y 1 + m × (xx 1 ), given two known data points (x 1, y 1 ) and (x 2, y 2 ). 

Taking the rough surface of the JRC_28 as an example, a flow channel size of 500 mm × 100 mm × 1 mm rough fracture is modeled based on the actual size of the fractures produced. The polygon curve is first created from the point cloud data. The interpolated 3D function is then obtained using linear cell interpolation. Thereafter, by combining this function to create a parametric surface to obtain the autocorrelated rough surface of the diagram, the corresponding rough fractured flow channel model is obtained by Boolean operations on this surface, as shown in Figure 9. Taking the rough surface of the JRC_28 as an example, a flow channel size of 500 mm × 100 mm × 1 mm rough fracture is modeled based on the actual size of the fractures produced. The polygon curve is first created from the point cloud data. The interpolated 3D function is then obtained using linear cell interpolation. Thereafter, by combining this function to create a parametric surface to obtain the autocorrelated rough surface of the diagram, the corresponding rough fractured flow channel model is obtained by Boolean operations on this surface, as shown in Figure 9. As shown in Figure 10, several models of the flow channel with different joint roughness coefficients in unconventional reservoirs were established, based on the above method. The models for these five different roughnesses are of identical length and width. As shown in Figure 10, several models of the flow channel with different joint roughness coefficients in unconventional reservoirs were established, based on the above method. The models for these five different roughnesses are of identical length and width.
section_name: Establishment of Flow Channels with Different Roughnesses

section_name: Validation of the Model

section_content: On the premise that the boundary conditions of the numerical simulation are consistent with the parameters of fluid in the field, the plane injection method is used in the study of mixed phases. The simulation defines the left surface of the near-well section with the velocity inlet; the outlet is set on the whole right side with the pressure outlet; all but
section_name: Establishment of Boundary Conditions

section_content: On the premise that the boundary conditions of the numerical simulation are consistent with the parameters of fluid in the field, the plane injection method is used in the study of mixed phases. The simulation defines the left surface of the near-well section with the velocity inlet; the outlet is set on the whole right side with the pressure outlet; all but these two surfaces are set as walls; and the negative direction of the y-axis is set as the direction of gravity. In the DEM, a particle factory of uniform size was set at the inlet of the flow channel. The particle parameters include the recovery coefficient and static friction coefficient for different sizes of proppant particles with homogeneous particle size, which was set to 0. 5, the rolling friction coefficient was set to 0. 4 and the generation rate (generation rate) was set to 0. 02 kg/s. Since the transport of proppant in fractures is affected by fluid viscosity, and the size of the drag coefficient acting on particles is related to the shape of particles, the particles in this paper were set as the same circular particles. The boundary conditions are shown in Figure 11. 

these two surfaces are set as walls; and the negative direction of the y-axis is set as the direction of gravity. In the DEM, a particle factory of uniform size was set at the inlet of the flow channel. The particle parameters include the recovery coefficient and static friction coefficient for different sizes of proppant particles with homogeneous particle size, which was set to 0. 5, the rolling friction coefficient was set to 0. 4 and the generation rate (generation rate) was set to 0. 02 kg/s. Since the transport of proppant in fractures is affected by fluid viscosity, and the size of the drag coefficient acting on particles is related to the shape of particles, the particles in this paper were set as the same circular particles. The boundary conditions are shown in Figure 11. Typically, the anisotropy of the surface of a rough model is particularly complex, which increases the difficulty of meshing. Hence, to better ensure computational accuracy, the structured meshing method is adopted in this study. Ensuring iterative convergence in the coupled calculation (CFD-DEM) requires that the compiled coupling interface meets the condition where the volume of the mesh elements in the DEM is greater than 2. 5 times the volume of the particles. Based on the above division method, research set the minimum cell size length of the mesh grid of JRC_28 to approximately 2. 3 mm. To optimize the results of the selection, the number of grid was chosen to be 10,500 by monitoring the turbulent kinetic energy and absolute pressure at different positions along the x-axis at the same height with several different numbers of grids to reduce the simulation time while ensuring the accuracy of the calculation. The analysis in Figure 12 shows that the results guarantee the accuracy of the calculation while satisfying a grid quality of over Typically, the anisotropy of the surface of a rough model is particularly complex, which increases the difficulty of meshing. Hence, to better ensure computational accuracy, the structured meshing method is adopted in this study. Ensuring iterative convergence in the coupled calculation (CFD-DEM) requires that the compiled coupling interface meets the condition where the volume of the mesh elements in the DEM is greater than 2. 5 times the volume of the particles. Based on the above division method, research set the minimum cell size length of the mesh grid of JRC_28 to approximately 2. 3 mm. To optimize the results of the selection, the number of grid was chosen to be 10,500 by monitoring the turbulent kinetic energy and absolute pressure at different positions along the x-axis at the same height with several different numbers of grids to reduce the simulation time while ensuring the accuracy of the calculation. The analysis in Figure 12 shows that the results guarantee the accuracy of the calculation while satisfying a grid quality of over 95%. 

which increases the difficulty of meshing. Hence, to better ensure computational accuracy, the structured meshing method is adopted in this study. Ensuring iterative convergence in the coupled calculation (CFD-DEM) requires that the compiled coupling interface meets the condition where the volume of the mesh elements in the DEM is greater than 2. 5 times the volume of the particles. Based on the above division method, research set the minimum cell size length of the mesh grid of JRC_28 to approximately 2. 3 mm. To optimize the results of the selection, the number of grid was chosen to be 10,500 by monitoring the turbulent kinetic energy and absolute pressure at different positions along the x-axis at the same height with several different numbers of grids to reduce the simulation time while ensuring the accuracy of the calculation. The analysis in Figure 12 shows that the results guarantee the accuracy of the calculation while satisfying a grid quality of over 95%. To reveal the transport and placement of proppant in the flow channel under different roughnesses and particle sizes, the test parameters set in this paper are combined with the field production conditions. The viscosity was set at 0. 005 Pa•s for general fracturing fluids of slickwater. Simultaneously, the initial particle diameter was set at 0. 6 mm and the JRC_28. The properties of particle and fracturing fluid are shown in Table 2 To reveal the transport and placement of proppant in the flow channel under different roughnesses and particle sizes, the test parameters set in this paper are combined with the field production conditions. The viscosity was set at 0. 005 Pa•s for general fracturing fluids of slickwater. Simultaneously, the initial particle diameter was set at 0. 6 mm and the JRC_28. The properties of particle and fracturing fluid are shown in Table 2.
section_name: Validation of the Model 4.1. Establishment of Boundary Conditions

section_content: To verify the accuracy of the model, the study adopts the same experimental parameters as in the literature. Comparing a smooth plate model of the same size (380 × 190 × 2 mm), the results validate our model by monitoring the average settling velocity of particles within 4 s, compared to experimental data from the literature. The average error between the model results and the experimental data in Figure 13 is 8. 48%, which indicates a strong match between our model and the experimental results. The error is primarily due to the joint fracture established by Bo Zhang et al., which resulted in a small fraction of the particles entering the branch fracture, leading to a low particle settling velocity [32]. Altering the particle injection method from a fixed entry point to random injection using external programming can enhance alignment with real production and reduce errors to some extent. In addition, the study also tracked the particle deposition velocity at different times with high roughness (JRC_48). By comparison, the results show that the roughness of the surfaces hampers the movement of particles, which leads to an increase in the irregular movement of particles. Moreover, the particles are in the state of suspension for a longer period of time, a trend that facilitates the transport of particles during the hydraulic fracturing process. 

due to the joint fracture established by Bo Zhang et al., which resulted in a small fraction of the particles entering the branch fracture, leading to a low particle settling velocity [32]. Altering the particle injection method from a fixed entry point to random injection using external programming can enhance alignment with real production and reduce errors to some extent. In addition, the study also tracked the particle deposition velocity at different times with high roughness (JRC_48). By comparison, the results show that the roughness of the surfaces hampers the movement of particles, which leads to an increase in the irregular movement of particles. Moreover, the particles are in the state of suspension for a longer period of time, a trend that facilitates the transport of particles during the hydraulic fracturing process.
section_name: Verification of the Numerical Simulation Method

section_name: Transport Analysis of Proppants in Rough Fractures

section_content: As previously mentioned, the study employs two transport equations from the existing K-epsilon model to describe the turbulent characteristics within rough fractures. Turbulent kinetic energy (K) governs energy in turbulent motion and notably affects the energy dissipation of multi-phase particle motion in turbulent flow, a crucial consideration in numerical simulations [33]. Given the viscous nature of fracturing fluid, wall shear stress significantly shapes flow behavior at rough boundaries and controls wall friction resistance, exerting substantial influence on frictional drag reduction in two-phase flow. 

Figure 14 depicts the wall shear force and turbulent kinetic energy distribution curves along the x-axis for the same cross-section of JRC_28 (rough tight sandstone) and smooth fractures under the same operating conditions. Furthermore, the figure illustrates an explicit resemblance between the turbulent kinetic energy of the fluid against the rough wall and the wall shear intensity distribution. Both turbulent kinetic energy and wall shear exhibit large fluctuations compared to smooth walls. The turbulent kinetic energy near rough surfaces (7. 08 m 2 /s 2 ) is 354 times greater than that near smooth wall surfaces (0. 02 m 2 /s 2 ). Consequently, the magnitude of fracture roughness significantly impacts proppant transport within fractures. The process of proppant transport within rough fractures is also subject to greater wall and fluid resistance, as well as greater particle energy dissipation, than that within smooth walls. 

The volume fraction of the proppant per unit grid volume is calculated using the discrete element method. Figure 15a shows the contours of the volume fraction distribution for five discrete phases with different roughnesses, which clearly reflect the phenomenon when the fracture wall roughness is large (JRC_48,52) and the amount of proppant deposited at the bottom is small. Compared to the premise of a dense proppant distribution at the bottom of a tight sandstone fracture (red dashed line), greater roughness corresponds to a greater probability of regions with proppant volume fraction below 0. 55. This discovery implies fewer particles per unit area, indicating that under flow conditions, a relatively greater amount of proppants can be trapped near the inlet. Moreover, varying shear-stress curves on the upper surfaces of different fracture planes P90, as shown in Figure 15b, indicate that rougher fracture planes significantly hinder both the proppants and fluids. The trends clearly reflect the distributions of deposited and suspended proppant within the smooth fractures are homogeneous. In conclusion, with increasing roughness, the proppant deposition at the bottom becomes more chaotic and irregular. At a roughness of JRC_52, proppants are widely distributed within the fracture, becoming trapped in larger curvatures. The proppant bed exhibits a high near-wellbore and low far-wellbore effect, severely impeding the backflow of proppants and causing sand loss. 

Figure 14 depicts the wall shear force and turbulent kinetic energy distribution curves along the x-axis for the same cross-section of JRC_28 (rough tight sandstone) and smooth fractures under the same operating conditions. Furthermore, the figure illustrates an explicit resemblance between the turbulent kinetic energy of the fluid against the rough wall and the wall shear intensity distribution. Both turbulent kinetic energy and wall shear exhibit large fluctuations compared to smooth walls. The turbulent kinetic energy near rough surfaces (7. 08 m 2 /s 2 ) is 354 times greater than that near smooth wall surfaces (0. 02 m 2 /s 2 ). Consequently, the magnitude of fracture roughness significantly impacts proppant transport within fractures. The process of proppant transport within rough fractures is also subject to greater wall and fluid resistance, as well as greater particle energy dissipation, than that within smooth walls. The volume fraction of the proppant per unit grid volume is calculated using the discrete element method. Figure 15a shows the contours of the volume fraction distribution for five discrete phases with different roughnesses, which clearly reflect the phenomenon when the fracture wall roughness is large (JRC_48,52) and the amount of proppant deposited at the bottom is small. Compared to the premise of a dense proppant distribution at the bottom of a tight sandstone fracture (red dashed line), greater roughness corresponds to a greater probability of regions with proppant volume fraction below 0. 55. This discovery implies fewer particles per unit area, indicating that under flow conditions, a relatively greater amount of proppants can be trapped near the inlet. Moreover, varying shear-stress curves on the upper surfaces of different fracture planes P90, as shown in Figure 15b, indicate that rougher fracture planes significantly hinder both the proppants and fluids. The trends clearly reflect the distributions of deposited and suspended proppant within the smooth fractures are homogeneous. In conclusion, with increasing roughness, the proppant deposition at the bottom becomes more chaotic and irregular. At a roughness of JRC_52, proppants are widely distributed within the fracture, becoming trapped in larger curvatures. The proppant bed exhibits a high near-wellbore and low farwellbore effect, severely impeding the backflow of proppants and causing sand loss.
section_name: Analysis of the Flow Field

section_content: The process of proppant deposition is influenced by different factors, including the proppant properties, physical characteristics of fractures, and fracturing fluid parameters, etc. Consequently, as the focus of this paper is on the influence of wall characteristics on proppant transport, three representative trials were selected for analysis. 

Figure 16 illustrates the characteristics of proppant transport behavior within rough tight sandstone fractures at different times. The background represents the contour of the height for fracture aperture distribution along the z-axis. Due to their higher density compared to fracturing fluid, proppant particles tend to settle at the fracture bottom, forming a proppant bed once injected from the wellbore. A combination of monitoring and tracking the proppants reveals that the injected early particles will remain largely immobile at the bottom. The proppant bed becomes higher as more proppant particles settle until a balance is reached between proppant settlement and friction-induced scouring. Ultimately, three zones are formed: the deposition zone, the sliding of particles atop the deposition area (sliding zone), and suspended movement within the fracture (suspension zone). As a result, the newly injected particles will spread farther, creating a longer proppant bed. However, in this flow scenario, proppants tend to halt and accumulate near the wellbore within rough fractures, hindering the movement of subsequently injected proppant toward the distant end.
section_name: The Process of Proppant Transport

section_content: The process of proppant deposition is influenced by different factors, including the proppant properties, physical characteristics of fractures, and fracturing fluid parameters, etc. Consequently, as the focus of this paper is on the influence of wall characteristics on proppant transport, three representative trials were selected for analysis. 

Figure 16 illustrates the characteristics of proppant transport behavior within rough tight sandstone fractures at different times. The background represents the contour of the height for fracture aperture distribution along the z-axis. Due to their higher density compared to fracturing fluid, proppant particles tend to settle at the fracture bottom, forming a proppant bed once injected from the wellbore. A combination of monitoring and tracking the proppants reveals that the injected early particles will remain largely immobile at the bottom. The proppant bed becomes higher as more proppant particles settle until a balance is reached between proppant settlement and friction-induced scouring. Ultimately, three zones are formed: the deposition zone, the sliding of particles atop the deposition area (sliding zone), and suspended movement within the fracture (suspension zone). As a result, the newly injected particles will spread farther, creating a longer proppant bed. However, in this flow scenario, proppants tend to halt and accumulate near the wellbore within rough fractures, hindering the movement of subsequently injected proppant toward the distant end. 

ing the proppants reveals that the injected early particles will remain largely immobile at the bottom. The proppant bed becomes higher as more proppant particles settle until a balance is reached between proppant settlement and friction-induced scouring. Ultimately, three zones are formed: the deposition zone, the sliding of particles atop the deposition area (sliding zone), and suspended movement within the fracture (suspension zone). As a result, the newly injected particles will spread farther, creating a longer proppant bed. However, in this flow scenario, proppants tend to halt and accumulate near the wellbore within rough fractures, hindering the movement of subsequently injected proppant toward the distant end. Figure 17 reveals the transport of the particle proppant for large values of shale (JRC_52) within different moments. Based on the contour at three different times, it can be detected that the effect of plugging occurs where the particles pass from the start of the timing. The proppant is observed by magnification to be clustered in the grooves with greater surface curvature (larger absolute value of z coordinate). Surprisingly, a significant portion of these particles accumulate in the lower half of the grooves, resulting in plugging that obstructs the downstream movement of subsequent proppants, thus confirming the previous findings. To observe the trajectory of the particles, the substitution of color was implemented for the particles. The black dashed line at t = 10. 5 shows that as the proppant bed rises, the movement of the particles gradually creates a high-speed channel due to the rough surface of the fractures. The proppant will move farther along the channel and will remain in this state of motion for an extended period of time. Compared with Figure 15, higher roughness elevates the wall shear force within the fractures, consequently causing greater disturbance and velocity of proppant transport. With the fluid entering the fractures, the proppant settles rapidly at the entrance to the fracture and forms a proppant dike, which increases in height as the fracturing fluid is injected and gradually reaches a relatively equilibrium height. Particularly, the analysis revealed that the maximum value of particle velocity at t = 15. 2 s in conventional proppant transport occurs near the tip of the proppant bed build-up above the same location as the smooth wall. Upon reaching the equilibrium height, the proppant bed reduces the flow area of the fracturing fluid, resulting in heightened velocity, thereby enhancing the proppant-carrying capacity towards the far side of the fracture. Combining the results of particle volume fractions at different roughness models, a clear conclusion emerges: the particle arrangement within the suspension and deposition regions of rough fractures is uneven. Due to the pronounced curvature, gaps form between particles in the deposition region. 

the equilibrium height, the proppant bed reduces the flow area of the fracturing fluid, resulting in heightened velocity, thereby enhancing the proppant-carrying capacity towards the far side of the fracture. Combining the results of particle volume fractions at different roughness models, a clear conclusion emerges: the particle arrangement within the suspension and deposition regions of rough fractures is uneven. Due to the pronounced curvature, gaps form between particles in the deposition region.
section_name: The Process of Proppant Transport

section_content: The proppant deposition height and coverage area in the fracture provide lateral indications of the proppant clogging and fracture support. In addition, the profile of the suspension zone can also be applied to characterize the distance the proppant travels.
section_name: Deposition in Fractures with Different Roughnesses

section_content: The proppant deposition height and coverage area in the fracture provide lateral indications of the proppant clogging and fracture support. In addition, the profile of the suspension zone can also be applied to characterize the distance the proppant travels. 

The relationship between the deposition height and the transport distance was established in our research to characterize the coordinates of the proppant movement position for different JRC premises, as shown in Figure 18. Corresponding to the results obtained in the previous section, the height of particle deposition at the inlet will be higher for particles with larger roughness. However, the farther away from the inlet, the more significantly the height of the deposited particles decreases. Considering the solid lines in the figure illustrating the farthest suspension distance of proppant at various rough fractures, we can deduce that rougher fractures lead to proppant traveling deeper into the flow region and with broader distribution. The conclusion stems from the fact that the properties of the concave and convex rough surfaces affect the flow behavior of the fracturing fluid, resulting in vortices in the flow near the wall. Consequently, greater wall roughness results in stronger vortex flow within the fracture, which is consistent with the findings of Guo et al. [34]. As previously analyzed, when the value of JRC is large, the particles have the effect of plugging. Later, as the proppant bed becomes higher, these plugging proppants are incrementally eliminated by the larger y coordinate. 

By monitoring the number of particles through the inlet to the flow region, the total of which is recorded as C. The derived particle velocity data are then divided into intervals and the number of suspended particles in the flow region is counted as C 0. Comparing the value of C 0 /C at different moments in the two types of fractures, JRC_52 and smooth, the values reflect the probability of the suspended solid phase being retained as it passes through the fractures. The greater the total number of proppants that exits out of a rough fracture compared to a smooth fracture, the more pronounced the proppant suspension will be and the better the effect of transport will be. Accordingly, the smooth model is less capable of transporting and carrying sand, and the effective length of proppant in the resulting fracture is relatively short. In Figure 19a, compared to smooth fractures, rough fractures exhibit significant liquid turbulence, which increases the total efflux of proppant and enhances the suspension of proppant. This results in improved proppant migration efficiency. Figure 19b shows the statistical analysis of proppant collision frequencies on two surfaces with different roughnesses, which supports this observation. Some particles can remain within fractures for extended periods, increasing the likelihood of collisions with the walls, particularly on rough surfaces. Due to the limited influence of turbulence, smooth fractures exhibit a weaker proppant-carrying capacity, leading to a relatively shorter effective propped fracture length. 

for different JRC premises, as shown in Figure 18. Corresponding to the results obtained in the previous section, the height of particle deposition at the inlet will be higher for particles with larger roughness. However, the farther away from the inlet, the more significantly the height of the deposited particles decreases. Considering the solid lines in the figure illustrating the farthest suspension distance of proppant at various rough fractures, we can deduce that rougher fractures lead to proppant traveling deeper into the flow region and with broader distribution. The conclusion stems from the fact that the properties of the concave and convex rough surfaces affect the flow behavior of the fracturing fluid, resulting in vortices in the flow near the wall. Consequently, greater wall roughness results in stronger vortex flow within the fracture, which is consistent with the findings of Guo et al. [34]. As previously analyzed, when the value of JRC is large, the particles have the effect of plugging. Later, as the proppant bed becomes higher, these plugging proppants are incrementally eliminated by the larger y coordinate. By monitoring the number of particles through the inlet to the flow region, the total of which is recorded as C. The derived particle velocity data are then divided into intervals and the number of suspended particles in the flow region is counted as C0. Comparing the value of C0/C at different moments in the two types of fractures, JRC_52 and smooth, the values reflect the probability of the suspended solid phase being retained as it passes through the fractures. The greater the total number of proppants that exits out of a rough fracture compared to a smooth fracture, the more pronounced the proppant suspension will be and the better the effect of transport will be. Accordingly, the smooth model is less capable of transporting and carrying sand, and the effective length of proppant in the resulting fracture is relatively short. In Figure 19a, compared to smooth fractures, rough fractures exhibit significant liquid turbulence, which increases the total efflux of proppant and enhances the suspension of proppant. This results in improved proppant migration efficiency. Figure 19b shows the statistical analysis of proppant collision frequencies on two surfaces with different roughnesses, which supports this observation. Some particles can remain within fractures for extended periods, increasing the likelihood of collisions with the walls, particularly on rough surfaces. Due to the limited influence of turbulence, smooth fractures exhibit a weaker proppant-carrying capacity, leading to a relatively shorter effective propped fracture length.
section_name: Deposition in Fractures with Different Roughnesses

section_content: Over the years, most scholars have found that the rough surfaces of fractures have important influences on proppant particle transport. The d/w (ratio of particle size to fracture width) affects the ability of the proppant to pass through the fracture. Yet, there is no specific research on d/w in relation to the degree of proppant transport on rough surfaces and the position of final placement. In contrast, the location and degree of proppant filling within fractures in practical operations significantly influence the final fracturing effect [35]. 

For the defined fracture width with a median JRC value of 28 and keeping other boundary conditions consistent, a computational time of 10 s is configured due to the larger particle size expediting deposition. Varying the d/w to 0. 3, 0. 4, 0. 6, 0. 7, and 0. 8, Figure 20 reveals the particle volume fraction distribution. Notably, when the particle size is small (d/w_0. 3), the figure obviously reflects that some of the proppant settles during transport, forming a thin, multi-segmented bed of raised proppant. Accordingly, in the enlarged section of the contours, the velocity of the fluid was maximal at the top of the bed within the proppant bed, which was primarily attributed to the narrowing of the fluid channels caused by the elevated proppant bed. The proppant that had settled remained mobile under the disturbance of sand-carrying fluid. Meanwhile, with a further injection of fracturing fluid, the proppant concentration gradient at the leading edge of the deposited zone soon reached a steady state, but most of the proppant remained suspended in the flow region and entered the branch fractures or moved deeper into the main fractures. This phenomenon indicates that the small proppant can greatly improve the carrying efficiency of the proppant system.
section_name: Deposition Behavior of Different Sizes for Uniform Proppants

section_content: Over the years, most scholars have found that the rough surfaces of fractures have important influences on proppant particle transport. The d/w (ratio of particle size to fracture width) affects the ability of the proppant to pass through the fracture. Yet, there is no specific research on d/w in relation to the degree of proppant transport on rough surfaces and the position of final placement. In contrast, the location and degree of proppant filling within fractures in practical operations significantly influence the final fracturing effect [35]. 

For the defined fracture width with a median JRC value of 28 and keeping other boundary conditions consistent, a computational time of 10 s is configured due to the larger particle size expediting deposition. Varying the d/w to 0. 3, 0. 4, 0. 6, 0. 7, and 0. 8, Figure 20 reveals the particle volume fraction distribution. Notably, when the particle size is small (d/w_0. 3), the figure obviously reflects that some of the proppant settles during transport, forming a thin, multi-segmented bed of raised proppant. Accordingly, in the enlarged section of the contours, the velocity of the fluid was maximal at the top of the bed within the proppant bed, which was primarily attributed to the narrowing of the fluid channels caused by the elevated proppant bed. The proppant that had settled remained mobile under the disturbance of sand-carrying fluid. Meanwhile, with a further injection of fracturing fluid, the proppant concentration gradient at the leading edge of the deposited zone soon reached a steady state, but most of the proppant remained suspended in the flow region and entered the branch fractures or moved deeper into the main fractures. This phenomenon indicates that the small proppant can greatly improve the carrying efficiency of the proppant system. The proportion of proppant in the suspension zone is relatively large when d/w is 0. 3. The proppant slowly forms some deposition at the entrance and improves the filling effect near the well. Additionally, for larger proppant sizes (d/w_0. 8), the flow of particles form a stable channel between the fractures. In comparison to the previous study on roughness, the channel accelerates the movement of particles to depth, and the results are in agreement with the conclusions of the previous section. Excessively large proppant sizes, however, lead to sand clogging at the entrance of the fracture channel for a short period. By visualizing the different particle-size distributions within the flow channel in the coupled simulation DEM, we have established the particle size distribution schematic shown in Figure 21. The figure illustrates that the proppant is staggered in the fracture and not neatly distributed, which ensures the stability of the proppant bed. The proportion of proppant in the suspension zone is relatively large when d/w is 0. 3. The proppant slowly forms some deposition at the entrance and improves the filling effect near the well. Additionally, for larger proppant sizes (d/w_0. 8), the flow of particles form a stable channel between the fractures. In comparison to the previous study on roughness, the channel accelerates the movement of particles to depth, and the results are in agreement with the conclusions of the previous section. Excessively large proppant sizes, however, lead to sand clogging at the entrance of the fracture channel for a short period. By visualizing the different particle-size distributions within the flow channel in the coupled simulation DEM, we have established the particle size distribution schematic shown in Figure 21. The figure illustrates that the proppant is staggered in the fracture and not neatly distributed, which ensures the stability of the proppant bed. The proportion of proppant in the suspension zone is relatively large when d/w is 0. 3. The proppant slowly forms some deposition at the entrance and improves the filling effect near the well. Additionally, for larger proppant sizes (d/w_0. 8), the flow of particles form a stable channel between the fractures. In comparison to the previous study on roughness, the channel accelerates the movement of particles to depth, and the results are in agreement with the conclusions of the previous section. Excessively large proppant sizes, however, lead to sand clogging at the entrance of the fracture channel for a short period. By visualizing the different particle-size distributions within the flow channel in the coupled simulation DEM, we have established the particle size distribution schematic shown in Figure 21. The figure illustrates that the proppant is staggered in the fracture and not neatly distributed, which ensures the stability of the proppant bed.
section_name: Deposition Behavior of Different Sizes for Uniform Proppants

section_name: Sequence of Proppant Injection

section_content: Research on the transport and deposition pattern of unified particles in rough fractures in the last section has established a rule that large particles can promptly form supports near the well and small particles are exceptionally transportable. Nonetheless, mixed particle size is also a common way of adding sand during production. Thus, under the consistent boundary conditions, the paper examines the deposition proppant for mixed particle sizes of 0. 4, 0. 6 and 0. 8 mm. The injection model is simplified with a consistent displacement to study the effects of different proppant injection sequences: small to large (S:M:L), large to small (L:M:S), and mixed injection (nix), on proppant transport and placement. 

Figure 22a presents the velocity distribution of the three types of particles between the fractures within 15 s under mixed particle-size conditions. Combined with the deposition curve in Figure 17, this vividly illustrates that under mixed particle-size injection conditions, the maximum proppant deposition height (52. 1 mm) and coverage area are notably greater than those of uniform particles (41. 5 mm). For dynamically deposited proppant beds, the same pattern of greater particle velocities at the top than elsewhere emerges at increasing heights. 

Research on the transport and deposition pattern of unified particles in rough fractures in the last section has established a rule that large particles can promptly form supports near the well and small particles are exceptionally transportable. Nonetheless, mixed particle size is also a common way of adding sand during production. Thus, under the consistent boundary conditions, the paper examines the deposition proppant for mixed particle sizes of 0. 4, 0. 6 and 0. 8 mm. The injection model is simplified with a consistent displacement to study the effects of different proppant injection sequences: small to large (S:M:L), large to small (L:M:S), and mixed injection (nix), on proppant transport and placement. 

Figure 22a presents the velocity distribution of the three types of particles between the fractures within 15 s under mixed particle-size conditions. Combined with the deposition curve in Figure 17, this vividly illustrates that under mixed particle-size injection conditions, the maximum proppant deposition height (52. 1 mm) and coverage area are notably greater than those of uniform particles (41. 5 mm). For dynamically deposited proppant beds, the same pattern of greater particle velocities at the top than elsewhere emerges at increasing heights. A comparative analysis between the S:M:L and L:M:S methods reveals that the initial injection of larger proppants results in an excessively elevated sand accumulation zone near the wellbore. Consequently, the transport distance of smaller proppants is significantly reduced. Specifically, with 0. 6 mm proppant, the transport distance decreases by 9%. Furthermore, the maximum height of the proppant bed near the wellbore increases by 17%. Conversely, adopting the S:M:L injection sequence ensures an overall enhancement in proppant transport efficiency, promoting deeper proppant transport while improving support closer to the wellbore during later stages. Additionally, when comparing the uniform particle size injection method, the utilization of a mixed particle-size composition Consequently, the transport distance of smaller proppants is significantly reduced. Specifically, with 0. 6 mm proppant, the transport distance decreases by 9%. Furthermore, the maximum height of the proppant bed near the wellbore increases by 17%. Conversely, adopting the S:M:L injection sequence ensures an overall enhancement in proppant transport efficiency, promoting deeper proppant transport while improving support closer to the wellbore during later stages. Additionally, when comparing the uniform particle size injection method, the utilization of a mixed particle-size composition intensifies the risk of localized fracture plugging. Detailed examination reveals that both plugging and gaps primarily stem from the larger particle sizes accumulating at high-curvature fracture regions (solid black line in the enlarged figure), progressively leading to blockage formation. Thus, maintaining consistent proppant particle sizes throughout production process is advised to mitigate the occurrence of sand plugging. Consequently, maintaining small proppants (d/w < 0. 6) throughout the production process is recommended for offshore shale reservoirs to minimize the incidence of plugging.
section_name: Sequence of Proppant Injection

section_content: Existing studies typically focus on the transport process of proppants in smooth fractures or do not systematically quantify fracture surface roughness, leading to inaccurate proppant transport results. The proppant transport patterns in unconventional reservoir rocks have been systematically introduced. Additionally, based on the roughness distribution of fracture surfaces across different rock samples, a comprehensive method for the quantitative characterization of fracture surface roughness was systematically established. This method can be applied to the quantitative calculation of any rough rock surface. By reconstructing fracture surfaces or generating point cloud data with the same roughness in MATLAB, Boolean operations were used to create 3D fracture models for simulating proppant transport. However, the simulation model developed in this study is specific to the reservoirs in block L, and it has certain limitations. Future research, with more hydraulic post-fracture surface data, whether from unconventional or conventional reservoirs onshore, could expand our findings. This would enable the recommendation of optimal proppant injection methods and other engineering parameters tailored to specific reservoirs. Further exploration of this challenging case will be conducted in future work.
section_name: Discussion

section_content: The reconstruction of post-fracturing fracture surface roughness in offshore unconventional reservoirs with different lithologies was carried out, in combination with true triaxial tests. The roughness was identified by the box dimension method, which is quantified by fractal dimension. Bi-directional coupling between finite element and discrete element software (Fluent-DEM, Ansys19. 2 and EDEM2018) was used to evaluate the transport and placement of proppant under different JRC and proppant size conditions for fracture channels after Boolean operation.
section_name: Conclusions

section_content: The average roughness of shale (JRC_48) is greater compared to tight sandstone (JRC_37) after hydraulic fracturing. Compared with smooth fractures, high JRC values can cause blockages between fractures, increasing transport irregularity. This results in the rapid formation of high-velocity channels near the inlet, fundamentally altering the flow dynamics of oil or natural gas within the fractures. 2. 

In conventional proppant transport processes, fracturing fluid entering fractures causes rapid settling of proppant near the inlet, forming an equilibrating bed. As equilibrium is reached, the flow area of fluid decreases, enhancing its sand-carrying capacity and accelerating proppant transport into deeper fractures.
section_name: 1.

section_content: When d/w is very small (d/w < 0. 4), the formation of a stable proppant bed within the fractures is not achievable. Simultaneously, the small proppants can be transported further under the effects of the fluid due to their low gravity. Large proppants are influenced by gravity, driving rapid deposition and short-distance transport, making them suitable for near-well propping. 4. 

The injection of mixed proppants will produce separation and form a cascading distribution for offshore shale (JRC_48, 52), which limits the migration of small proppants, resulting in a 25. 5% increase in the distribution height of the near-well proppant bed compared to the case of uniform proppants.
section_name: 3.

section_content: Injecting larger proppant particles followed by smaller ones reduces the migration distance by 9% and increases the near-well height by 17%. The optimal solution to enhance effective proppant support in the near-well zone of offshore shale fractures involves using an S:M:L combination, with d/w less than 0. 4 in the early stages of injection.
section_name: 5.

section_name: Conflicts of Interest:

section_content: Data Availability Statement: Data will be made available on request.

section_content: The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
section_name: Conflicts of Interest:

section_name: Funding:

section_content: Author Contributions: B. Y. : conceptualization, methodology, writing-original draft preparation. Y. L. : supervision, validation, project administration, funding acquisition. S. L. : investigation, data curation. P. X. : writing-review and editing. All authors have read and agreed to the published version of the manuscript.

section_content: The authors are grateful to the Hubei Key Laboratory of Oil and Gas Drilling and Production Engineering (Yangtze University) (YQZC202410) and CNOOC Research Institute Co., Ltd. (KJGG2022-0701) for the financial support of this work.
section_name: Funding:

section_content: Informed Consent Statement: Not applicable.
section_name: Institutional Review Board Statement: Not applicable.