On how asymmetric stimulated rock volume in shales may impact casing integrity

Microseismic data and production logs in our study area have confirmed asymmetric developments of the stimulation rock volume in shales with respect to the wellbore, while severe casing deformation problems have been reported frequently in this area. Here, we propose a systematic methodology to investigate the possibility of casing failure due to strong shear stresses induced by asymmetric stimulated zones. A mechanical earth modeling (MEM) is initially performed to determine the in situ stress field in the target layer before fracturing by incorporating the existing geological features, logging data, and rock anisotropy. Then, we provide a computationally cheap and efficient estimation for stimulated rock volume of each stage by considering the possible overlaps in adjacent stages based on the clustered microseismic clouds. Using this approach, a reservoir‐scale 3D coupled model tied to a more detailed near‐wellbore part incorporating the casing string and the cement sheath is established to simulate the development of stimulation zones, stress redistribution, and their impacts on casing deformation as each stage fracturing treatment chronologically goes on. Our numerical results indicate that continuous redistribution and re‐orientation of stress field near the borehole are tracked during pumping the treatment which reveals formation of some pockets of tensile stresses along and around the wellbore. Asymmetric stimulations are observed to generate strong shear stress on the suspended casing. These shear forces result in deflection and S‐shape deformations accompanied with cross‐sectional ovality. Some regions receive repeating treatments, which results in intensifying formation stress heterogeneity and worsen casing deformation severity. The calculation results are compared with measurement of multi‐finger imaging tool (MIT) to validate the accuracy. Our analysis has indicated that simply increasing the flexural strength by increasing thickness of casing cannot radically mitigate casing deformation problems. This paper presents a novel workflow for a coupled modeling of casing deformation during hydraulic fracturing operations, while current modeling efforts assume symmetric bi‐wing fracture geometries.


