Numerical Modeling of the Conductivity of the Particle Monolayer with Reduced Size

Fractures filled with a proppant monolayer play an important role in the hydraulic fracture network. Predicting the conductivity of these fractures is the basis of fracture network optimization. However, little attention has been paid to the conductivity of the proppant monolayer. The change of conductivity under various conditions is currently not fully understood. Therefore, in this paper, the conductivity variation under different conditions are simulated. The reduction of particle size was calculated by existing analytical models. The permeability variation was calculated through computational fluid dynamics (CFD) combined with COMSOL Multiphysics. The controlling factors of conductivity under a proppant monolayer were identified. Simulation results indicate that elastic parameters, closure pressure, and proppant distribution have significant influence on conductivity, while creep parameters, such as rock viscosity and time, have limited influence on conductivity. Moreover, the changes in permeability, porosity, and tortuosity with variation of embedment were analyzed. Results indicated that with an increase in embedment, the permeability and porosity decrease as expected. The main reduction (nearly half) emerges in the first 20% of proppant embedment. Furthermore, the permeability of a single particle deviates largely from the prediction of Carman-Kozeny (CK) equation. The tortuosity of proppant particle increases with a decrease in particle size due to embedment. A modification of the Carman-Kozeny equation is proposed to address this influence.


Introduction
The reserve of oil and gas in conventional and unconventional resources is abundant [1,2].Unconventional oil and gas reservoirs have become the focus of exploration and exploitation in petroleum industry.Horizontal drilling and multistage fracturing are widely used in the exploitation of these reservoirs [3].Accurately measuring and understanding these technologies are critical for ensuring successful reservoir exploitation [4].For hydraulic fracturing, it consists of initiating, propagating, and maintaining a fracture from the wellbore to pay zones.Proppants are pumped into the created fracture along with the fracturing fluid.When the hydraulic pressure is removed from a well, the surfaces of fractures compress onto the proppants, creating a highpermeability pathway in deep rock formations through which natural gas and oil flow more freely.
The performance of fractured wells is controlled by two types of parameters.One type of parameters is formation parameters, such as porosity, permeability, and mechanical properties [5].The second type of parameters is fracture parameters, such as fracture length and fracture conductivity [6,7].Fracture conductivity reflects the transport capacity of the high permeable channel, which limits the amount of stimulation achieved from the fracturing treatments of many wells [8].It is defined as the product of fracture width and the apparent permeability of the proppant pack.Results summarized from over 80 field studies show the benefit of increasing fracture conductivity [9].As low-quality reservoirs become the increasingly popular focus, increasing the conductivity of hydraulic fracture is one of the possible strategies in the future to improve the fracturing response [10].Any measures that intended to increase the fracture conductivity will increase the cost.The International Energy Agency (IEA) indicates that the new worldwide oil and gas infrastructure is projected to cost US $20 trillion between 2011 and 2035 [11].However, no matter whether the project is onshore or offshore, many oil and gas projects experienced significant cost overruns [12].Under the condition of low oil and gas prices, improving the evaluation efficiency of oil and gas reservoir quality and optimizing the technical parameters are two effective ways to cut the cost.Therefore, it is worth optimizing the fracture conductivity in order to reduce the costs and improve postfrac performance.
Slick water fracturing is widely used in the development of tight reservoirs [13].A complex network of narrow secondary fractures is generated to extend the contact area between the network of flow channels and formation.Some of these narrow fractures are only propped by a proppant monolayer [14].Initially, the industry believes that a monolayer is difficult to achieve.However, the situation changes in unconventional reservoirs; many hydraulically activated shear fractures propagate deeply in the reservoir.These fractures could only open at a limited width and would only be poorly propped by a monolayer of proppants.Therefore, the conductivity of the proppant monolayer is an important factor in the design and optimization of hydraulic fracturing.
But little attention has been paid in previous studies.The dominant factors that affect the conductivity of the proppant monolayer are not clear.Besides, the prediction of fracture conductivity is the foundation of fracture optimization.In order to predict the fracture conductivity with the proppant monolayer, a good understanding of the changes of fracture width and permeability of the proppant pack is required.Proppant embedment is the main cause of fracture width reduction, which can be solved by the models available in the literatures.On the calculation of permeability of proppant pack, Carman-Kozeny (CK) equation was widely used [15][16][17][18].The semiempirical Carman-Kozeny (CK) equation is the most famous model that describes the permeability-porosity relation.However, this relation has many limitations.It assumes that the porous medium is a packing of identical spheres and the tortuosity is an empirical constant.Then, when there is only a proppant monolayer, whether the equation is still applicable is unclear.Furthermore, after proppant embedding into rock, the geometry of the proppant is no longer a circular shape, as shown in Figure 1.How will this affect the permeability and its deviation from the prediction of CK equation?Therefore, in this paper, the conductivity of the proppant monolayer was studied.First, computational fluid dynamic theory is used to study the influencing factors controlling conductivity under a proppant monolayer.Then, the permeability of the proppant monolayer with different levels of embedment was studied.A modification of CK equation is proposed based on numerical results.The models and results can provide an insight into the features and factors controlling the conductivity of the proppant monolayer.

