Hydraulic Fracture Conductivity in Shale Reservoirs

Optimum conductivity is essential for hydraulic fracturing due to its significant role in maintaining productivity. Hydraulic fracture networks with required fracture conductivities are decisive for the cost-effective production from unconventional shale reservoirs. Fracture conductivity reduces significantly in shale formations due to the high embedment of proppants. In this research, the mechanical properties of shale samples from Sungai Perlis beds, Terengganu, Malaysia, have been used for computational contact analysis of proppant between fracture surfaces. The finite element code in ANSYS is used to simulate the formation/ proppant contact-impact behavior in the fracture surface. In the numerical analysis, a material property of proppant and formation characteristics is introduced based on experimental investigation. The influences of formation load and resulted deformation of formation are calculated by total penetration of proppant. It has been found that the formation stresses on both sides of fractured result in high penetration of proppant in the fracture surfaces, although proppant remains un-deformed.


Introduction
The purpose of injecting proppants in shale reservoirs is to maintain the fracture conductivity for a longer period and to prevent the fracture from closure due to subsurface stresses.On the other hand, the proppants themselves can be a problem in the case where they develop surface penetration in the formation.As a result, the proppant is embedded into the formation and decreases the fracture conductivity of the reservoir as shown in Figure 1.Due to inhomogeneous stress distributions between quartz grains and proppants, high tensile stress concentrations beneath the area of contact between quartz grains and proppants are observed even at small external stress applied to the rock-proppant system.These high-stress concentrations are responsible for the early onset of damage at the fracture face and determine the type of proppant failure [1].
Water imbibition and some other tests on saturated shale were carried out to observe the crack generation process and compare the failure patterns as well as damage resistance of saturated shale kernels and unsaturated shale kernels.The average damage resistance of saturated kernel water is found to be around 11.69 MPa compared to 30.57MPa of unsaturated shale kernels, which implies that water can decrease the resistance to shale damage and helps in generating fractures [3].Fracture networks created during the process of hydraulic fracturing usually have a complex pattern.Most of these fractures are kept open by the incorporation of proppants in the form of proppant packs, as shown in Figure 1.
In the case of secondary fractures, other than bi-wing fractures, proppants are unable to enter into the fractured surface due to narrow apertures and thus, these fractures cannot maintain conductivity for a longer period.The effective vertical and horizontal stresses are responsible for the decrease in hydraulic fracture conductivity and an estimated 60% decrease in propped fracture conductivity occurs by increasing effective stress from 6.2 to 34.48 MPa [4].
Considering the narrow apertures of secondary hydraulic fractures, a partial monolayer of proppant, that is, a single layer of proppant having uneven distribution of proppants over the fracture surface, can be introduced instead of multilayer proppant to maintain the maximum possible conductivity for production improvement [5].The variation in the aperture and surface roughness of the hydraulic fractures are considered as main reasons for the uneven distribution of proppants.In this regard, the study of the conductivity of fractures with narrow apertures, filled with a monolayer of proppant, can be used for the optimization of hydraulic fracturing and the analysis of production in shale reservoirs.In the past, various types of compression, such as long-term and short-term compression on a single proppant, have been studied in-depth by diametric compression tests and DEM/ FEM simulations.Most proppants have shown creep behavior under long-term compression [6]. Figure 2 shows that the embedment potential is related to many factors especially the proppant material, shape, concentration, and ability of the proppant to resist sinking in the fracture zone [9].During hydraulic fracturing treatment, high fluid velocities in the fracture are generated by the small contact area between the wellbore and fracture, which results in the erosion of the proppant and fracture connectivity [10].
A computational fluid dynamics study with Eulerian granular modeling (EGM) that is based on solid pressure model and kinetic theory indicates that the transport of the proppant in complex fracture geometries is significantly affected by the dynamics of the fracturing fluids and the properties of the proppant [11].According to parametric studies, a higher injection rate and lightweight proppants are beneficial for the transport of the proppant through the fracture junctions and to carry proppant in hydraulic fractures and natural fractures [11].A DEM-CFD (discrete element method-computational fluid dynamics) and the experimental study indicate that during the closure period, the height of the proppants pillar decreases and diameter increases [12].The proppant flowback could occur easily with a large proppant pillar height or a large fluid pressure gradient.However, the higher bonding strength of the fibers results to improve the stability of the proppant pillar [12].Proppant pillar is defined as concentrations of proppant in the form of pillars that maintain the aperture of the hydraulic fractures.The change in the optimum distance which is defined as the distance between proppant packs that has the potential to maintain the maximum conductivity after proppant embedment under a sparse distribution condition is primarily controlled by closure pressure, the rock's elastic modulus, and the proppant elastic modulus.It also states that the proppant concentrations and the poroelastic effect do not influence this optimum distance [13].Studies based on analytical and discrete element method (DEM) have led to the understanding of the effects of various factors such as proppant size combination, concentration, time ratio, elastic modulus-to-stress ratio, and looseness coefficient [14,15].In these studies, deformation was considered elastic; however, actual phenomena can be captured by considering the intermediate states of elasticity and plasticity such as elastoplastic behavior of rock as well as proppant.In the case of monolayer proppant distribution, the embedment depth and contact stress decrease with the increase in proppant concentration [16].In the past, machine learning and computational fluid dynamics approaches have been used to explore the well operation and the transport of sand particles by the injection of foam [17][18][19][20][21][22][23].
The production performance of fractured wells depends on two factors, that is, formation parameters and fracture parameters [24].Formation parameters include porosity, permeability, and geo-mechanical properties of the formation, while the fracture parameters comprise a length, aperture, and conductivity of fractures [25,26].Hydraulic fracture conductivity reflects the transport capacity of the permeable channel through the reservoir and any alteration to this permeable channel will directly impact the stimulation achieved from the fracturing treatment [27].The experimental study performed on shale samples with fluids shows that the reduction in the elastic modulus can lead to a significant reduction in the effective fracture conductivity [28].Zhang et al. reported an 88% reduction in fracture conductivity by injecting water at 27.58 MPa closure stress [29].In this study, water as a fracturing fluid has been injected to find the excessive proppant embedment caused by the interaction of water with shale matrix, altering the hydraulic fracture conductivity.Water injection increased of local pore pressure and reduction of bonding strength of mineral in clay-rich shale that led to the softening of shale.The effects of rock stiffness, the roundness of proppant, and the effective stresses on the conductivity of the fracture were studied by a geomechanics-fluid mechanics-coupled numerical workflow considering the interaction between rock matrix and proppant as well as fluid flow in a hydraulic fracture during the process of the reservoir depletion [30].Compared to the weak shale, less embedment of proppant is observed in sandstone, having high stiffness, which indicates that the rock matrix with higher stiffness is helpful in maintaining the fracture aperture and conductivity [30].The correlation between the fracture conductivity and the corresponding production performance was quantitatively analyzed using the finite element method [31].The proposed research can provide valuable information on the unconventional maximization of resource recovery [31].For future extensions, a network of fractal fractures with a stochastic-based fractal fracture network combined with micro-seismic events can be coupled to quantify the complex fractures of the network to improve fracture conductivity and production performance [31].
The hydraulic fracturing process is a costly job; therefore, improving the reservoir quality evaluation mechanism and optimization of the technical parameters are important.As the low-quality hydrocarbon (HC) shale reservoirs are also gaining increasing attention, to this end, optimization of the hydraulic fracture conductivity is of utmost importance to make the job profitable.To estimate the actual hydraulic fracture conductivity in shale reservoirs, the computational contact analysis of proppant between fracture surfaces has been carried out in this study.In the numerical analysis, a material property of proppant and formation characteristics is introduced from the experimental analysis.The influence of formation load and resulting deformation of formation is calculated by the total penetration of proppant.The deformation mechanism and proppant embedment in shale rocks, saturated with fracking fluid, are then simulated.The finite element code in ANSYS is used to simulate the shale reservoir/proppant contact-impact behavior in the fracture surface.The embedment depth of the shale samples was obtained by numerical as well as experimental methods and the permeability was calculated by the Kozeny-Carman correlation.

