DEM-CFD Modeling of Proppant Pillar Deformation and Stability during the Fracturing Fluid Flowback

In this study, proppant pillar deformation and stability during the fracturing fluid flowback of channel fracturing was simulated with DEM-CFD- (discrete element method-computational fluid dynamics-) coupling method. Fibers were modeled by implementing the bonded particle model for contacts between particles. In the hydraulic fracture-closing period, the height of the proppant pillar decreases gradually and the diameter increases as the closing stress increases. In the fracturing fluid flowback period, proppant particles could be driven away from the pillar by the fluid flow and cause the instability of the proppant pillar. The proppant flowback could occur easily with large proppant pillar height or a large fluid pressure gradient. Both the pillar height and the pillar diameter to spacing ratio are key parameters for the design of channel fracturing. Increasing the fiber-bonding strength could enhance the stability of the proppant pillar.


Introduction
Channel fracturing is a relatively new stimulation technique first proposed by Gillard et al. in 2010 [1].It is mainly composed of three technical parts: the pulse pumping technique, the multicluster perforating process, and the injection of fracturing fluid mixed with fibers.The major difference between channel fracturing and conventional fracturing lies in the pattern of proppant placement.In channel fracturing, the fibers are expected to keep proppant pulses cohesive and prevent them from spreading when traveling through the fracture slots.As a result, the proppant after channel fracturing could assemble as clusters in the fracture.When the fractures close during fracturing fluid flowback or the production process, those proppant clusters are analogous to pillars which resist the closing process.The void space between proppant pillars form the fluid channels which greatly improve the fracture conductivity during the production [2].An investigation based on the results of more than 1000 times channel fracturing found that more than 99.9% of the channel-fracturing jobs fully completed the proppant placement and by average channel fracturing could save 43% of the proppant compared with the conventional technique implemented in adjacent wells [3,4].The oil production with channel fracturing could also be greatly enhanced.For instance, the channel-fracturing operation of tight oil and gas reservoirs in Ordos Basin, China, has produced 2.4 times as much oil as conventional fracturing, and the gas well production is 4-5 times as high as that of conventional fracturing [5].
After the hydraulic-fracturing operation, the fracturing fluid flows back to the wellbore which results in additional drag force to the proppant particles.Consequently, the proppant pillar might lose its stability due to the flow of proppant, resulting in sanding and significant decrease of fracture conductivity.Therefore, how to avoid the flowback of proppant has attracted great attention.Channel-fracturing technology adds fibers to the fracturing fluid by which the rheological properties of the proppant are changed.The fracturing fluid mixed with low concentration short fibers can be affected by the formation temperature and exhibit some adhesion [6,7].Figure 1 shows the network structure of fibers and proppant mixture.The bonding strength of fibers could effectively glue the proppant particle together, and thus, the fiber-proppant cluster is able to better support the hydraulic fractures, prevent the proppant flowback, and also reduce the risk of sanding.In the channel-fracturing tests of Eagle Ford and other shale fields, sand production was improved and fibers showed better proppant flowback resistance [5,8,9].
The network structure of fibers and proppant mixture and the adhesion of fibers make the mechanical mechanism of the proppant pillar greatly different from the one for the traditional proppant placement scheme.The existing literatures mainly focus on the mechanical process of the flowback of fiber-free proppant.Vreeburg et al. [10] performed a laboratory study on two types of back-production to help clarify the effect of curing temperature, water production rate, proppant size, and stress cycling on the integrity of resin-coated proppant (RCP) packs.The number of applied stress cycles and the initial RCP pack strength appear to be the dominant factors that govern proppant back-production.Goel and Shah [11] investigated proppant flowback phenomena using a large-scale fracturing simulator.Experimental results show that flowback initiates at lower cleanup rates when the closure stress increases or when the fracture width increases and flowback initiates at higher cleanup rates when the sand size increases.Smith et al. [12] investigated the effects of various factors on the proppant flowback, including the differential-fracture closure, leak-off of fracturing fluid, and fracturing fluid rheology, and claimed that those effects should be combined in the proppant transport analysis.Aidagulov et al. [13] presented a quantitative model to predict proppant flowback based on treating both the proppant pack and the reservoir as poro-elasto-plastic media.But with adoption of such a continuum model, interactions between flowback proppant particles would not be observed directly.Daneshy [14] suggested that the main cause of the proppant flowback was the presence of shear fractures which led to the formation of many randomly distributed tight proppant packs.The three requirements for proppant flowback were motion initiation, motion maintenance, and infinite conductivity along the return path, and gravity played a very important role on proppant flowback.Hu et al. [15] established a proppant mechanical model to predict critical flowback velocity.They concluded that adopting a lower velocity before the fracture closes and a higher velocity after the fracture closes helps to discharge fracturing fluid in time and effectively prevent proppant flowback.Qi and Jiao [16] presented a model to predict proppant flowback in a fractured gas well.However, with the hypothesis of a one-dimensional flow and ignoring the influence of gravity, this model cannot demonstrate real interactions between liquid and proppants.McLennan et al. [17] described the experimental results of a radially convergent flow and a linear flow into a wellbore through an initially packed fracture.Results of a radial flow through a transverse fracture suggested the creation of flow channels, while results of a linear flow through a longitudinal fracture suggested more proppant removal.The above studies have pointed out that effects such as fluid viscosity, fluid pressure gradient, and closing stress are crucial to the proppant flowback for the conventional proppant placement scheme.They laid an important foundation for the research of proppant flowback with channel fracturing.
Other studies explored the evolution of fracture width and conductivity in channel fracturing.The deformation characteristics of the proppant pillar are described using numerical and analytical methods though some of them are not quite reasonable.Gomaa et al. [18] described experimental and numerical models to generate pillar-propped fractures based on fingering phenomena observed in fluid injection.Experimental and numerical results confirmed that increasing the injection rate reduced the main channel width and increased branching.Neto and Kotousov [19] developed a simplified semianalytical method for calculating the residual opening of fractures partially filled with proppant.They used the one-dimensional Terzaghi's consolidation model to describe proppant response.Both pore water and particles are assumed to be incompressible and changes of void volume are linked to deformation.Zheng et al. [20] derived an analytical model to calculate permeability and conductivity in channel fracturing.The proppant pillar is treated as a cylindrical indenter which is rigid.Guo et al. [21]  Later, Deng et al. [23] adopted the discrete element method to study interactions between shale and proppant.Effects of factors such as shale modulus and proppant size on fracture aperture were numerically modelled and learned.Fan et al. [24] developed a DEM-LB workflow to understand the interaction between reservoir depletion, proppant compaction, and single-/multiphase flows in a hydraulic fracture.2 Geofluids Furthermore, Zhang et al. [25] developed an integrated DEM-CFD modeling workflow to model proppant embedment and fracture conductivity.Fracture conductivity after fracture closing was evaluated by modeling fluid flow through the proppant pack by use of DEM coupled with CFD.Results showed that the fracture conductivity increased with the increase of proppant concentration or proppant size and decreased with the increase of fracture-closing stress or degree of shale hydration.Shale-hydration effect was confirmed to be the main reason for the large amount of proppant embedment.With increasingly frequent application of DEM-CFD, it was believed that DEM-CFD was the most suitable computational method for modeling two-way solid-fluid interactions (e.g., Tomac and Gutierrez [26]).This paper utilizes the DEM-CFD-coupling method, by which the cementation between the fibers and the proppant particles is mimicked by the tensile and shear bonding strength between discrete element particles, to set up the hydromechanically coupled model for the interactions of proppant-fiber-fracturing fluid.The model aims to reveal the mechanism of proppant pillar deformation and stability during fracturing fluid flowback by combing the analysis at both micro-and macroscales.The effects of pressure gradient, fluid viscosity, pillar height, pillar diameter to spacing ratio, and bonding strength of fibers are investigated, and the amount of flowback proppant and the spreading area of proppant are studied.This work could provide a potential guideline and theoretical background for the design of channel fracturing and also the optimization of field fracturing operation.[27] first studied the stability of a proppant pack reinforced with fibers and presented a fiber-reinforced failure criterion in their study.They set up a laboratory model to demonstrate the effect of fiber properties on proppant pack stability.Two modes of failure were distinguished in their experiments: collapse of the arch at the perforation and total failure of the proppant pack creating a channel within the fracture.The fracturing fluid used in the channel fracturing contains fibers that have wrapping and restraining effects on proppant, which can enhance the sand-carrying capacity of the fracturing fluid.When proppant and fibers are injected into the fractures with fracturing fluid, the bonding and friction between the fibers and proppant make the proppant cluster as a whole, which achieves the balance under the action of closing stress and lateral stress.During the flowback, the proppant pillar is subjected to shear deformation by the fluid force; however, the existence of fibers helps to prevent the collapse of the proppant pillar.