| INTRODUCTION
Shale gas reservoirs with ultra-low permeability have progressively drawn the industry attention for exploration and production activities to meet the rapid increase of the global energy consumption. These reservoirs generally do not possess natural productivity, and conventional stimulation techniques cannot easily turn these reservoirs into highly productive ones. Hence, extensive multistage hydraulic fracturing stimulations in recent years have been a renewed technique to address these limitations by developing complex fracture networks in the rock to accelerate fluid production. [1][2][3] In recent years, there have been some efforts to improve ultimate recovery by enhancing fracture network complexity, while often well integrity has been compromised due to the focused attention on improving stimulation outcomes. [4][5][6][7] Consequently, in some areas, severe casing deformation problems have occurred during hydraulic fracturing. During the period from 2005 to 2015, over fifty thousand multi-fractured horizontal wells for unconventional resource plays have been performed in North America with approximately one million total fracture stage treatments. 8 However, the high-volume fracturing operations have also exacerbated casing failure problems, for instance in Texas, casing failures related to various aspects of multistage fracturing operations have caused 8 blowout incidents. 9 In China, casing failure problem is even worse, over 40 percent shale gas wells suffered from this problem to different extents, [10][11][12][13][14] which has caused downhole tools like milling shoe or gauge cutter to get stuck at the deformed parts of the casing when running into the borehole after the frack job. Even in the worst cases, the well had to be plugged and abandoned. Therefore, it is important to anticipate potential casing issues during multistage fracturing with minimum costs using simulation tools and revise the pumping schedule and stage design accordingly to avoid these problems.
The serious hazards of the excessive casing deformation to the subsequent production especially in shale gas reservoirs are already highlighted in the petroleum engineering literature. However, the driving mechanisms behind the casing failure have still remained unknown due to the complexity of the events involved in multistage hydraulic fracturing jobs. Most of the current studies about casing deformation have discussed various hypotheses such as fault slippage, 15,16 temperature change, 17 poor cement quality, 18 and perforation damage. 19 Daneshy 20 was one of the first references to document the impact of hydraulic fracturing on borehole stability and casing failure. He claimed that asymmetric fractures created by off-balance fracturing with respect to the borehole might result in the differential compaction of stimulation zones, which could develop strong shear stresses on casing and cause its failure along the weakest interface. However, this argument is remained to be verified by experiments or simulations. Furui et al 21 conducted a series of triaxial tests on selected reservoir rocks and found large localized inelastic deformations occurs in presence of high-confining stresses during acidization and fracturing operations. They found the local compaction of the reservoir could increase the axial-compression load on casing and result in its yielding or buckling. Sugden et al 22 discussed different factors such as aggressive trajectories, cement quality, high pumping rate, and temperature differences on augmenting the load that casing is experiencing during multistage fracturing of shale wells. Through simplified three-dimensional modeling of a subsurface formation for hydraulic fracturing and sub-model of a selected portion for casing failure, Shen et al 23 investigated the impact of high injection pressure on casing integrity during multistage fracturing. They simulated various factors including the value of injection pressure, distribution of natural fractures, cementing quality, and initial geostress, which may influence the upper bound of the injection pressure window (IPW) and discussed their impacts on casing integrity. From the near-wellbore perspective, Lin et al 24 introduced a voronoi tessellation modeling method to obtain a quantitative relationship between rock strength and induced fractures, based on which they illustrated how the complexity of fracture network may influence casing deformation during and after the hydraulic fracturing. Haghshenas et al 25 found that significant stress concentration developed in severe doglegs during multistage fracturing operations can cause casing connection failure and suggested to reduce localized dogleg severity to less than 8°/ft that could prevent casing failure, effectively. Yan et al 18 argued the self-weighted of casing itself and complex well trajectory in shales result in the casing touching the borehole wall on one side which leave some drilling mud pockets and lead to cement voids. The severe eccentricity of cement can exacerbate load nonuniformity on casing to make it exceed the yield stress under the combined effect of internal and external pressure. Actually, the impact of cement eccentricity on casing deformation could be minimal because the mechanical properties of shale and cement are very similar. From the perspective of fatigue cumulative damage, Barreda et al 26 conducted a case study in US onshore unconventional wells and analyzed the impact of cyclic loading during pumping on casing distortion specially in multistage fracturing operations. They speculated the magnitude of von Mises stress on casing fluctuates up and down continuously during fracturing cycles, which makes the casing prone to fatigue failure, although the number of cycles is much less than typical number of cycles required for fatigue appearance. By accounting for the thermal-seepage-solid coupling interaction and mechanical anisotropy of shale, Xi et al 13 developed several finite element models to discuss the potential impacts of fault slippage during hydraulic fracturing on casing failure. They conducted a series of sensitivity analyses to determine the weight of the main influence factors on casing failure including sliding distance, casing inner pressure, casing thickness, and cement quality. Their results showed that decreasing sliding distance and maintaining high inner pressure could help improve the loading conditions around casing. Increasing the thickness of casing and the elasticity of cement sheath could help enhance the resistance to the nonuniform load induced by fault slippage.
Despite several studies have been conducted about casing failure during hydraulic fracturing, current models only consider rock deformations very close to the wellbore. None of these models try to tie further changes occurred to the formation during multistage fracturing to assess casing deformations at the reservoir scale. However, modeling such large complex fracture networks induced in each stage of hydraulic fracturing and how the stimulation zones may overlap in adjacent stages introduce major challenges in terms of computational costs. Therefore, a crucial step is to locate fractured zones characterized with elevated pore pressure and deteriorated rock strengths due to hydraulic fracturing. To pursue this objective, an integrated methodology is proposed in this paper which consists of using the location of microseismic events to estimate dimensions of the stimulated zones. Following this method, a 3D coupled poroelastic model of the reservoir is established for fracture network growth and stress redistribution due to hydraulic fracturing to obtain a more realistic condition of the casing deformation. Eventually, transverse and axial deformations of casing during a multistage fracturing job are calculated through numerical simulations to investigate potential reasons behind casing failure, and finally, some suggestions to avoid such situations in the future are discussed.