Procedure
In this section, steps to measure embedment and fracture conductivity of fractured shale have been presented.Governing equations of numerical simulation to study the impact of embedment on the reduction of fracture conductivity have been presented in Section 2.3.The conductivity reduction due to embedment was modeled with a computational fluid dynamics (CFD) approach.We used the CFD software package CFX (ANSYS Inc.) in this work and applied the boundary conditions as shown in Figure 3 and Table 1.Initially, geotechnical characteristics of shale formation were calculated to define these properties in the software.

Elastic geo-mechanical properties
Dynamic elastic properties of the shale lithofacies, that is, Young's modulus and Poisson's ratio were calculated using compressional and shear velocities measured on shale core samples.The following equations [32] are used to obtain the respective values.

Dynamic Young Modulus
where V s and V p represent the shear and dynamic wave velocities in Km/s and P b is the bulk density of the shale in gm/cc.The results (Table 2) show that massive siliceous shale has a high value of Young's modulus and a low Poisson's ratio in Turb-KE which is turbulence kinetic energy, it is defined as half the sum of the variances (square of standard deviations) of the velocity components and Turb-Freq which is turbulent frequency here infers to the vortex shedding frequency (the frequency with which the vortices are shed behind the proppant during fluid flow).
Assumptions (1) The proppants are rigid body; (2) monolayer/single proppant is used in the embedment and conductivity analysis.In addition, deformation is calculated in the formation, while no deformation proppant only penetration in the formation is considered; and (3) the entire simulation is isothermal.