Numerical Methodology
A computational fluid dynamics (CFD) approach was utilized to model the conductivity reduction due to proppant embedment.A series of numerical simulation were performed to study the variation of fracture conductivity with the proppant monolayer.CFD has been widely used in studying fluid flow [14,17,19,20].In this work, COMSOL Multiphysics was used.If the structure of the medium is determined, CFD can solve the steady-state Navier-Stokes equations to obtain the velocity field and then compute the permeability.Incompressible fluid and laminar flow model were used, as shown in Table 1.
Actually, CFD simulation of fluid flow inside a field scale geometry will require large computational efforts, which exceeds the capacity of most computers.In order to improve the computational efficiency, three computational domains were used as the physical model, as shown in Figure 2. One is a simple regular lattice called facecentered cubic (FCC) packing.Fluid flow in a FCC packing was demonstrated to be consistent with the Carman-Kozeny equation [21].Therefore, it was used to verify the model and numerical solution of COMSOL Multiphysics.The second one is a single particle, which is the focus of this paper.The third one includes multiple particles,  2 Geofluids which are used to study the influence of proppant distribution on conductivity.The Carman-Kozeny equation is a well-known formula and a relatively simple way to predict permeability.The permeability k can then be expressed as where k is the permeability, m 2 ; ϕ is the porosity, dimensionless; τ is the tortuosity, dimensionless; and s is the specific surface area per pore volume, 1/m.If the porous medium is a packing of identical spheres, then s is given by where d is the diameter of particles, m. s in equation ( 2) does not consider the influence of embedment.In later calculation, the definition of s is used to calculate its value.Equation ( 2) can be substituted into equation (1) to yield For random packing of spheres with the same radius, τ 2 = 2 5 is commonly used and results in the following well-known equation: For the numerical simulations, it is assumed that the flow was three-dimensional, single-phase, and laminar.Each domain also contains an entrance and an exit on each side.A constant pressure boundary was used at the inlet.A constant zero static pressure was used at the outlet.Wall boundaries were set on the other surfaces.Grid size affects the accuracy of the numerical simulation.We conducted a series of calculations with different grid sizes to check the computational stability and accuracy.Each set of comparison processes (data in the same figure) used the same mesh to ensure the convergence and computational speed.Simulation results indicate that the numerical permeability obtained by this mesh is 919 μm 2 , the analytical permeability calculated by CK equation is 928 μm 2 , and the difference is 0.9%.

Numerical Results and Discussions
In the reference [15], analytical models were derived to calculate proppant embedment, proppant deformation, change in fracture aperture, and fracture conductivity.And then, this model has been improved by Zhang and Hou [16] with consideration of the rock creep deformation.Here, we would like to give a brief introduction of their models.
Proppant embedment, proppant deformation, and the change of fracture width were derived as follows: Then, considering creep deformation, these parameters change to  3 Geofluids Poisson's ratio of rock, mm; E 1 is the elastic modulus of the proppant, MPa; E 2 is the elastic modulus of rock, MPa; h is the value of embedment, mm; α is the change in fracture aperture, mm; η 1 is the proppant viscosity, mPa•s; and η 2 is the rock viscosity, mPa•s.
This paper focuses on the effect of different factors on the conductivity of the proppant monolayer.Therefore, analytical models [15,16] were used to calculate the embedment depth and fracture residual width under different conditions.Then, fracture conductivity was evaluated by a CFD method.

