Numerical simulation of proppant embedment in rough surfaces based on full reverse reconstruction

In the hydraulic fracturing process of shale reservoir, proppant will embed and deform under the action of high closure stress. Thus, it is necessary to analyze the mechanical process and influencing factors of proppant embedding in shale. This study reproduced rough fracture surface of real rock slabs based on the reverse reconstruction method and obtained the mechanical parameters of the rock slab and proppant. In this study, a numerical model for the elastoplastic deformation of the proppant embedding in the rough fracture surface, along with mechanical test experiments, is proposed. The reliability of the numerical model is verified by the proppant embedding simulation experiment. Based on this model, the process of proppant embedding in a rough fracture under the action of closing stress and the factors influencing the width of supporting fracture were analyzed. The results show that the embedment degree of the proppant is different in different areas of the rough fracture surface. Furthermore, a stress concentration effect is apparent. The proppant embedment becomes significant after the slab enters the plastic deformation stage under high closure stress. The slabs with high roughness and low Young's modulus or proppants with small particle sizes and high Young's moduli result in a smaller width of the propped fracture.


3-D Three-dimensional FD
Fractal dimension NURBS Non-uniform rational B-splines D-P criterion Drucker-Prager criterion C3D8R Hexahedral reduced integral element Change in width of the propped fracture 1 3 ΔW3 Variation of the piston and rock in the diversion chamber X i (k) Dimensionless mean value of the ith parameter i Parameters of the serial number k Data serial number The ith parameter value x i Arithmetic mean