Table 1.
Basic input parameters, conditions, and assumptions.
comparison with massive argillaceous shale facie.Static Young's modulus and static Poisson's ratio were measured on cylindrical core samples of 2.0 in diameter and 4.0 in length.The specimens were failed in a triaxial setup and the deformations (axial and radial strains) were measured on the installed strain gauges.The static elastic parameters were calculated using the following equations.

Static Young Modulus
where σ a represents axial stress, while ε a and ε r represent axial and radial strains measured on the samples under deformation.Like dynamic parameters, the static Young's modulus of the massive siliceous shale is also greater than the static Young's modulus of argillaceous shale facie.Overall, values obtained for the static moduli are less than the dynamic moduli values measured on the same samples (Table 2).This is per previous findings done on producing shale formations.The correlation between dynamic and static moduli is shown in Figure 4. Strength parameters of the lithofacies are calculated using the equations below.

Calculation of embedment depth
The factors responsible for the change in fracture width and conductivity after hydraulic fracturing in shale reservoirs include proppant embedment and proppant deformation.The embedment of the proppant involves the penetration of the proppant inside the fracture surface, while proppant deformation is directly related to the strength of the proppant [33].The deformation of the proppant can be described by Eq. ( 7), as the proppant can be assumed as an elastic body while the penetration of proppant into a body can be solved using contact mechanics [34].The contact problem can be formulated as a constrained minimization problem, where the objective function to be minimized is the total potential energy (Π) of the bodies in contact.The energy for this system can be written as where k is the stiffness matrix, u is the displacement field, and f is the external force.Several constrained minimization algorithms can be used to solve the problem of the equation, such as the penalty method, the Lagrange multipliers method, and the augmented Lagrangian method.The results presented in this paper are based on the augmented Lagrangian method according to the ANSYS implementation.Augmented Lagrangian methods are a certain class of algorithms for solving c onstrained optimization problems.They have similarities to penalty methods in that they replace a constrained optimization problem with a series of unconstrained problems and add a penalty term to the objective; the difference is that the augmented Lagrangian method adds yet another term, designed to mimic a Lagrange multiplier.The augmented Lagrangian is related to, but not identical with the meth od of Lagrange multipliers.A general discussion of these techniques can be found in the literature on contact mechanics [35][36][37].In this study, a numerical model is developed based on the experimental investigation carried out on the embedment of 20/40 mesh proppant (size between 20 and 40 μm) on shale samples from Sungai Perlis beds, Terengganu, Malaysia.The core samples were subjected to uni-axial stress of 20 MPa to find the proppant embedment in the formation, as shown in Figure 4.Such a high compression force was applied to investigate the embedment under reservoir conditions.
The embedment cell consists of a transparent cylindrical tube where a shale sample is placed.A metal loading ram is used to load the shale-proppant stack and deformation is measured as the axial load is increased.The deformation in shale, assuming elastic behavior, is quantified using Young's modulus and the applied load.Figure 5 shows the universal testing machine (UTM) used for measuring  In this study, the contact behavior of proppant and rock was carried out using structural mechanics.Experimental boundary conditions are implemented to find the impact of load on proppant penetration in the rock surface.The resulting properties from experiments are introduced in a finite element model to find the effect of fracture surface on the embedment of proppant, as shown in Figure 6.A finite element study with ANSYS Workbench has been performed for the computational contact analysis [38].In a subsurface reservoir, proppants experiences compression from both sides in the formation; therefore, biaxial test should be carried for precise estimation of embedment.Proppants inside the fracture also contact each other due to subsurface stresses.Due to external force, the stressdisplacement relationship is as follows: where ρ is the density, kg/m 3 ; u is the displacement, mm; σ is the stress tensor; F v is the external force per unit volume; n = 1, indicating the proppant, n = 2, indicating the rock matrix.The stress-strain relationship and strain-displacement relationship are shown in Eqs. ( 9) and (10): where ε is the strain tensor; E is Young's modulus, GPa; υ is the Poisson's ratio.As we assume that Young's modulus and Poisson's ratio of the reservoir rock would change from the outer surface of the fracture toward the inside as shown in Figure 7 and the specific correlation is expressed as follows: where l is the distance from the surface of the fracture toward the inside as shown in Figure 7.The closure pressure causes the proppant to embed into the fracture surface and the porosity of the propped fracture to change.The embedment of the proppant (h em ) can be expressed as follows: where u m is the displacement of the no-contact part between the fracture and the proppant under closure pressure; u p is the displacement of the contact part between the fracture and the proppant.As the focus is proppant embedment into the shale formation, therefore, proppant is assumed as a rigid body.Thus, no deformation occurs in the proppant.However, the deformation of the shale surface is simulated with the penetration of proppant under uni-axial and bi-axial loads.