Theoretical Background and Numerical Model Setup
In this study, we apply the bonded particle model in DEM to model the interaction between proppant particles and fibers [28,29].DEM is an ideal tool to model the mechanical behaviors of a granular assembly such as a proppant pack.The macroscale mechanical response emerges from the interaction of particles or grains through their contact behaviors at the particle scale, either frictional or with cementation.Since the shape of fibers is irregular, it is technically difficult to model fibers together with proppant particles.Thus, we indirectly model the fibers by implementing the bonded particle model which adds bonding strength between proppant particles but without explicitly representing the fibers.The theory of DEM and bonded particle model are well documented in the literatures and hence are not presented here for brevity.

DEM-CFD Coupling for Proppant Particle-Fracturing
Fluid Interaction.By considering force balance for particles within a cell, the driving force applied to a single particle is given by (1).The driving force consists of two kinds of forces caused by the flow.One is the viscous force induced by fluid viscous properties and characterized by shear force on the particle surface.The other is the force acting on the particles by pressure gradient and characterized by a normal force on the particle surface [30].
where f d i is the driving force applied to particle i, f int is the interaction force per unit volume between particles and fluid, d pi is the diameter of particle i, ϕ is the porosity of this cell, and ∇p is the pressure gradient.
In the first term of (1), f int represents the interaction force per unit volume between particles and fluid.The minus sign means that the force applied to the fluid is opposite in direction with f int .The sign of f int depends on the relative moving direction of particles in flow.If the relative moving direction is opposite between particles and the flow, f int is negative and the driving force is positive, which means particles are pushed forward by the flow.f int equals 0 when particles are static relative to the flow.In the second term of (1), the minus sign means that pressure decreases with positive flow when ∇p is negative.The force acting on the particles by the pressure gradient pushes particles forward if ∇p is negative.
When fracturing fluid begins to flow back, the flow direction is opposite to the relative moving direction of proppants and pressure decreases with positive flow.Thus, both viscous force and the force acting on the particles by the pressure gradient work as positive forces to particles.
Interaction force per unit volume between particles and fluid in the j direction is given by where u j is fluid velocity in the j (j = x, y, z) direction, v j is average velocity of particles in the j direction, and β int j is the fluid-particle friction coefficient.β int j is defined by 3 Geofluids different expressions according to different porosity values.When porosity is lower than 0.8, β int j is given by where d p is the average diameter of particles, μ f is the kinematic viscosity of fluid, and ρ f is the density of fluid.When porosity is higher than 0.8, β int j is given by where C D is the drag coefficient (dimensionless).
Pressure gradient ∇p j also has two expressions accordingly.When porosity is lower than 0.8, ∇p j is given by Ergun's equation: When porosity is higher than 0.8, ∇p j is given by the Wen-Yu equation: 2.2.Numerical Model Setup. Figure 2 shows the numerical model setup.The model consists of one proppant pillar which is sandwiched by two shale plates.The rock plate and proppant pillar are made of blue particles with a facecentered cubic (FCC) structure and yellow particles with random distribution, respectively.Taking the 6 mm high proppant pillar as an example, the simulation procedure can be summarized in the following steps.
(1) Generate a cubic-shape rock sample with the FCC structure.The initial model size is 20 mm × 20 mm × 30 mm.The microscale properties of rock are assigned to all contacts.
(2) A small confining stress of 0.5 MPa is applied to the sample by a servocontrol procedure.
(3) A 6 mm thick layer in the middle of the sample is deleted to generate a fracture slot perpendicular to the z direction.
(4) A cylindrical proppant pillar with a height of 6 mm and diameter of 10 mm is generated in the middle of the fracture slot.
(5) The two rock plates are gradually loaded with the compression stress of 41.4 MPa while the proppant pillar is deformed.
(6) After the loading, the fluid grid is set up for the proppant layer and fluid flow calculation is carried out.
For the calculation of fluid-particle coupling, 10 fluid cells are set up along the two horizontal directions and 1 fluid cell is set up along the vertical direction of the proppant pillar height.In order to apply the fluid boundary condition, one additional layer of fluid cells is needed at each of the six faces of the calculation domain.It means that in each of the three directions (x, y, z), there are two additional boundary layers.The final cell pattern is 10 + 2 × 10 + 2 × 1 + 2 = 12 × 12 × 3 given that the initial cell pattern is 10 × 10 × 1.Therefore, the gird mesh has 12 × 12 × 3 = 432 fluid cells in total.With the two additional boundary layers, the final model size in both the x and y directions is 10 + 2 × 2 mm = 24 mm (Figure 3).Cells outside the red line in Figure 4 are boundary cells.Only cells within this "boundary" ran the fluid-particle interaction computation.
In fluid flow calculation, a negative pressure gradient is applied along the x direction by applying positive fluid pressure on the left boundary at the inlet and applying zero fluid pressure on the right boundary at the outlet.The flowback pattern of the proppant pillar at different times can then be simulated.(1) Put proppant particles and fibers (weight ratio 1 : 0.004) into a beaker.