| MODEL INITIALIZATION
Analysis of casing deformation during multistage fracturing requires knowledge of in situ stresses prior fracturing as the initial condition. Hence, we initially perform a mechanical earth modeling (MEM) to determine in situ stress in the target layer and along the well trajectory more accurately (see Figure 1).
To provide a meaningful description of our analysis and comparison of the outcomes with the field condition that stress differential is large and the stress distribution is nonuniform due to the existing geological features, we build a model for a well drilled in the Longmaxi shale of the Lower Silurian (S 1 l) located in Sichuan Basin, China. The main body of the geological structure is a major anticline called XDZ whose direction is NE-SW. To avoid the boundary effects on our stress analysis, upper layer-P 2 l and lower layer-O 3 W are also incorporated into the model. Figure 1C  XY-H1 is the analyzed well (marked as red star in Figure 1C) that suffered from severe casing deformations after fracking. Stress gradients of four wells (marked as circles in Figure  1C are measured through Kaiser-DRA method using samples from different depths. The memory effects in acoustic emission (Kaiser effect) and strains through deformation rate analysis (DRA) can be observed and recorded while the rock specimens are cyclically loaded, based on which the in situ stress state can be obtained when applied stress reaches the previous peak value. 27,28 Stress data of two wells (XY-H2 and XY-H7) are used as the boundary inputs for MEM to calculate and determine the in situ stress distribution in the model. Stress data of the other two wells, that is, XY-H3 and XY-H6, are used to verify the accuracy of simulation results. Rock properties including elasticity modulus and Poisson's ratio are determined by calibrating with the logging data for different layers, as shown in Figure 2. Figure 3 shows the anisotropy of rock elasticity and its compressive strength with respect to the bedding plane. It can be seen that the rock elasticity in the direction of the bedding plane is approximately 1.3 times larger than that normal to the bedding plane. 29,30 That is consistent with previous experimental reports in the literature ( [31][32][33] in which the modulus of elasticity in the direction parallel to the bedding plane is larger (stiffer) than the modulus perpendicular to the bedding plane. However, due to the weakness along the bedding plane, shale is prone to shear failure along the bedding plane under compression loadings. Our experimental results show that compressive strength in the direction of the bedding plane is only 70% of that normal to the bedding plane. Hence, rock anisotropy is incorporated into the model to achieve more realistic outcomes. Figure 4 shows in situ stress field calculated numerically for the reservoir from 2200 m to 2400 m depth based on existing geological features found in seismic data and anisotropic rock properties derived from logging data and laboratory triaxial tests of different samples with various angles with respect to bedding planes. It can be found the minimum and maximum horizontal stresses for XY-H1 well that suffered severe casing deformation are 34.7 MPa and 68.5 MPa. XY-H1 well is located on the anticline structure, and the elastic modulus in the reservoir is relatively low but the Poisson's ratio is high, which makes the maximum principal stress parallel to the direction of the anticline (NE-SW direction).

FRACTURING FOR CASING FAILURE
Based on the in situ stress field obtained for our case study, some well-defined analyses at a reservoir scale are conducted to mimic the situation led to casing deformation during hydraulic | 1527 fracturing. XY-H1 well is a horizontal shale gas well, whose measured depth (MD) at the bottom is 3250 m and its lateral section is 1050 m. The diameter of the horizontal section is 215.9 mm (production casing is Ф139.7 mm × 9.17 mm with P110 steel grade). Ten-stage fracturing operations were performed from 2250 m to 3250 m MD, and uniform spacing of 100 m. During each stage of the treatment, about 2200 m 3 fluid carrying 80 tons of quartz sand was injected into the formation at an average pumping pressure of 70 MPa and average pumping rate of 25 m 3 /min. Before fracturing, completion tools could pass through the casing very smoothly without any sticking issues. While after the tenth stage, a milling shoe with a diameter of 116.5 mm was run into the casing to drill bridge plugs but stuck at the measured depth (MD) of 2450.6 m (corresponding to 9th stage interval between eighth and ninth stage perforations) that indicates occurrence of severe casing deformations. It can be guessed that a chain of events happened during the fracturing operation are the key factors to cause casing failure. Hence, the following modeling and discussion are closely combined with fracturing practices in the field. position of casing failure as these stimulation zones may intensify the load nonuniformity on casing. Different dot colors in the figure represent the seismic signals recorded in the adjacent well (X-2 well) during each stage to give a rough estimate of the hydraulic fracture geometry or to be more specific stimulation volume vs time. It can be seen that stimulated zones are developing on either side of the well. It is notable that the shape of stimulation regions shows strong asymmetry with respect to the borehole. However, microseismic clouds are recorded on both sides of the borehole, which confirms that asymmetry is not imposed by the receivers' position with respect to the treated well. Asymmetry of stimulated zones could induce strong shear stresses on the casing string. This asymmetry could be due to the lack of stability in hydraulic fracture growth in presence of natural fractures. 34 Additionally, some regions due to the presence of extended natural fractures 5,35 may have received fluids during two or more stages of the treatment, which may weaken its overall stiffness and mechanical strength especially the tension and shear strengths. According to the field exercise described in the previous section, the position of the casing failure coincides with the location of the most asymmetric fracture zones.

| Estimating stimulated rock volume
Since the research objective of this paper is focused on casing deformation, it is not necessary to accurately simulate the detailed process of fracture network propagation and evolution that can be a long and computationally expensive simulation. 36,37 Hence, we propose to estimate SRV of each stage using the corresponding cloud of clustered microseismic events (not sparse events) as inclusions and determine its implications on the stress field near the borehole. Some simple geometries like rectangular cubes and ellipsoids are used to represent the effective coverage of microseismic events by removing dispersed events to obtain an almost continuous volume for modeling and meshing SRV and surrounding rocks. To accelerate recognition and highlighting SRVs, clustering algorithms can be applied on microseismic events, and then, minimum square error are used to find the best size that fits the data. These regions can be recognized as the zones enriched with complicated networks of fractures; therefore, the overall rock stiffness and strength should be deteriorated while pore pressure is augmented. The overlapping regions are assumed to be distinguishable with repeating events and perhaps denser population of microseismic events due to formation of more fractures.
In order to describe the mechanical property of stimulated region, we define the dimensionless fracture spatial density 3d inside stimulation zones using the effective medium theory, [38][39][40] and it can be expressed as where V is the volume of stimulated zone; n is the number of fractures in the region; l i is the average radius of the i th facture. Then, the effective elastic modulus of each stimulation zone containing n randomly oriented, three-dimensional fractures 38 can be calculated as where E 0 and 0 are the undamaged elasticity modulus and Poisson's ratio of the intact rock. It can be known that the overall stiffness of stimulation region is gradually degrading with the increasing number of fractures per certain volume. For those overlapping regions, the stiffness may further decrease due to the presence of more fractures.

| Finite element analysis
Following the described methodology, a reservoir-scale 3D finite element model coupled hydro-mechanic effect is built to simulate multistage fracturing and shown in Figure 6. As shown in Figure 6, the dimension of FE model is 700 m long (ranging from 2200 m to 2900 m measured depth), 500 m wide and 200 m thick. The middle layer is the target formation that contains the stimulation zones from 8th to 10th stage of fracturing. This reservoir-scale model is tied to a more refined mesh for the near-wellbore part including the casing and cement sheath to simulate casing deformations. In other words, SRVs are considered as inclusions with different pore pressure and degraded elastic moduli, an approach similar to what is assumed in micro-poroelasticity, 41,42 however, analytical solution of micro-plasticity solutions is mostly limited to one or two interacting inclusions, so it can be only used as a benchmark of the numerical solution for the simple cases. Figure 7 demonstrates our approach in converting microseismic clouds to stimulation zone inclusions in the FE model. All the events from 8th to 10th stages of fracturing are presented. Then, we partition the target formation using the basic clustering techniques of data points illustrated in Figure 5 for the stimulation regions of the corresponding fracturing operations. The repeating-fracturing regions are also marked to show zones stimulated several times. Before fracturing, these stimulation regions were sharing the same property as the rest of the formation. Once effect of each stage is chronologically considered, the corresponding zones will be updated with new pore pressure and averaged damaged mechanical properties.

| Governing equations
The reservoir rock is considered as a porous media consisting of solid matrix and fully saturated micro-pores. Hence, the fluid flow inside the formation obeys Darcy's law as following F I G U R E 7 Modeling stimulation zones based on microseismic inversion. Pink regions represent events from 8th stage. Gray regions represent events from 9th stage. Yellow regions represent events from 10th stage. Other color regions like dark blue, green, and red represent repeating-fracturing events

Stimulated reservoir v olume
where v w is fluid velocity vector relative to the solid phase; n w is the void ratio; g is gravitational acceleration vector; w is the mass density of fracturing fluid; k is the permeability tensor. The fracturing fluid is assumed to be an incompressible fluid. The mass balance equation for the porous media based on the principle of virtual work can be written as The continuity equation of fluid flow inside the formation can be expressed as where X is spatial coordinate vector; is effective stress tensor; p w is pore pressure; is virtual strain rate tensor; f is surface force vector; v is virtual velocity vector; is the gravity vector; J is formation volume change rate, which is the ratio of the formation's volume in the current configuration to its volume in the original configuration.
During each stage of fracturing, fracturing fluid is injected into the reservoir through the perforation holes along the casing. According to the fluid mass balance, a part of the injected fluid at a certain period fills the stimulation zones enriched with a network of fractures and the rest will be lost to the surrounding rock matrix as described in the following equation where V f is fracture volume; q is the mass flow rate along each fracture; q t and q b are normal flow loss rates into the top and bottom layers of the reservoir, respectively; q inj is the injection rate at a known location denoted by the Dirac delta function, (x,y). Then, the governing equation for rock matrix involves coupling fluid flow and rock deformation as following Six faces of the model are constrained with respect to normal displacements, and in situ stresses are assigned as initial stresses into the whole model before fracturing using the results from the MEM model. A three-step analysis is established to simulate the 8th to 10th stages of fracturing. During each stage, the corresponding stimulation regions will be activated and rock stiffness degrades gradually according to Equations (1) and (2). In the meantime, pressure inside the casing build up in the parts after the previous bridge plug, and the fracturing fluid begins to flow into the stimulation zones and leak off into the rest of the formation based on Equations (3)- (7). After a period of pumping, pore pressure in the reservoir is expected to build up and the stress field of the whole model is brought to another equilibrium status. In summary, when 8th stage is analyzed, stimulation zones corresponding to 9th and 10th stages are treated as intact rock or the deactivated status. The stimulation regions of 9th and 10th stage would be activated at the end of the prior stage and consequently the above-mentioned procedures would be repeated, that is, simulating the pumping and adjusting rock properties for representing damage in the stimulation zones. Figure 8 demonstrates number of times each part of the formation has received the treatment in different stages. It can be seen that the stimulation regions are continuously expanding stage by stage, while there are some overlappings. The degree of stimulation for each part of the reservoir after the 10th stage of fracturing can be classified into four different categories: no stimulation or fracturing, stimulated or fractured once, twice or three times stimulated.

RESULTS AND DISCUSSION
During 8th stage of the treatment, high-pressure fracturing fluid flows into the identified SRV along the wellbore and increases pore pressure in these regions. As a result, the  increasing fluid pressure would couple with the in situ stress induced by tectonics and other stages to significantly redistribute and re-orient it with respect to the initial principal stresses. The rock cracking and fragmentation inside stimulated regions would also intensify stress heterogeneity inside the stimulated zones. From Figure 9, it can be found that the maximum values of σ H , σ v, and σ h have been increased to 45.8 MPa, 53.8 MPa, and 99.2 MPa around stimulation regions (as marked in the arrowed position of Figure 9). However, stress field inside stimulation zones has decreased to small values or even positive values (tensile stress) due to the offsetting effect of the pore pressure, as highlighted in the areas enclosed by red curves. The zero or tensile stress zones are beneficial to maintain fractures and their intersections open for proppant placement. But at the same time, these zones will cause the production casing temporarily to lose its support and suspend in the formation in some extent. Hence, this section of casing is prone to transverse or axial deformation due to the shear stress induced by the asymmetric development of stimulation zones. Although the stress field may gradually restore its initial status with the flowback of fracturing fluid, the damage on the casing caused by the fracturing treatment can be irreversible. After experiencing 9th and 10th stage of fracturing, extensive tensile stress areas continue to spread and coalesce with each other, as shown in Figures 10 and 11. The stimulation pattern is becoming more and more asymmetric as well, in which some regions have received twice or more treatments. All these have resulted in a drastic redistribution and re-orientation of the stress field. The maximum values of σ H , σ v, and σ h have been increased to 57.7 MPa, 69.7 MPa, and 132.6 MPa around the stimulated regions. Large tensile stress areas indicate that the purpose of Stimulated Reservoir Volume (SRV) has been realized. However, it should be noted that the induced shear force has also increased, which makes the production casing suffer more significant shearing and therefore more prone to failure. According to the field reports, milling shoe was stuck at the severely deformed casing right after the tenth stage, which confirms the impact of induced shear force on casing failure. Figure 12 shows the displacement vector contours after several stages. It can be seen that the formation movement is becoming stronger after each stage and the maximum displacement has increased to 276.8 mm after 10th stage of fracturing. Additionally, we can find that not only the regional displacement magnitudes are different from each other, but also the direction of displacements in different regions are not the same and in some cases completely opposite, which may induce strong shearing forces on the casing.

F I G U R E 8
To illustrate how the external load on casing may change after multistage fracturing, we obtained the contact pressure contour on casing and corresponding curve of casing along well depth after 10th stage fracturing, as shown in Figure  13. It should be noted that face 2 (green section in Figure  13A) is the plane where casing failure occurs. The distribution of contact pressure on casing is even more uniform than that in the adjacent faces, which means nonuniformity of external load on casing is not the immediate cause for large deformation. Figure 13B demonstrates a clear curve of contact pressure along the casing. It can be seen that the overall magnitude of contact pressure within the scope of SRV has decreased compare to the initial status (section from 2200 m to 2250 m is not subject to hydraulic fracturing and can represent the initial contact pressure before fracturing). That is because the elevated pore pressure in SRV offsets the in situ stresses and reduce the effective stresses, which is consistent with the previous discussions (as demonstrated in Figures  8-11). However, the nonuniformity of contact pressure has been intensified obviously, especially on both sides of casing failure location (as shown in tinge green and orange areas). Combined with Figure 12, it can be known that directions of force in these two regions are completely opposite which could have verified the strong shearing load on casing string.  Figure 14 shows the Mises stress distribution of casing after different stages of hydraulic fracturing. It can be seen that maximum stress is gradually increasing with multistage fracturing due to further formation shearing induced by more asymmetric stimulated zones (as illustrated in Figures 11 and  12), and the maximum of von Mises stress reaches 773.8 MPa by finishing the 10th stage of fracturing that indicates parts of P110 grade casing has entered into the yield conditions and suffered some irreversible plastic deformations. In order to quantitatively evaluate the casing deformation in the axial direction, the corresponding displacements along the measured depth are extracted from numerical simulations and shown in  Figure 15B. To gain a better understanding on the bending severity of casing in a short distance, the curvatures along the measured depth are also calculated to reflect the amount of deviation from the straight-line direction, which is like the concept of dogleg severity for well trajectory (as shown in Figure 15A,C). The calculation for bending curvature of sectional casing string can be expressed as where K is the curvature at some point of deformed casing; ∆S is the arc length of sectional casing string; Δ is the tangent angle corresponding to the arc; ̇y is the first derivative of curvilinear equation; ÿ is the second derivative of curvilinear equation. Here, we used interpolation method to process the data of discrete nodes and calculate the corresponding derivations. From Figure 15, it can be seen that the displacement and curvature of casing are very small before the fracturing, which indicates that the casing is mechanically undamaged. As multistage fracturing goes on, the severity of casing deformation and curvature keeps increasing due to the shearing effect. After the 10th stage of fracturing, the maximums of displacement and curvature have increased to 129.3 mm and 1341.3/10000, which indicates large casing deformations along a short section. Furthermore, there is even a reversed bending phenomenon occurring near 2450 m, which makes the casing deformed in the S-shape. In this situation, rigid completion tools such as milling shoe can hardly pass through the deformed casing. Figure 16 illustrates transverse deformations of the casing along the measured depth after 10th stage of fracturing. It can be found that plastic deformations occur on parts of casing due to the large axial deformation, and the maximum of equivalent plastic stain is 0.214% near 2750 m. The maximum and minimum inner-diameters of deformed casing are 130.2 mm and 111.3 mm occurring at the same location in the proximity of 2451 m, which makes casing deform like an ellipse in the transverse direction. Hence, after the tenth fracturing stage, milling shoe with a diameter of 116.5 mm could not pass such a narrow pipe smoothly and stuck there. In the field practice, multi-finger imaging tool is used to measure the casing deformation and diameter change of casing, as shown in Figure 17A,B. It can be seen that obvious casing deformations occur from 2448 m to 2452 m. Lead impression block is subsequently run into the casing to print the shape of cross section, as shown in Figure 17C. The deforma-good consistency with numerical simulation results, which verify the reliability of the modeling in this paper.
Based on the above studies about the casing deformation during a multistage fracturing treatment, we found that the failure mechanism is the large deformation due to the shearing effect induced by asymmetry of the stimulation zones. The repetitive intense fracturing would intensify the resultant shear forces to exacerbate casing deformation; however, we cannot deny the potential impacts of geological uncertainties such as pre-existing discontinuities. Increasing the flexural strength of casing by increasing steel grade and wall thickness could enhance the resistance to bending to a certain extent, but it cannot solve casing deformation problem substantially. Therefore, the best solution is to optimize the spacing design of multistage fracturing in situations that there is a risk of developing asymmetric SRVs. Therefore, the spacing design should be adjusted dynamically based on the real-time microseismic monitoring to make the stimulation zones as symmetric as possible on the premise of achieving industrial production of SRV.

| CONCLUSIONS
The significant impact of asymmetric stimulation zones induced by multistage fracturing on casing failure in shales is already highlighted in the petroleum and natural gas engineering literature; however, current models only consider rock deformations very close to the wellbore based on various hypotheses due to the complex fracture network. Though, none of these models have discussed the in situ stress field before fracturing at a geological scale and evaluated stress redistribution after fracturing and its impact on casing deformations at the reservoir scale. To address these problems, this paper provides a comprehensive workflow to investigate the potential impact of asymmetric stimulated rock volumes on casing deformations during multistage fracturing treatment. Main outcomes can be drawn as follows: We initially developed a three-dimensional mechanical earth model incorporating geological features and rock physics features including rock anisotropic properties to