Calculation of fracture conductivity
CFD-CFX is used to calculate pressure drop and velocity across the inlets as shown in Figure 3.The finite volume method is adopted to solve the threedimensional Navier-Stokes equations.Consistent with the experimental conditions for conductivity measurements, the flow was the steady state at 25°C.The continuity and momentum balance equations for the steady-state flow are shown below.
Continuity equation where x, y, z are the dimensions; u, v, w are the velocity directions in x, y, z directions; p is pressure inside fracture.Fracture permeability was determined according to Darcy's law provided in Eq. (18).
where K is the permeability, Q is the flow rate of the injected fracturing fluid, μ is the viscosity, L is the length of the fracture around the proppant, A is the crosssectional area of the fracture zone, and ΔP is the differential pressure between inlet and outlet across the proppant.The conductivity of hydraulic fracture is generally defined as the maximum ability of the fracture to transmit a reservoir fluid through it.The conductivity is measured in μm 2 .cmbased on the propped fracture width (cm) and permeability (μm 2 ).

Fracture conductivity
where Kf is the proppant permeability, μm 2 and Wf is the fracture width, cm.

Experimental measured embedment
The main purpose of this experiment is to find the depth of the proppant embedded into the shale rock; the resulted embedment depths were used to validate the numerical model.Figure 8 shows the images after the proppant embedment test with some embedment occurring on the shale core at the (a) bottom and (b) top of the shale core sample with the 20/40 mesh proppants, that is, proppant having a size between 20 and 40 μm.The embedment depths of proppant in shale formation soaked underwater are recorded in the experiments, which are 76 μm at the top and 64 μm at the bottom surface.Table 3 shows the proppant embedment depth in soaked and unsoaked shale samples.

Numerical measured embedment
Contact analysis with finite element method has been performed to quantify the conductivity loss due to proppant embedment based on computational contact mechanics.A similar contact analysis has been also investigated between rail/wheel [26].Initially, the model was developed and validated according to the boundary load as applied in the experimental procedure.As shown in Figure 9, the load is  applied at the top surface, while the lower surface has been kept static to validate the results with a uni-axial compression test.In the case of the uniaxial load test, the resulted embedment depth in this study is found different for the top and bottom of the fracture surface, similar to an earlier experimental study performed using different types of proppants [39].The significant difference between embedment profiles is the result of different proppant types being tested.20/40 Ceramic proppants are rounder, more spherical, uniform, and stronger than 20/40 Ottawa proppants.This difference in proppants makes 20/40 Ottawa sand more prone to proppant embedment as well as any other damage mechanism caused by mechanical and chemical factors in fractures [39].The present study shows that even though the proppant is strong enough to get deformed, and there is a penetration of proppant in the rock surface due to fracking fluid flow that has reduced Young's modulus of the fractured surface.Initially, ceramic proppant and rock properties are introduced based on the shale formation.A pair of contact between surfaces and proppant is defined.Then, the external loads are applied to the rock surfaces to obtain the stress transfer across the proppant of contact surfaces.The penalty-based method is used to simulate the contact behavior [27,40].The finite element method is a numerical method that can be successfully used to generate solutions for problems belonging to a vast array of engineering fields: stationary, transitory, linear, or nonlinear problems.For the linear case, computing the solution to the given problem is a straightforward process, and the displacements are obtained in a single step and all the other quantities are evaluated afterward.When faced with a nonlinear problem, in this case with a contact nonlinearity, one needs to account for the fact that the stiffness matrix of the systems varies with the loading, the force vs. stiffness relation being unknown before the beginning of the analysis.Modern software using the finite element method to solve contact problems usually approaches such problems via two basic theories that, although different in their approaches, lead to the desired solutions.One of the theories is known as the penalty function method [40].The penalty method is simple to implement in practice.The penalty is a sort of friction coefficient, and one can specify a friction model that defines the force resisting the relative tangential motion of the surfaces in a mechanical contact analysis.By selecting a penalty, one can use a stiffness (penalty) method that permits some relative motion of the surfaces when they should be sticking.By applying the penalty method, the penetration of the proppant has been achieved higher at the top which is 75 μm, whereas the penetration at the bottom surface is recorded 60 μm.The penetration of proppant in numerical model has been achieved almost the same as recorded in the experiments, which was 76 μm at the top and 64 μm at the bottom surface.The similarity of the results shows that the developed model has satisfactory results and parametric study can be carried out for further analysis.Once the model has been validated with the experimental results, then the external force was applied on both sides of the proppant to represent the actual condition in the fracture formation.