Calibration of the
(2) Use a glass rod to mix the proppant and fibers (Figure 5(a)).
(3) Add glue into the proppant-fiber mixture and stir the sample until it is solidified.(9) Repeat step (1) to (8) for a different closing stress.
To calibrate the microscale properties of proppant pillars, a numerical sample of the proppant pillar was first made and then compressed with constant stress from the two ends (Figure 6(c) and 6(d)).Numerical tests were carried out with  5 Geofluids different compression stresses, and the pillar height and pillar diameter after the tests were compared with the experimental results for the model calibration, as shown in Figures 7 and 8.The final calibrated microscale properties of proppant particles are listed in Table 2.

Proppant Pillar Deformation and Damage
Laws during Fracture Closing and Fracturing Fluid Flowback      8 Geofluids and form a single layer.After 0.0105 s, the proppant particles gradually escaped from the right side of the flow field and start to flow back.The number of flowback particles in this paper is defined as the number of proppant particles escaping from the flow field.
We select the first row of fluid grids at the entry (shown in Figure 11(b)) to monitor the change of flow velocity.Fluid velocities reach stable values at about 0.04 s flow time (Figure 12).Before 0.005 s, fluid velocities show fluctuation which indicates that the flow is subjected to unstable change of proppant ahead due to the spreading of particles.Between 0.025 s and 0.035 s, fluid velocities show more drastic fluctuation caused by fast flowback of a number of proppants.The results above indicate that a larger pillar diameter is prone to cause proppant particle spreading as well as unstable fluid flow.On the other hand, a smaller pillar diameter will lead to wider fluid channel and favor the stable fluid flow; however, the pillar-supporting capacity might not be met to resist the fracture-closing stress.Therefore, a proper size of proppant pillar is needed which can satisfy both the stability requirement and conductivity enhancement.10 Geofluids different fluid pressure gradients is minor.At the second stage, the difference of the number of flowback proppant particles gets larger and larger.Note that the difference of spreading area for different fluid gradients is minor.The reason is that the calculation of spreading area is based on the farthest particle in the domain.Figure 14 plots the evolution of the fluid velocity at the monitoring fluid cell (red cell in Figure 11(b)) for different pressure gradients.As the calculation starts, the fluid velocity gradually increases and reaches the plateau value which is proportional to the pressure difference.The difference of fluid velocity for the first stage of the flowback is small which accounts for the small difference of flowback proppant particles as shown in Figure 13.
Figure 15 plots the evolution of the proppant pillar profile for different pressure gradients.In order to make better comparison for pillar profiles, proppants which run outside the flow field are not shown and similar figures are treated in the same way.For different pressure gradients, the profiles of proppant pillars are similar at the first stage.At the second stage, more proppants flow back in the high-pressure gradient case than in the low-pressure gradient case.