Introduction
In recent years, with the gradual increase in shale gas exploration and development, the volumetric stimulation of "horizontal well + subsection + multi-cluster + large displacement + high sand volume" has become the leading technology for shale gas reservoir development (Ma et al. 2020;Zheng et al. 2020). The key to effective volumetric modification is the ability of the proppant to create a high conductivity channel during formation, and the embedment of proppant in the fracture under closure stress can substantially reduce the conductivity of the fracture. Therefore, the research on proppant embedment during proppant fracture is important for ensuring an efficient seepage of the shale gas reservoir (Bandara et al. 2019;Katende et al. 2021;Lei et al. 2018). Scholars mainly study proppant embedment through experimental and theoretical methods. The experimental methods rely on statistical theory of proppant embedment experimental results for constructing empirical models (Masłowski & Labus, 2021). Lacy et al. found that the embedment of high-strength proppant into a low-permeability reservoir rock reduced the formation fracture width significantly (Lacy et al. 1998(Lacy et al. , 1997. Wen et al. (2005) investigated the embedment degree of different formation cores with different sand paving concentrations and proppants. Lu et al. (2008) and Guo et al. (2008) conducted an embedment experiment with a self-developed proppant embedment tester and observed the embedment on the rock surface microscopically. Guo and Zhang, (2011) tested the embedment of resin-coated sand, ceramic and quartz sand in carbonate, and clastic rocks using an FCES-100 fracture conductivity instrument. Singh et al. (2019) investigated the effect of different types of fluid immersions on proppant embedment. In addition, other experimental studies of proppant embedment have highlighted the basic principles (Hu et al. 2014;Huang et al. 2019;Li et al. 2016;Alagoz, et al. 2021). However, the proppant embedment experiment is unrepeatable because the results are sensitive to test conditions. Thus, the influence of various factors cannot effectively be analyzed. Hence, theoretical methods are used for quantitative characterization of proppant embedment. The method establishes an analytical or numerical model based on the appropriate constitutive equation and necessary equilibrium conditions. Volk et al. (1981) derived an empirical formula for proppant embedment based on experimental results and studied the influence of proppant size and concentration on embedment depth. Guo and Liu, (2012) established a long-term proppant embedment analytical model considering rock elasticity and creep deformation through the viscoelastic theory. Li et al. (2018) established an analytical model considering the influence of proppant multilayer placement on fracture width based on elastic mechanics and the Hertz contact theory. Chen et al. (2018) used the elastic-plastic deformation stage to describe the embedment process and established an elastic-plastic analytical model of proppant embedment. Li et al. (2017) established a numerical model considering the influence of shale creep characteristics on proppant embedment using an ABAQUS software and showed that rock creep aggravates proppant embedment. Zhang et al. (2017) studied the numerical model of rock embedment with different proppant sizes using the computational fluid dynamics-discrete element method and highlighted that small proppant sizes aggravate the embedment. Mehmood et al. (2022) analyzed the embedment of rod-shaped proppant in the formation and found that higher recovery can be achieved using rod-shaped proppants than spherical proppants. Some of the assumptions in the theoretical models mentioned above are incomplete-the proppant is considered rigid, the fracture surface is considered smooth, and the actual roughness of the fracture surface is ignored.
In this study, the problem of proppant embedding in rough fracture surfaces of shale is discussed. Moreover, using the Longmaxi formation shale in southern Sichuan as an example, the mechanical properties of shale and proppant are tested through rock mechanics experiments. Based on the generalized Hooke's law and D-P criterion, a corresponding elastic-plastic constitutive equation is proposed, and the numerical model of proppant embedding in rough fracture surfaces is established. Furthermore, proppant embedding rough process and its influencing factors are analyzed.

Reverse construction of the rough fracture surface
The complete reverse reconstruction of a rough fracture surface has four stages: point-cloud data collection, point-cloud data processing, polygon processing, and surface processing.
Point-cloud data collection: Brazilian splitting method was used to shear rock slabs to obtain the surface topography of shale-supported fractures. A 3-D laser scanner was used to obtain the point cloud data of rough fracture surfaces (Lu et al. 2017) (Fig. 1a). The sampling spacing was 0.06 mm in the X direction and 0.05 mm in the Y direction.
Point-cloud data processing: The scanning interval of the 3-D laser scanner is inconsistent in the X and Y directions. The collected grid data points are not arranged in a regular grid and are easily lost in rough rock slab scanning. The kriging interpolation method was used to interpolate and regularize the fractured-point cloud data (Cai et al. 2013;Jin et al. 2003;Lu et al. 2021). The regional-variable value of the unsampled points in the region was calculated using the values of sampling points. A fracture plane with dimension of 35 mm × 130 mm was obtained by shearing, and the number of nodes in the X and Y directions were 25 and 100, respectively. Subsequently, the structured-point cloud data (Fig. 1b) were obtained.
Polygon processing: The structured point cloud data were 3-D coordinate data without a topological relationship; therefore, it was necessary to further reconstruct the surface features of point metadata. The ball-rolling method was used to generate triangular surfaces from point metadata (Vida et al. 1994) (Fig. 1c). The rolling-ball method is based on the principle that a sphere passes through a path that is simultaneously tangent to two surfaces, which form an envelope, and one of the envelope surfaces is a transition tangent contact line; thus, a transition surface is formed.
Surface processing: The surface roughness of the supporting fracture is high. Thus, it is necessary to achieve precise regional division on triangular surfaces. For areas with high surface roughness, the surface pieces are divided more precisely to generate non-uniform rational basis spline (NURBS) surfaces with high precision in full reverse reconstruction software (Fig. 1d) and realize a fully reverse 3-D reconstruction of rough fracture surfaces.

Slab and proppant modeling
The NURBS surface was projected on the X and Y planes and stretched along the Z-axis, and the redundant part was removed by a Boolean operation to obtain the processed rough fracture model (Fig. 2a) A single proppant particle was treated as an ideal sphere, ignoring the influence of sphericity and surface smoothness. The distribution of proppant particles in rough fractures was assumed to be tightly layered. The breakage of the proppant particles was not considered (Fig. 2b).

Model constitutive equation
The mechanical constitutive equation of proppant embedment in rough fractures includes the elastic-plastic  The elastic stage of the proppant embedment and proppant deformation processes are described using the generalized Hooke's law (Liu and Luo, 2004). In finite element analysis, a matrix form is generally used for characterization. The slabs and proppant are treated as isotropic materials, and breakage is not considered. The elastic constitutive equation can be as follows.
Where [D] e is the three-dimensional elastic matrix, { } is the element stress matrix, { } is the element strain matrix, E is the Young's modulus of the material, and v is the Poisson's ratio of the material.
The plastic stage of the proppant embedment process is described by the D-P criterion (Fan et al. 2021;Öztekin et al. 2016): where I 1 is the first principal invariant of the stress tensor, J 2 is the second principal invariant of the deviatoric stress tensor, α and k are constants related to the cohesion and internal friction angle, φ is the angle of internal friction, and c is the cohesion. (1) Different values of θ correspond to different values of α and k, representing different circles reflected in the π plane and are applicable to different rock states.

Model material parameters
The numerical simulation of proppant embedment represents the mechanical process of contact deformation between the proppant particles and the rough rock slab. It is necessary to set the rough-rock slab and proppant material mechanical parameters, which are derived from the uniaxial rock mechanics test of the shale. The rock samples are taken from the shale core of the first member of the Longmaxi Formation, Weiyuan (No. 1#~4#). The experimental results are listed in Table·1. The deformation parameters in the plastic stage adopt the corresponding values in the stress-strain curve (Fig. 3).
The proppant mechanical parameters are derived from the radial compression experiment of a single-particle proppant conducted by Man and Chik-Kwong Wong, (2017) (Fig. 4). The Young's modulus and the Poisson's ratio of the proppant materials were calculated to be 13 GPa and 0.2, respectively.

Meshing and boundary conditions
The proppant unit was dominated by a hexahedral linear reduction integral unit, which is meshed by an advanced algorithm and sweep. The mesh type of the rock plate was C3D8R, and the cell shape was hexahedral. The high-precision sweep mesh division technique and the neutral axis algorithm were adopted. Furthermore, the linear seed was distributed on the interface perpendicular to the rock plate to simplify the calculation.
For the upper slab, symmetrical displacement constraints were applied in the X and Y directions, and a normal stress in the Z direction was applied in the upper section. A symmetrical displacement constraint was applied to the bottom face to ensure that it is completely fixed. Lastly, X and Y horizontal constraints were applied to all proppant particles.

Numerical simulation
To obtain the variation of the fracture width after proppant embedment in the rough slabs under the loading conditions, this study adopted the method of setting selected points on the rough fracture surface to output the displacement of the process. Then, the statistical fracture width between the whole upper and lower slabs was obtained.
(1) After meshing, a 100-pinpoint set of shale plate contours with rough surfaces and proppant contacts was established for the model.
(2) The displacement values of the output-point sets under different load increments were calculated statistically using integration.
where W f is the width of the propped fracture; n is the nth pair of computation nodes, n = 1, 2, 3, 100; and h n is the spacing of the nth pair of computation nodes.
Using equation (9), the change in fracture width after proppant embedment in the slab under different loads was obtained.

Laboratory experiment and comparison of results
In the propped fracture embedment modeling experiment, a fracture conductivity test device was used to test the variation of proppant embedment fracture width under different closure stresses (Fig. 5). The results were compared with the calculation results of the propped fracture embedment numerical model to verify the accuracy of the model.
In this study, the gaps between the pistons and the upper and lower shale slabs are completely closed by cyclic compression. The variation of the piston and the rock in the guide chamber is represented by the change in displacement of the loaded shale slabs. The fracture width is calculated using the following equation where ΔW is the displacement change in the sensor; ΔW 1 is the change in the displacement of the shale rock, ΔW 2 is the change in width of the propped fracture, and ΔW 3 is the variation of the piston and rock in the diversion chamber.
When the load is stabilized at 10, 20, 30, 40, 50, and 60 MPa, the displayed value of the displacement sensor is recorded to obtain the overall displacement variation of the diversion chamber under different loads. The variation curves of fracture width obtained from numerical simulation and experiment are shown in Fig. 6. The variation trends of propped fracture width obtained by numerical simulation and physical model experiment become approximately similar with increasing closing stress. When the closure stress is lower than 40 MPa, the calculation error of the model is approximately 10%. The error at this stage is caused by the micro-bumps on the rough fracture surface under experimental conditions being gradually compacting under the closure stress. When the closure stress is higher than 40 MPa, the error occurs because the numerical model does not consider the proppant crushing, resulting in a change of the propped joint width tested in the experiment being larger than the numerical simulation results. The calculation error is less than 16 %, indicating that the numerical model is accurate and reliable.

Process analysis of proppant embedding in rough fractures
The stress cloud of the proppant embedment model under different closure stresses is shown in Fig. 7. Due to the roughness of the fracture surface, the contact between the proppant particles and the fracture surface is different. Thus, the stress in different regions of the fracture surface is different. With an increase in the closure stress, the stress region is gradually transferred to the interior of the rock plate.
At a closure stress of 60 MPa, the autoshapes of different proppant particles differ due to the different contact stresses applied to the particles (Fig. 8a). The total proppant shape was calculated by recording the coordinates of the upperand lower-end points of the proppant particles. Combined with the analysis of the proppant embedment depth, the proppant shape contributed a small amount of 22% to the proppant width reduction. The stress distribution on the fracture surface is uneven due to different proppant embedment amounts in different zones (Fig. 8b). The roughness of the contact profile before and after embedding is characterized by the joint roughness coefficient. The roughness of the fracture surface increases after embedding.
The proppant under different closing stresses was embedded into the mesh model for materialized processing. A rectangular entity with the same length, width, and height was rebuilt to perform a Boolean operation of the two volumes to obtain the embedded fracture space (Fig. 9). The change in propped fracture void volume was obtained by loading different closing stresses (Fig. 10). With an increase in the closure stress, the deformation of the proppant and embedment of proppant into the rough fracture surface led to the decrease of flow channels in the upper and lower slabs, resulting in the continuous reduction of fracture conductivity (COOK and Gw, 1992). When the closure stress was 60 MPa, the propped fracture void volume decreased by 60% compared with that when the closure stress was 10 MPa.

Analysis of influencing factors
The process of coarser proppant embedded in shale is mainly affected by rock properties (fracture surface roughness, Young's modulus, and Poisson's ratio) and proppant properties (Young's modulus, particle size, and arrangement). In this study, the control variable method is adopted for analysis, and the simulation parameter settings are listed in Table·2.
The roughness of the shale fracture surface is complex. In this study, a fractal dimension is used to quantitatively characterize the fracture surface, and a control experiment is carried out with a smooth rock. From Fig. 11a, the larger the fractal dimension of the fracture surface, the larger the initial fracture width and the greater the variation of the propped fracture width after embedding. This is because as the protrusions on the rough fracture surface contact with the proppant reach a stress concentration, the protrusions yield deformation, resulting in a significant decrease in the propped fracture width.f Figure 11b shows the influence of the Young's modulus of the rock on the width of the supporting joints. When the closure stress is lower than 40 MPa, the rate of change of the propped fracture width decreases with an increase in the Young's modulus. When the closure stress is higher than 40 MPa, the width of propped fracture varies and the closing stress has the same variation, indicating that the Young's modulus of the rock affects the elastic deformation stage of proppant embedment.
The influence of the rock's Poisson's ratio on the propped fracture width is shown in Fig. 11c. The decreasing trend of the propped fracture width with the increase in closing stress is approximately the same under the three Poisson's ratio conditions. The final variation in propped fracture width is small, which indicates that the influence of the rock's Poisson's ratio on the propped fracture width is small.
The influence of the proppant's Young's modulus on the propped fracture width is shown in Fig. 11d. With an increase in the Young's modulus of the proppant, the variation in proppant fracture width decreases. This is because at low proppant Young's moduli, the change of proppant fracture width is affected by proppant deformation.  The maximum size of the proppant particle size was obtained. The 20/40 mesh, 30/50 mesh, and 40/70 mesh diameters were 0.85 mm, 0.60 mm, and 0.424 mm, respectively. The proppant fracture width also decreased as the proppant particle size decreased. Because the contact area of the smaller-particle proppant is smaller, the contact stress increases, increasing the proppant embedment and reducing the proppant fracture width (Fig. 11e).

Sensitivity analysis
The influence degree of each influencing factor on the support seam width is different, and it is necessary to conduct a sensitivity analysis on the influencing factors. In this paper, the grey relational analysis method is used to determine the influence weight of each factor on proppant embedding (Hu, 2003). The specific implementation steps are as follows: (1) Perform dimensionless processing: (2) Solving for the grey relational coefficient: (3) Solve the correlation degree: (4) Solve for the weight coefficient: where X i (k) is the dimensionless mean value of the ith parameter, i is the parameters of the serial number, k is the data serial number, x i (k) is the ith parameter value, − x i is the arithmetic mean, i (k) is the grey relational coefficient, is the distinguishing coefficient, i is the correlation degree, W i is the weight coefficient. Figure 12 shows that the factors affecting proppant embedment are proppant particle size, proppant Young's modulus, surface roughness, rock Young's modulus and rock Poisson's ratio in descending order of weight. Therefore, various construction parameters can be optimized in the construction process to alleviate the damage of proppant embedding to propped fractures and improve the fracture conductivity.

Conclusions
In this study, a 3-D numerical model of proppant embedment in rough fractures was established based on a fully reverse reconstruction method. The model was verified by proppant embedment experiments, and the effects of rock roughness, mechanical properties, and proppant mechanical properties on the embedment process were analyzed. The following conclusions are drawn from this study.
(1) Based on the full reverse reconstruction method to replicate the real surface topography of shale fractures, the slab and proppant models of rough surface were established.
(2) Based on the generalized Hooke's law and D-P yield criterion, the proppant deformation and contact deformation models between rough slabs and proppant particles were established, and the validity of the numerical i models was verified by proppant embedment experiments.
(3) Under closure stress, the degree of proppant embedment in different areas of the rough fracture surface of shale was different, there is an apparent stress concentration effect, and the embedment made the fracture surface coarser. At a closure stress of 60 MPa, proppant deformations contributed approximately 22% to the fracture width variation, reducing the fracture space by approximately 60% compared with initial contact. (4) The proppant fracture width with a high roughness, high rock's Young's modulus, high proppant's Young's modulus, and large proppant particle size has little variation, and the rock's Poisson's ratio has little influence on proppant embedment process.
Author contribution CL was involved in the conceptualization and methodology; JL helped in the data curation and writing-original draft; FH contributed to software and visualization; JG was involved in the conceptualization; JW contributed to the methodology and supervision; GZ was involved in the validation and formal analysis; JW contributed to the methodology and software; XM was involved in the formal analysis and writing-review & editing; YL helped in writing-review & editing; XW contributed to the literature research and resources; XS helped in the literature research; HL was involved in drawing graphs.
Funding This work was supported by National Natural Science Foundation of China [Grant No. 52022087].

Conflict of interest
The authors declare that they have no known competing financial interest or personal relationships that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.