Impact of embedment on fracture conductivity
In this section, the change in pressure and velocity of backflow of fluid across the proppant has been presented.Embedment has a profound impact on the pressure drop as well as velocity profile as shown in Figure 10.In this study, three different embedment cases have been considered (0, 60, and 80%).The percentage of embedment is defined as the proportion of the total proppant that is embedded through the fracture surface.Without embedment, a slight difference between inlet and outlet pressure has been recorded; however, a significant difference between inlet and outlet pressures can be seen at 60 and 80% embedment as shown in Figure 10(b) and (c).Inlet velocity in all cases is 0.5 m/s but around the proppant, the flow velocity is recorded around 2 m/s and at outlet, the velocity is achieved 1.5 m/s.
Based on the different embedment depths, the velocity of the injected fluid varies significantly as shown in Figure 11.In all cases, injection velocity is constant, that is, 0.5 m/s.A sudden increase in the fluid velocity is recorded around the proppant and a decrease in velocities is presented at the end of the proppant.The results show a significant decline in the velocity at 80% embedment; therefore, fracture conductivity is recorded significantly low at high embedment (see Figure 12).As fluid flowing continues around the exit sides of the proppant, it begins to slow down due to eddy generated at the outlet/backside of the proppant.The Lagrangian analysis is capable of revealing the underlying structure and complex phenomena in unsteady flows [41].
Distance is measured from one side (fracture inlet) of the proppant until the other side (fracture outlet) of the proppant.All the three positions of proppants have been presented in Figure 10(a), (c), and (e).Numerical analysis is conducted by finite volume method to obtain pressure drop across the proppant and resultant fracture permeability.Table 4 shows the fracture conductivity based on embedment percentage.
Finally, fracture conductivity is achieved based on fracture permeability and fracture width.Figure 12 shows that fracture conductivity was measured based on  embedment depth obtained with experimental and finite element methods.A dramatic decrease in fracture conductivity has been obtained with the increase of embedment depth.The reason for a significant decrease in the conductivity is the significant pressure drop across the embedment.The results show that pressure loss at 60 and 80% embedment is 29,716 and 64,721 pa, respectively.

Conclusion
A numerical model is developed for contact analysis of proppant embedment in the formation based on experimental investigation.Initially, the model was developed based on the experimental design of proppant embedment in the laboratory where the load is applied uni-axially from the top.Then, the study is extended by applying load from top as well as from the bottom side of the proppant in the fractured surface to simulate the actual reservoir condition.The amount of proppant embedment has been computed on both sides of the proppant in the fracture surface.Also, the deformation and normal stress profile have been plotted along with the formation and proppant.The total penetration of the surfaces has been recorded 141 μm on each side as the equal loads have been applied on both sides of formations around the proppant.This shows that actual proppant embedment is very high if stresses are present on both sides of the proppant in the fracture.The computational contact mechanics analyses have been able to capture the actual conductivity of fracture showing that the finite element method can be used to estimate embedment depth and has comparable results with experimental measurements.Long-term production of hydrocarbon from shale reservoirs is directly related to fracture conductivity in the hydraulically stimulated reservoir volume.This study shows that the uncertainty and reduction in hydrocarbon production profile with time can be mimicked by exact estimation of proppant embedment and fracture closure with finite element method since it relates to fracture conductivity.The presented method can serve as a valuable criterion to effectively reduce the loss of hydraulic fracture conductivity in shale reservoirs with time.Based on this numerical model, the required fracture conductivity can be achieved by keeping the extra width of fracture in the design criteria to reduce the conductivity loss in the formation.