Effect of Elastic Parameters.
To investigate the effect of elastic parameters on fracture conductivity, a normalized parameter can be defined as follows: where E is Young's modulus (1 is proppant; 2 is rock), Pa.Simulation tests were carried out with six different closure pressure values under different E n .
Figure 3 shows the relationship between E n and normalized conductivity at different stress conditions.The conductivity after proppant embedment and deformation is normalized by the conductivity of fracture with a width that equals the proppant diameter (0.85 mm).As shown in the figure, the conductivity decreases when E n increases for all the cases.With an increase in the proppant elastic modulus or a decrease in the rock elastic modulus, the fracture conductivity shows an approximately linear decrease.When the closure pressure is 10 MPa, the conductivity decrease 6.66% with an increase in E n .However, when the closure is 60 MPa, the conductivity decrease 26.91% with an increase in E n .Under high closure pressure, the conductivity is more sensitive to E n .For some reservoirs, such as shale, the rock elastic modulus will decline when encountered with fracturing fluid [22].If the proppant with high elastic modulus is used, the E n will become higher too, which would undoubtedly result in sharp decrease in fracture conductivity.4 and 5 show the relationship among rock viscosity, closure pressure, and normalized conductivity.As the rock viscosity increases, the conductivity increases.It is found that the conductivity increases relatively obvious when the rock viscosity increases from 1000 mPa•s to 4000 mPa•s and gradually levels off when the rock viscosity is above 4000 mPa•s.

Geofluids
As shown in the figure, the conductivity increases by 4.22% as the rock viscosity increases from 1000 mPa•s to 4000 mPa•s under 10 MPa closure pressure and by 16.62% when the closure pressure is 60 MPa.In the case of low closure pressure, the increase of rock viscosity does not make a big difference to the increment of conductivity.
However, the influence of closure pressure on conductivity is significant.As shown in Figures 4 and 5, the conductivity decreases by 44.1% as the closure pressure increases from 10 MPa to 60 MPa with 1000 mPa•s rock viscosity and by 35.78% when the rock viscosity is 10000 mPa•s.
Besides, as shown in Figures 6 and 7, the conductivity decreases by 11.93% as the time increases from 100 d to 600 d at rock viscosity of 1000 mPa•s and by 8.07% when the rock viscosity is 10000 mPa•s.The time does not have a significant influence on conductivity.

Effect of Diameter and Distribution of the Proppant.
Figures 8 and 9 show the relationship among proppant diameter, closure pressure, and normalized conductivity.As the proppant diameter increases, the conductivity increases sharply.The relationship between the particle size and conductivity satisfies the power function.
As shown in Figures 8 and 9, the conductivity increases by 49.00% as the closure pressure increases from 10 MPa to 60 MPa with a proppant diameter of 0.2 mm and by 28.87% when the proppant diameter is 0.85 mm.  5 Geofluids Several proppant distribution cases were simulated as shown in Figure 10.When considering proppant embedment, the changes of physical models are illustrated in Figure 11.Simulation results are shown in Figure 12.Embedment ratio is defined as the fraction of the proppant grain embedment to the initial fracture width (proppant diameter).As the proppant distribution changes, the conductivity varies significantly.The patterns of A and E are better than the other ones.The results demonstrate that the discontinuous proppant placement is better than continuous placement under different embedment ratios.

Modification of Carman-Kozeny Equation.
Carman-Kozeny equation can be used to calculate the conductivity as evidenced in many analytical models, such as reference [15].However, as described before, the Carman-Kozeny    6 Geofluids equation does not consider the particle size change resulted from embedment.A modification of the CK equation is proposed based on numerical results.Three particle sizes were simulated, ranging from 212 μm to 840 μm (70 mesh to 20 mesh).The results are shown in Figure 13.The permeability is normalized against the permeability of a channel with the maximum initial fracture width, which is the initial permeability of 20-mesh proppantsupported fracture.The permeability can be affected by the proppant size and the embedment ratio.With an increase in proppant size and a decrease in the embedment ratio, the permeability increases.For the proppant monolayer, the permeability reduces to its half when the embedment ratio increases to nearly 20%.Proppant has two embedment mechanisms, elastic embedment and creep embedment [16].The fracture will keep closing under creep deformation.
Furthermore, results show that with an increase in the embedment ratio, the porosity decreases and the specific surface area remains constant, as shown in Figure 14.When the porosity decreases to a certain value, the decreased amplitude reduces.After the proppant embedment increases to 90%, the porosity reduces to 20%.From equation (1), the decreased porosity will also lead to a decrease in permeability.However, the prediction of the CK equation is still larger than the numerical results.
Therefore, tortuosity induces the difference.Previous researches showed that tortuosity will change with porosity [23].Under each different situation, they meet a specific relationship.Four equations are used to fit the relationship between porosity and tortuosity.Their results are shown in Figure 15.
Results are summarized in Table 2.It can be seen that all the equations' calculations fit the numerical results well with a R 2 > 0 9. Equation (11) achieves the best fitting with a R 2 = 1.However, when these results were used in the permeability calculation, only equation (11) matches well with the results, as shown in Figure 16.