Effect of Fluid Viscosity.
To study the effect of fracturing fluid viscosity on the flowback of proppant particles, the proppant pillar with 6 mm height and 10 mm diameter was run with different fluid viscosities, while the pressure difference between the inlet and outlet of the flow field was fixed as 64 kPa. Figure 16 shows the evolution of the number of flowback proppant particles and proppant spreading area for different fluid viscosities.It can be seen that when the pressure boundary condition is used, the amount of flowback particles and the spreading area decreases with the increase of fluid viscosity.According to (5), with the constant pressure gradient, increase of the fluid viscosity causes the decrease of the velocity difference v j − u j .However, from (3), the size of the β int j is difficult to determine (one part of (3) increases and the other decreases).Therefore, the force exerted on the particles by the fluid is difficult to determine.Figure 17 plots the evolution of fluid velocity for the monitoring fluid cell (blue cell in Figure 11(b)) for different fluid viscosities.It shows that when the viscosity of the fluid is large, the overall flow rate is small, which indirectly explains the phenomenon that the number of particles and the spreading area are smaller when the viscosity is larger.Figure 18 plots the evolution of the proppant pillar profile for different fluid viscosities.For lower viscosity cases, proppant pillars break down into separated packs, while for higher viscosity cases, proppant pillars are able to retain integrity during the flowback period.
It should be noted that the trends in Figure 16 will be reversed if velocity boundary is applied at the inlet.Figure 19 plots the evolution of the number of flowback proppant particles and proppant spreading area for different fluid viscosities, with a fixed fluid velocity of 3.64 m/s applied at the inlet.The results show that larger fluid viscosity causes more flowback proppant particles and larger spreading area, opposite with the results in Figure 16.At the beginning of calculation, forces applied by the fluid on particles varied due to different particle sizes.Particle movement was further complicated by collision, and thus, average particle velocity might change dramatically.This kinematical disorder in turn caused fluctuation at the beginning of the blue line by particle-fluid interaction term in (2).
Figure 20 plots the evolution of the proppant pillar profile for different fluid viscosities, with velocity boundary applied at the inlet.Since the flowback process is different for the pressure boundary case and velocity boundary case, pillar profiles at different flow times are used.For the velocity boundary case, proppant pillars are more possible to break down when fluid viscosity is higher, which shows an exactly reverse law to the pressure boundary case.11 Geofluids 4.3.Effect of Proppant Pillar Height.With the given proppant pillar diameter, an initially tall proppant pillar is expected to provide larger closing aperture and more conductive channels for the fluid flow after fracture closing than the initially short proppant pillar.However, the tall proppant pillar might also be subjected to the large amount of flowback proppant particles since during the fracture-closing process, particles at the edge of the proppant pillar are prone to be stripped off. Figure 21 shows the proppant pillar profile after the flowback for two different proppant pillar heights, 6 mm and 8 mm.The case with the 8 mm pillar height has larger spreading area and more flowback proppant particles than the case with the 6 mm pillar height.
The lab experiments by Gillard et al. [1] and Nguyen et al. [31] show that the conductivity of proppant pillars can be up to two orders of magnitude higher than the one with the traditional proppant placement scheme.The proppant pillar profile after the conductivity measurement (Figures 22 and 23 flowback of proppant particles cannot be clearly identified.
Our modeling results however show that the flowback could occur with a large proppant pillar height or large fluid pressure gradient.As a result, the conductivities measured in the numerical models could be largely dropped due to the occupancy of particles in the fluid channels.The field tests with channel fracturing shows the increase of oil/gas production [3][4][5] compared with the traditional proppant placement scheme, but the improvement is much less remarkable than the experimental results by Gillard et al.
Due to the flowback of proppant particles, the proppant particles move into the channels between pillars, which results in the single layer or a few layers of proppant particles between pillars (Figure 21).Part of a proppant pillar could still remain stable after flowback, and the initially open fluid channels might be filled with loosely packed proppant particles.The work by Fredd et al. [32] shows that the single-layer proppant in the fracture could enhance the fracture conductivity.Therefore, the fluid channels between pillars, though filled with a single layer or a few layers of proppant, may still effectively increase the fluid conductivity and make channel fracturing preferable than conventionalfracturing methods.13 Geofluids Among them, d is the diameter of the proppant pillar, x is the pillar spacing, L is the length of the fracturing fluid without proppant, and D is the length of proppant-laden fluid.
The pillar diameter to spacing ratio (d/x) is therefore a key parameter for the design of channel-fracturing treatment.If the pillar diameter to spacing ratio is small, the pillar column might not be stable which causes small fractureclosing apertures and the fracture conductivity is reduced.On the contrary, large d/x leads to a narrow fluid channel and increasing amount of flowback particles.To optimize the value of d/x and determine the appropriate volume ratio of two fracturing fluids, simulation series by fixing the pillar spacing but varying the pillar diameter is carried out.
Figure 25 plots the percentage of flowback proppant particles for different d/x by using the pillar height of 6 mm and diameter of 10 mm.The fluid viscosity is 0.001 Pa•s, and the pressure difference is 64 kPa.There is nearly no flowback particles when the d/x ratio is smaller than 0.43.The percentage of flowback proppant particles reaches 57% when the d/x = 0 70.If the d/x ratio is more than 0.70, the percentage of flowback proppant particles is close to 100%.The results suggest that the optimized d/x ratio is between 0.43 and 0.70, with the assigned input parameters in this case.However, it is worth mentioning that the analysis above is based on the assumption that the pillars are placed in arrays, while the real distribution of proppant pillars in field operation is expected to be less regular, given the fact that there could be multiple circumstances which cannot all be considered in the numerical models.
Figure 26 plots the evolution of the proppant pillar profile for different pillar diameter to spacing ratios.With increasing  4.5.Effect of Fiber-Bonding Strength.The effect of fiberbonding strength on the stability of the proppant pillar was also studied.Figure 27 shows the evolution of the number of flowback proppant particles and proppant spreading area for different fiber-bonding strengths, with the proppant pillar height of 6 mm and the diameter of 10 mm.The fluid viscosity is 0.001 Pa•s, and the pressure difference is 64 kPa.The greater the bond strength, the stronger the cohesive effect of the simulated fibers on the proppant particles.With the enhancement of cohesion, the number of flowback particles from the proppant pillar is reduced and the spreading area is also decreased.Figure 28 plots the evolution of the proppant pillar profile for different fiberbonding strengths.Comparing with bonding strengths of 3E + 3 and 3E + 6 Pa, the proppant pillar remains quite stable during the flowback period when the bonding strength is 3E + 7 Pa.It can be seen that reasonably adjusting the physicochemical properties of the fibers increases the cohesion between the proppant particles and helps to enhance the stability of the proppant pillar in the highpermeability fracture.

Conclusions
In this study, proppant pillar deformation and stability during the flowback of channel fracturing was simulated with the DEM-CFD-coupling method.The effect of fibers was considered by implementing the bonded particle model for contacts between particles.The major conclusions of this study can be summarized as follows.
In the fracturing-closing period, the height of the proppant pillar decreases gradually and the diameter increases  as the closing stress increases.In the flowback period, proppant particles could be driven away from the pillar by the fluid flow and cause the instability of proppant pillar.
The flowback could occur easily with a large proppant pillar height or large fluid pressure gradient.Increasing the initial proppant pillar height gives more conductive channels for the fluid flow after fracture closing, but on the other hand, a tall proppant pillar might also be subjected to the large amount of flowback proppant particles.
Besides the proppant height, the pillar diameter to spacing ratio is another key parameter for the design of channel fracturing.A small pillar diameter to spacing ratio causes the instability of the proppant pillar and small fractureclosing aperture; on the contrary, a large pillar diameter to spacing ratio leads to narrow fluid channels and increasing amounts of flowback particles.
The profile of proppant placement inside the fracture in the field should be somehow between the ideal competent proppant pillar arrangement and the uniform layer of the proppant based on the traditional proppant placement scheme.Adjusting the physicochemical properties of the fibers increases the cohesion between the proppant particles and helps to enhance the stability of the proppant pillar in the high-permeability fracture.
established an analytical model to describe fracture aperture change and conductivity for cuboid shaped pillars.Yan et al. [22] developed an analytical model to represent the physical deformation of channel-fracturing fractures, and the Darcy-Brinkman equation was applied to simulate the flow in pillars and fluid channels.Both Guo et al. and Yan et al. calculate the axial deformation of the proppant pillar using Hooke's law.

Figure 2 :
Figure 2: (a) Model perspective view and (b) coupled DEM-CFD analysis for proppant pillar deformation and stability.The blue and yellow particles represent the rock matrix and proppant pillar, respectively, and the green meshes represent the CFD fluid grids.

( 4 )( 6 )( 7 )( 8 )
Fill the sample into a metal mold (hollow cylinder with inner diameter of 10 mm and height of 10 mm) and make proppant pillars (Figure 5(b)).(5)Heat the proppant pillars in an oven for 1 hour with the temperature of 60 degrees, then cool them down under ambient condition.Fill the array of proppant pillars in the testing chamber designed based on the API standard for proppant conductivity testing; load the pillars at the given normal loading stress and keep the loading for 1 minute until reaching the steady state (Figure 5(c)).Calculate the change of pillar height based on the recorded variation of vertical displacement.Disassemble the experimental setup and record the profiles of proppant pillars after the testing (Figure 6(a) and 6(b)).

Figure 3 :
Figure 3: Illustration of model size: (a) grid mesh without particles; (b) size of boundary cells and flow area.

3. 1 .
Fracture-Closing Period.The calibrated microscale properties (Tables1 and 2are used as input parameters in the model shown in Figure2.Numerical modeling of fracture closing and fracturing fluid flowback is then calculated.Figures 9 and 10 plot the profiles of a proppant pillar at two different closing stresses, 20.7 MPa and 41.4 MPa, respectively.The height of the proppant pillar is 6 mm, and the diameter is 10 mm.With the increase of the closing stress, the height of the proppant pillar decreases gradually and the diameter increases.The pillar is compressed into a pancake shape.The compression from the shale plate changes the stress state in the proppant pillar and crushes down the pillar.Some particles near the edge of the pillar falls down and rolls away from the pillar.3.2.Fracturing Fluid Flowback Period.The fluid flowback calculation with coupled DEM-CFD analysis was carried out after the fracture closing.Figure11shows the profiles of the proppant pillar at different flow times.The initial height and diameter of the proppant pillar are 6 mm and 10 mm, respectively.The fluid viscosity is 0.001 Pa•s and the pressure difference is 64 kPa.The black arrows in the plot represent the fluid velocity vector at each fluid cell, and the velocity magnitude is indicated by the length of the arrow.The particles for rock plates are not shown in the plot.With the increase of flow time, the proppant pillar starts to spread and particles around the proppant pillar diffuse outward

Figure 7 :
Figure 7: Comparison of numerical modeling results and experimental results for pillar height at different closing stress.

3 Figure 8 :
Figure 8: Comparison of numerical modeling results and experimental results for pillar diameter at different closing stresses.

Figure 11 : 10 Figure 12 :)Figure 13 :Figure 14 :
Figure 11: Profiles of the proppant pillar at different flow times: (a) right after fracture closing, (b) flow time of 0.00015 s, (c) flow time of 0.0005 s, and (d) flow time of 0.0105 s.The red square in (b) represents the location of a monitoring fluid cell.

2 )Figure 16 :Figure 17 :
Figure 16: Evolution of the number of flowback proppant particles and proppant spreading area for different fluid viscosity.

Table 1 :
Summary of microscale properties of rock.