Figure 3 .
Figure 3.The direction of flow at inlet and outlet across the proppant.

Figure 4 .
Figure 4. Reduction in the fracture conductivity by embedment in the fracture surfaces, modified from Zhang et al. [29].

Figure 5 .
Figure 5. Embedment by compression (a) embedment test by bi-axial compression and (b) sample under UTM machine.
compression and tensile strengths of materials and samples used in the proppant embedment tests.

Figure 6 .
Figure 6.Illustration of embedment of the surface due to external load.

Figure 7 .
Figure 7. Change in Young's modulus due to fracturing fluid interaction with the shale.

Figure 8 .
Figure 8. SEM images of proppant (20/40 mesh) embedment on the original core soaked in fracking fluid: (a) bottom of the core and (b) top of the core.

Figure 9 .
Figure 9. Proppant under vertical principal stress: (a) total deformation and penetration due to an external load, (b) directional deformation along the Z axis, one surface under load while another surface in static, (c) plot of directional deformation profile, (d) normal stress profile along with the formation and proppant, and (e) plot of normal stress profile along with the formation and proppant under uni-axial load.

Figure 10 .
Figure 10.Velocity and pressure profiles in the fracture zone around the single proppant with and without embedment.(a) Pressure profile with no embedment.(b) Velocity streamlines with no embedment.(c) Pressure profile with 60% embedment.(d) Velocity streamlines (a streamline is a line that is tangential to the instantaneous velocity direction (velocity is a vector, and it has a magnitude and a direction.Color represents velocity magnitude) with 60% embedment).(e) Pressure profile with 80% embedment.(f) Velocity streamlines with 80% embedment.

Figure 11 .
Figure 11.Velocity profile of injection fluid around the proppant under different embedment percentages.

Figure 12 .
Figure 12.Fracture conductivity obtained with finite volume method based on experimental and numerical measured embedment depths with finite element method.

Table 2 .
Elastic and strength properties of the lithofacies identified in this study.
18] Khan JA, Pao WK.Fill removal with foam in horizontal well cleaning in coiled tubing.Res.J. Appl.Sci.Eng.Technol.2013, 6, 2655-2661 [19] Khan JA et al.Optimization of coiled tubing nozzle for sand removal from wellbore.Journal of Petroleum Exploration and Production Technology.2020;10(1):53-66 [20] Khan J, Pao WK.Horizontal well cleanup operation using foam in different coiled tubing/annulus diameter ratios.American Journal of Applied Sciences.2013;14:3235-3241 [21] Khan JA, Pao WK.Effect of different qualities of foam on fill particle transport in horizontal well cleanup operation using coiled tubing.In: Advanced Materials Research.Trans Tech Publications; 2014 [22] Khan JA et al.Comparison of machine learning classifiers for accurate prediction of real-time stuck pipe incidents.Energies.2020;13(14):3683 [23] Khan JA et al.Quantitative analysis of blowout preventer flat time for well control operation: Value added data aimed at performance enhancement.27] Cooke C Jr. Conductivity of fracture proppants in multiple layers.Journal of Petroleum Technology.1973;25(09): 1101-1107 [28] Akrad OM, Miskimins JL, Prasad M. The effects of fracturing fluids on shale rock mechanical properties and proppant embedment.In: SPE Annual Technical Conference and Exhibition, Denver, Colorado, USA.Society of Petroleum Engineers; 2011 [29] Zhang J et al.Experimental and numerical studies of reduced fracture conductivity due to proppant embedment in the shale reservoir.Journal of Petroleum Science and Engineering.2015;130:37-45 [30] Fan M et al.Investigating the impact of proppant embedment and compaction on fracture conductivity using a continuum mechanics, DEM, and LBM coupled approach.In: 52nd US Rock Mechanics/Geomechanics Symposium.Seattle, Washington USA: 35] Andrei N. Penalty and augmented Lagrangian methods.In: Continuous Nonlinear Optimization for Engineering Applications in GAMS Technology.Springer.;2017.pp.185-201 [[[