Conclusions
We have conducted a numerical study on the fracture conductivity change caused by proppant embedment.Simulation results enhance the understanding of the conductivity variation under proppant monolayer condition.The following conclusions can be summarized according to the present study: (1) Simulation results indicate that elastic parameters, closure pressure, and proppant distribution have significant influence on conductivity, while creep parameters, such as rock viscosity and time, have limited influence on conductivity (2) The numerical results showed that the first 20% of proppant embedment can cause a significant permeability reduction (3) For a single proppant particle, with an increase in the embedment ratio, the porosity decreases and the 7 Geofluids specific surface area keeps constant.After the proppant embedment increases to 90%, there is still 20% porosity left.Meanwhile, tortuosity of the proppant monolayer increases with the increase in embedment (4) The numerical simulation was able to capture the flow characteristics around the proppant.The permeability of a single particle deviates largely from the prediction of the CK equation.A new correlation between porosity and tortuosity was proposed to solve this deviation (5) The placement of the proppant is very complex, which brings great difficulty to the prediction of the fracture conductivity.One should note that in a real fracture there is no perfect placement with uniform size throughout a complex fracture network.The models in this paper focus on the base case of the proppant monolayer.Any displacement of the     8) Equation ( 9) Equation ( 10) Equation ( 11) Numerical results Normalized permeability (numerical results) Normalized permeability (( 8), ( 9), (10), and ( 11 8), ( 9), (10), and (11).
8 Geofluids proppant monolayer could be a combination of single particle cases.However, the models consider single phase, laminar flow, single particle, and their combination.Any deviation from the assumptions can make the prediction of the models different from real cases.

Figure 2 :
Figure 2: Schematic of the computational domains.

Figure 3 :Figure 4 :Figure 5 :
Figure 3: Effect of E n on fracture conductivity at different closure pressures.

Figure 6 :Figure 7 :Figure 8 :
Figure 6: Effect of rock viscosity on conductivity at different times.

Figure 9 :
Figure 9: Effect of closure pressure on conductivity at different proppant diameters.

Figure 11 :
Figure 11: Illustration of different embedment ratios for distribution (a) in Figure 10.

Figure 12 :Figure 13 :
Figure 12: Effect of proppant distribution on conductivity at different embedment ratios.

Figure 14 :
Figure 14: Effect of embedment on porosity and specific surface area.

Table 1 :
Basic parameters for flow simulation.
ModelFluid type Viscosity (mPa•s) Density (kg/m 3 )Laminar flow Water 1 1000 η 2 /E 2 , 6where β is the deformation of the proppant, mm; D 1 is the diameter of the proppant, mm; D 2 is the thickness of rock, mm; K is the distance coefficient; p is the closure pressure, MPa; v 1 is Poisson's ratio of the proppant, mm; v 2 is

Table 2 :
Fitting results for different equations.
It is recommended to develop a more comprehensive model by the method given in this paper, such as multilayers and two-phase flow Diameter of particles, m D 1 : Diameter of the proppant, mm D 2 : Thickness of rock, mm e: Constant E 1 : Elastic modulus of the proppant, MPa E 2 : Elastic modulus of rock, MPa E n : Dimensionless Young's modulus h: Value of embedment, mm k: Permeability, m 2 K: Distance coefficient p: Closure pressure, MPa s: Specific surface area, 1/m τ: Tortuosity, dimensionless ϕ: Porosity, dimensionless α: Change in fracture aperture, mm β: Deformation of the proppant, v 1 : Poisson's ratio of the proppant v 2 : Poisson's ratio of rock η 1 : Proppant viscosity, mPa•s η 2 : Rock viscosity, mPa•s.