Optimizing Hydraulic Fracturing Treatment Integrating Geomechanical Analysis and Reservoir Simulation for a Fractured Tight Gas Reservoir, Tarim Basin, China

A comprehensive geomechanical study was carried out to optimize stimulation for a frac‐ tured tight gas reservoir in the northwest Tarim Basin. Conventional gel fracturing and acid‐ izing operations carried out in the field previously failed to yield the expected productivity. The objective of this study was to assess the effectiveness of slickwater or low-viscosity stim‐ ulation of natural fractures by shear slippage, creating a conductive, complex fracture net‐ work. This type of stimulation is proven to successfully exploit shale gas resources in many fields in the United States. A field-scale geomechanical model was built using core, well log, drilling data and experien‐ ces characterizing the in-situ stress, pore pressure and rock mechanical properties in both overburden and reservoir sections. Borehole image data collected in three offset wells were used to characterize the in-situ natural fracture system in the reservoir. The pressure re‐ quired to stimulate the natural fracture systems by shear slippage in the current stress field was predicted. The injection of low-viscosity slickwater was simulated and the resulting shape of the stimulated reservoir volume was predicted using a dual-porosity, dual-permea‐ bility finite-difference flow simulator with anisotropic, pressure-sensitive reservoir proper‐ ties. A hydraulic fracturing design and evaluation simulator was used to model the geometry and conductivity of the principal hydraulic fracture filled with proppant. Fracture growth in the presence of the lithology-based stress contrast and rock properties was com‐ puted, taking into account leakage of the injected fluid into the stimulated reservoir volume © 2013 Gui et al.; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. predicted previously by reservoir simulation. It was found that four-stage fracturing was necessary to cover the entire reservoir thickness. Post-stimulation gas production was then predicted using the geometry and conductivity of the four propped fractures and the en‐ hanced permeability in the simulated volume due to shear slippage of natural fractures, us‐ ing a dual-porosity, dual-permeability reservoir simulator. For the purpose of comparison, a conventional gel fracturing treatment was also designed for the same well. It was found that two-stage gel fracturing was sufficient to cover the whole reservoir thickness. The gas production profile including these two propped fractures was also estimated using the reservoir simulator. The modeling comparison shows that the average gas flow rate after slickwater or low-vis‐ cosity treatment could be as much as three times greater than the rate after gel fracturing. It was therefore decided to conduct the slickwater treatment in the well. Due to some opera‐ tional complexities, the full stage 1 slickwater treatment could not be executed in the bottom zone and treatments in the other three zones have not been completed. However, the posttreatment production test results are very promising. The lessons learned in the planning, design, execution and production stages are expected to be a valuable guide for future treat‐ ments in the same field and elsewhere.

A field-scale geomechanical model was built using core, well log, drilling data and experiences characterizing the in-situ stress, pore pressure and rock mechanical properties in both overburden and reservoir sections. Borehole image data collected in three offset wells were used to characterize the in-situ natural fracture system in the reservoir. The pressure required to stimulate the natural fracture systems by shear slippage in the current stress field was predicted. The injection of low-viscosity slickwater was simulated and the resulting shape of the stimulated reservoir volume was predicted using a dual-porosity, dual-permeability finite-difference flow simulator with anisotropic, pressure-sensitive reservoir properties. A hydraulic fracturing design and evaluation simulator was used to model the geometry and conductivity of the principal hydraulic fracture filled with proppant. Fracture growth in the presence of the lithology-based stress contrast and rock properties was computed, taking into account leakage of the injected fluid into the stimulated reservoir volume

Introduction
Following the success in exploiting shale gas resources by multi-stage hydraulic fracturing with slickwaters or low-viscosity fluid (i.e., linear gel) in horizontal wells in North America, there has been a lot of interest in applying this technique to other regions and other types of tight reservoirs. This is due in part to the fact that conventional gel fracturing treatments have been less successful in some naturally fractured reservoirs due to excessive unexpected fluid loss and proppant bridging in natural fractures, leading often to premature screen-outs. Additionally, the high-viscosity gel left inside the natural fractures causes the loss of virgin permeability of the reservoir in the case of inefficient gel breaking.
However, the challenge for doing this is that the physical mechanism responsible for this kind of stimulation is yet to be fully understood and a standard work flow for design and evaluation is yet to be developed. Furthermore, industry so far mainly relies on performance analogs to improve understanding of each shale play, and thus it usually takes years to advance up the learning curve for determining which factors best affect well production [1].
Currently, the general opinion on the mechanism leading to the success of waterfrac in shale gas reservoirs is that a complex fracture network is created by stimulation of pre-existing natural fractures. Although it is difficult to observe the processes acting during stimulation, microseismic imaging has enabled us to understand that both simple, planar fractures and complex fracture networks can be created in hydraulic fracture stimulations under different settings [2]. Fracture complexity is thought to be enhanced when pre-existing fractures are oriented at an angle to the maximum stress direction, or when both horizontal stresses and horizontal stress anisotropy are low, because these combinations of stress and natural fractures allow fractures in multiple orientations to be stimulated [3]. The result of stimulation therefore depends both on the geometry of the pre-existing fracture systems and on the in-situ stress state. It is now generally accepted that stimulation in shale gas reservoirs occurs through a combination of shear slip and opening of pre-existing (closed) fractures and the creation of new hydraulic (tensile) fractures [4][5][6]. In wells that are drilled along the minimum horizontal stress (S hmin ) direction, stimulation generally creates a primary radial hydraulic fracture that is perpendicular to S hmin . Then, pressure changes caused by fluid diffusion into the surrounding rock and the modified near-fracture stress field induced by fracture opening cause shear slip on pre-existing natural fractures. If the horizontal stress difference is small enough, new hydraulic fractures perpendicular to the main fracture can open. Each slip or oblique opening event radiates seismic energy, which, if the event is large enough, can be detected using downhole or surface geophones.
Founded on the idea that productivity enhancement due to stimulation results not just from creation of new hydraulic fractures but also from the effect of the stimulation on pre-existing fractures (joints and small faults), a new workflow dubbed "shale engineering", was established by combining surface and downhole seismic, petro-physical, microseismic, stimulation, and production data [7,8]. In this new workflow (Figure 1), the change in flow properties of natural fractures is predicted using a comprehensive geomechanical model based on the concept of critically stressed fractures [9][10][11]. Existing reservoir simulation tools can then be used to model the hysteresis of fracture flow properties that result from the microseismically detectable shear slip, which is critical to the permanent enhancement in flow properties and increased access to the reservoir that results from stimulation. The primary hydraulic fracture created and propped during the stimulation can be modeled using conventional commercial hydraulic fracture models by taking into account fluid leaked into natural fractures in the surrounding region. The propped conductivity is estimated using laboratory-based proppant conductivity data adjusted for the proppant concentration in the fracture. The propped main fracture model and the reservoir model with stimulated natural fracture properties can then be integrated into production simulators to predict production after the slickwater hydraulic fracturing treatment. When available, microseismic data can be used to help define the network of stimulated natural fractures that comprises the stimulated reservoir volume (SRV).
Although this new workflow was developed based on experiences in shale gas reservoirs, we believe it can also be applied to any unconventional reservoir requiring stimulation that has pre-existing natural fractures. Both Coal Bed Methane (CBM) and fractured tight gas reservoirs are examples of where this approach could be applied. In this paper, we will illustrate the workflow using the results of a study conducted in a fractured tight gas reservoir in the Kuqa Depression, Tarim Basin.

Project background
The project discussed in this paper was initiated to investigate various methods and practices to improve the economics of the field. Conventional gel fracturing had been tested in a few wells with disappointing results. One fault block (see Figure 2) was chosen as the target of a pilot study that included building a geomechanical model, optimizing hydraulic fracturing design, assessing the stability of faults near the target well, and (although it is not discussed here) analyzing wellbore stability for drilling horizontal wells. Three vertical wells were drilled. D2 and D3 are near the crest of the structure and D1 is ~ 2.5-3 km to the west. The main target is Cretaceous tight sandstone occurring at ~5300m to ~6000m depth. Reservoir rock is composed of fine sandstone and siltstones interlayered with thin shales. Average reservoir porosity is ~7% and average permeability is ~0.07 mD. The gross reservoir thickness is 180-220m in this fault block. Wells D1 and D2 were completed by acidizing and gel fracturing; test production was ~15-27 ×10 4 m 3 /d. The objective of this project was to optimize hydraulic fracturing design for Well D3 based on the geomechanical analysis and investigate whether it is better to conduct slickwater treatment in the D3 well to stimulate and create a complex fracture network or utilize conventional two-wing gel fracturing.
Comprehensive datasets were available for all three wells including drilling experiences, wireline logs, image data, mini-fracs and well tests. Laboratory tests were also conducted on cores from well D2 to estimate the rock mechanical properties of reservoir rocks.

Geomechanical model
A geomechanical model includes a description of in-situ stresses and of rock mechanical and structural properties. The key components include three principal stresses (vertical stress (S v ), maximum horizontal stress (S Hmax ) and minimum horizontal stress (S hmin )), pore pressure (Pp) and rock mechanical properties, such as elastic properties, uniaxial compressive strength (UCS) and internal friction. The relative magnitude of the three principal stresses and the consequent orientation of the most likely slipping fault or fracture define the stress regime to be normal faulting (S v >S Hmax >S hmin ), strike-slip faulting (S Hmax >S v >S hmin ) or reverse faulting (S Hmax >S hmin >S v ). The horizontal stresses are highest relative to the vertical stress in a reverse faulting regime and lowest relative to the vertical stress in a normal faulting regime. Hydraulic fractures are vertical and propagate in the direction of the greatest horizontal stress in a strike- slip or normal faulting regime. In a reverse faulting stress regime in which S v is the minimum stress, hydrofractures are horizontal. These different stress regimes also have consequences for the pressure that is required to open a network of orthogonal hydrofractures by stimulation. In places where the horizontal stresses are low and nearly equal, a relatively small excess pressure above the least stress may be required to open orthogonal fractures. Where the horizontal stress difference is larger, a larger excess pressure is required to open orthogonal fractures. Where the least stress is only slightly less than the vertical stress, weak horizontal bed boundaries and mechanical properties contrasts between layers may allow opening during stimulation of horizontal bedding ("T-fractures"). Except for the magnitude of S Hmax , other components of the geomechanical model can be determined using borehole data by reviewing a few representative wells in the field. Vertical stress is calculated by integrating formation density, which is obtained from wireline logs. The magnitude of S v across this fault block is in a similar range. Pore pressure was constrained, mainly by referencing direct measurement data and drilling experiences. This is due to the complex tectonic history. Conventional under-compaction approaches for pore pressure estimation may not apply in the study area. Evidence for this is the over-compacted density profile. In addition, due to the complex lithology changes the log response with depth may reflect lithology changes rather than pressure variation. Well test data from D1 and D2 showed that the reservoir pressure is ~88-90 MPa, an equivalent pressure gradient of ~1.6-1.7 SG, which is abnormally over-pressured.
Rock mechanical laboratory tests were conducted on cores from the sandstone reservoirs and the interlayered shales in the D2 well, and the results were used to constrain a log-calibrated range of UCS and other rock mechanical parameters. Figure 3 shows the match between logderived rock strength profiles and laboratory test results in D2. Dynamic Young's modulus was calculated from compressional and shear velocities and density and calibrated to static values using laboratory test results. The relationship between dynamic and static Poisson's Ratio was not obvious; the dynamic Poisson's Ratio computed from Vp/Vs matched reasonably well with the laboratory results, so it was used directly in the modeling. Young's Modulusbased empirical relationships were used to estimate the UCS for both sandstone and interlayered shales.
Minimum horizontal stress (S hmin ) at depth can be directly estimated from extended leakoff tests (XLOT), leak-off tests (LOT) or mini-frac tests. No extended leak-off tests were conducted in the field. LOTs and leak-off points from two reliable LOTs were used to constrain the upper limit of S hmin (~2.09 SG EMW at ~4000 m TVD). One mini-frac test was conducted in the sandstone reservoir in D2, with the interpreted fracture closure pressure (closest estimation to S hmin ) ~2.064 ppg EMW at ~5400 m TVD. Because LOTs are usually conducted in shaly formations while mini-frac tests are usually carried out in sandstone reservoirs, the LOTs and mini-frac tests are used to construct separate S hmin profiles in shales and sandstones, respectively using the effective stress ratio method (S hmin -Pp/S v -Pp). The effective stress ratio from LOT is ~0.725 and from mini-frac test is ~0.48, which indicates there is a dramatic stress difference between sandstones and shales (stress contrast). The contrast between different lithology significantly influences hydraulic fracturing design. The relative lower stress in sandstones indicates that a hydraulic fracture should be easily created in the tight sandstone, however, the interlayered shales which have higher stress act as frac barriers and pinch points, thereby complicating fracture propagation and the final fracture geometry and conductivity.
The azimuth and magnitude of maximum horizontal stress (S Hmax ) can be constrained through the analysis of wellbore failures such as breakouts and tensile cracks observed on wellbore images or multi-arm caliper data. Wellbore failure analysis allows constraining of the orientation and magnitude of the S Hmax because stress-induced wellbore failures occur due to the stress concentration acting around the wellbore once is drilled. The presence, orientation, and severity of failure are a function of the in-situ stress fields, wellbore orientation, wellbore and formation pressures and rock strength [12]. High-resolution electrical wireline image logs were available in all three study wells. Both breakouts and drilling-induced tensile fractures (DITFs) were observed in the reservoir sections in D2 and D3 wells. Only DITFs were observed in well D1, which could be due to the higher mud weights used during drilling and the poor quality of the image data in lower part of the reservoir. Figure 4 shows examples of the breakouts seen in the D3 well. The example shows the typical appearance of breakouts observed on images. Here, the average apparent breakout width is ~30-40 degree. The breakouts mostly occur in shales and more breakouts are observed in the lower part of reservoir where the formations become more shaly. The orientation of breakouts is quite consistent with depth and across the block. However, small fluctuations of breakout orientation can be observed locally while intercepting small faults (an example can be seen in the right plot in Figure 4). This may indicate that some of these faults are close to or at the stage of being critically stressed. This has important implications for the stress state in the area and the likelihood of stimulating fractures by injection. Breakouts usually develop at the orientation of S hmin and DITFs in the direction of S Hmax in vertical and near-vertical wells. In the left plot of Figure 4, DITFs can also be observed in the same interval as the breakouts with an orientation that is ~90 degrees from the breakout directions, consistent with this expectation. DITFs are seen more often in sandstone than in shale. Based on wellbore breakouts and DITFs interpreted from the image data in D3, the azimuth of S Hmax is inferred to be ~143° ±10 °. This is similar to the azimuth of S Hmax inferred from wellbore failures observed in the other two wells. It is also consistent with the regional stress orientation from the World Stress Map [13].
The magnitude of S Hmax is constrained by forward-modeling the stress conditions that are consistent with observations of wellbore failures observed on image logs, given the data on rock strength, pore pressure, minimum horizontal stress, vertical stress, and mud weight used to drill the well.  magnitudes in the different lithologies. However, both results are consistent with the magnitudes of S hmin inferred from LOTs and mini-fracs. Figure 5 shows that the magnitude of maximum horizontal stress is higher than the vertical stress in both cases, and higher in the shale than in the sand. Thus, the study area is in a strike-slip faulting stress regime (S hmin < S v < S Hmax ). The difference between the magnitudes of S Hmax and S hmin is ~0.8 SG in the reservoir section, suggesting high horizontal stress anisotropy. In such a condition, it is unlikely to open the natural fractures by tensile mode. However, the natural fractures might dilate in shear mode depending on their orientations and stress conditions. The final geomechanical model was verified by matching the predicted wellbore failure in these wells with that observed from image data and drilling experiences.

Natural fractures characterization and stimulation modeling
Natural fractures have been observed on cores and image logs in the study area. The fluid losses during drilling not only suggest the existence of natural fractures but also that some at least of these fractures are permeable in-situ. Based on the core photos shown in Figure 6, open high-angle tectonic fractures can be seen on cores from D2 and D3 wells near the crest of the structure. A fracture network consisting of a group of fractures with different orientations can be seen on the cores from the D1 well, and these fractures appear to have less apertures than high-angle fractures observed in D2 and D3.  Natural fractures were interpreted and classified using high-resolution electrical images in all three wells. Based on the appearance on image data, the natural fractures are classified as below: • Conductive: dark highly dipping planes on image logs  The plot on the right shows an example of drilling enhanced natural fractures for which the fracture trace is discontinuous. The fact that parts of these fractures can be detected on the electrical image is due to fluid penetration into the fracture at the orientation around where the rock at the borehole wall is in tension during drilling. The classification of the natural fractures indicates the relative strength of the fractures. For example, the resistive fractures are closed and mineralized. Active faults or critically stressed natural fractures might be open and conductive, even under the original conditions. During stimulation, these fractures are the most easily stimulated. However, it is important to note that the classification of natural fractures is purely based on their appearance on the electrical images, and cannot be used directly to quantify permeability or other flow properties. Figure 8 shows the fractures orientations on a crossplot of the strike and dip angles of all fractures observed in the three wells. The natural fractures observed can be divided into three groups. The first group is low-angle fractures (dip<20 °), which could be related to beddings. The second group is the major fractures seen in this block that have intermediate dip angles (~25-55°) and strike at an azimuth of ~155°N. The third group consists of fractures with strikes of ~355°N and ~100°N and dip angles ~ 35°-65° and 25°-35°, respectively. Because of their wide range of orientations and cross-cutting relationship, these three groups of fractures could be stimulated to form a well-connected grid with a major fracture azimuth (~155°N) aligned with the direction of maximum horizontal stress (~143°N). This direction is nearly perpendicular to the faults, defining the shape of fault block (see Figure 2). Because the structural trends and the stresses are aligned, it enabled us to create a reservoir model with a grid that is consistent with both. Effective stresses in the earth are always compressive, and natural processes tend to "heal" fractures through vein filling and other processes. Therefore, the intrinsic fracture aperture of most fractures is likely to be very small or even zero (cases where dissolution creates voids that prevent full closure are a notable exception). Thus, it is increasingly recognized that active processes are necessary to maintain fracture permeability. One such process is periodic slip along fractures that are critically stressed (i.e., those that are at or near the limiting ratio of shear to normal stress to slip). This process, and the influence of effective normal stresses on fracture aperture, can be modeled using a simple equation that describes the variation in aperture as a function of normal stress for a pure Mode I fracture. The same equation with different parameters can also be used to model the same fracture after slip has occurred [9][10][11].
Equation 1 is one example that describes aperture in terms of an initial aperture (A • a 0 ) and an effective normal stress at which the aperture is only 10% as large (B). A and B both increase due to slip, resulting in a larger "unstressed" aperture and a stiffer fracture caused by "selfpropping" due to generation during slip of a mismatch in the fracture faces and/or creation of minor amounts of rubble at the fracture face.
The contribution of fractures to the relative productivity of a well of any orientation can be computed by summing the contributions of all fractures, weighted by the product of their relative transmissivity (which is a function of aperture) and the likelihood of the well inter- secting the fracture (which is a function of the difference between the fracture and the well orientation). This relative productivity can be written as [10] P well = ∑ fracs {max (|ŵ • n î |, a) × P i } (2) where ŵ and n î are unit vectors along the axis of the well and normal to the i th fracture, a is a number representing the likelihood of a well intersecting a fracture if it lies in the plane of the fracture, and P i is the relative permeability of the fracture.
The fractures interpreted from image data are only those that intersect the logged wells that are a function of their orientations, and there is no information about the fracture distribution between the wells. To ensure the most meaningful representation of the fractures in the reservoir, the fractures interpreted from all three wells were combined and the distribution was corrected to account for the likelihood of each fracture intersecting the well at the point where it was observed. This combined fracture data set was then used to model the productivities of wells in their natural condition and the change in productivity due to the shear-slip of natural fractures. Effective and Sustainable Hydraulic Fracturing Figure 9 shows relative productivity for wells of all orientations based on the fractures observed in all three wells. Natural fractures are shown as poles to the fracture surfaces (black dots). Different apertures and strengths were assumed for the different types of fractures based on their classifications described above (Table 1). The plot on the left shows the relative productivity under pre-stimulation conditions, while the plot on the right shows the relative productivity calculated using equation 2 after the fractures were stimulated with a pressure 20 MPa above the original reservoir pressure. It can be seen that the maximum productivity increases by a factor of 5 if all fractures see the same 20-MPa pressure increase, which is obviously not the case during real stimulation. Superimposed on Figure 9 are the computed optimal orientations of wells based on the fracture and stress analysis (green circles). If none of the fractures is critically stressed, then the best orientation to drill a well is perpendicular to the largest population of natural fractures. If some fractures have enhanced permeability because they are critically stressed, the optimal orientation shifts in the direction of the greatest concentration of critically stressed fractures. Figure 9 shows there are some fractures already near or being critically stressed, even under ambient condition (left plot), and the maximum productivity is achieved by drilling highly deviated wells with ~20 °N hole azimuth. The optimum wellbore orientation after ~ 20-MPa stimulation is nearly horizontal and in the direction of ~228 °N.    The above relative productivity modeling of natural fractures shows the conductivity of natural fractures increases significantly if the stimulation pressure is at or above the minimum horizontal stress. This is because many of the natural fractures are non-optimally oriented. Assuming a connected fracture network exists, the conductivity increase could be a factor of five for the stimulated fracture network while stimulation pressure is ~130MPa or higher (assuming the pressure reaches all fractures).

Predicting the shape of the stimulated reservoir volume
Fracture stimulation modeling showed that the shear slip of natural fractures could be effective in improving reservoir properties. Next, we need to reproduce the affected productive volume in the reservoir using the "shear stimulation" concept to enable more accurate production prediction. At the present no commercial simulator can fully model this process in 3D, although some research simulators have been developed. It was decided to use two different commercial models to simulate both fracture network stimulation created by low-viscosity frac fluid and the growth of the main hydraulic fracture. A commercial dual-porosity, dualpermeability simulator is used to simulate the flow property changes of natural fractures due to the shear slip. A commercial hydraulic fracturing design and evaluation simulator is used to model the geometry and conductivity of the principal hydraulic fracture filled with proppant. The modeling in two separate simulators is coupled by the fluid volume used for stimulation. The fluid volume leaked off in the shear-dilated natural fracture network was estimated in the dual-permeability, dual-porosity flow simulator. By adjusting the pressuredependent leak-off coefficient, the fluid volume leaked off in the hydraulic fracturing simulator was matched with the fluid volume leaked into natural fractures networks estimated by the flow simulator. The prediction of the stimulated reservoir volume is discussed in the rest of this section and the hydraulic fracturing design will be discussed in next section.
To predict the extent and properties of the stimulated volume by a dual-permeability, dualporosity simulator, a finely gridded model (Model A) was created based on the original reservoir model. The main function of this model is to simulate the change in flow properties in every single frac stage during and immediately after injection. The model is initialized with average known reservoir characteristics such as matrix porosity and permeability, fracture permeability and initial pressure, characterized from core and log analysis. Although different cases have been tested in the study, only one of the most realistic cases will be discussed here: the average matrix porosity used in the initial model is ~7.4%, matrix permeability is 0.07 mD in all directions, and the initial fracture permeability is ~ 0.2 mD. The initial fracture permeability is set close to the lower bound of fracture permeability based on core and log analysis. The orientations of the principal flow directions were chosen to correspond to the principal directions of the fracture sets and of bedding, which also approximately corresponded to the principal stress directions.
The relative magnitudes of the permeability enhancements in different directions were constrained by the geomechanical analysis. A set of permeability-pressure tables for different directions were then used to describe the hysteretic rock behavior that results from shear fracture activation. Although the fracture properties during stimulation can be estimated as described in the previous section, it is better to calibrate and constrain the permeabilitypressure relationship based on real lab or in-situ tests, e.g., using a pre-stimulation injectivity test [4]. The injectivity test should ideally be conducted in the open hole using slow injection to evaluate the potential natural fractures being stimulated, as permeability changes could then be interpreted based on the flow-rate/pressure changes along with the reservoir pressure. Because the D3 well has already been cased it was impossible to conduct such a test in the field before the actual treatment is carried out. Consequently, it was decided to produce a permeability-pressure table based on experience from shale gas reservoirs. Based on this table, on fracture density in different directions and on the stress anisotropy, a composite transmissibility multiplier was produced for the prediction of properties and extent of the stimulated reservoir volume. Transmissibility multipliers were different for each of the I, J and K directions; those directions were aligned as discussed above with the primary structural fabric and stresses. The propagation of the pressure and fluid front in these directions can be controlled by modifying these multipliers. Figure 11 shows diagrammatically the relationship between the permeability multiplier and the pore pressure (green curve). A slow increase in the permeability multiplier with increasing pressure occurs until fractures begin to slip. Above this pressure, the injectivity increases rapidly as an increased number of fractures are stimulated. During decreasing injection pressure in the injectivity test, the injectivity should decrease more slowly, retaining behind a permanent injectivity increase. The post-stimulation response can also be extrapolated to pressures below the original reservoir pressure. This makes it possible to predict the reservoir's response to depletion, which could lead to improved predictions of production decline. When the pressure during stimulation exceeds the minimum horizontal stress, extensional hydrofracs are created, and the permeability-pressure relationship does not follow the green line. Three different flow paths (A, B, C) were assumed for conditions with pressure above S hmin , and the intermediate path, B was chosen to be used in the simulation.
The result of this modeling work is a 3D induced permeability map that describes the stimulated rock volume as discrete blocks, each with a unique permeability. The stimulated rock volume is therefore described not as a geometrical shape with identical flow properties throughout, but as a rock body with variable induced permeability, as shown in Figure 12.

Hydraulic fracturing design and reservoir simulation
As discussed earlier, a commercial simulator was used to model the hydraulic fracture created during the stimulation along with the stimulated natural fracture network using low-viscosity fluids. Stress profiles and other elastic rock properties estimated in the geomechanical analysis were used as input for the design. To achieve better proppant distribution, a low-viscosity linear gel was combined with slickwater in the treatment. The low-viscosity linear gel was optimized using different concentrations of ingredients for the high reservoir temperature (~126°C) using source water and local ingredients. Due to the high closure pressure and low viscosity of the fluid, high-strength small-mesh proppants were used in the design.
Modeling showed that four stages would be required for slickwater/linear gel treatment to cover the 160 m thick reservoir due to the high leak off of low-viscosity fluids ( Figure 13). A reasonable proppant distribution was achieved by using the low-viscosity linear gel.  To predict the production after the stimulation, the propped hydraulic fractures were imported into the reservoir model with flow properties enhanced by stimulated natural fractures (Model A). Because the natural fracture distribution between wells is unknown, the same stimulated Model A was used for all four stages. The left plot of Figure 14 shows a side view of the reservoir model combining four Model A's with stimulated reservoir volumes and four propped hydraulic fractures, which was used for production prediction.
To compare the prediction result from slickwater/liner gel treatment with conventional gel fracturing, a conventional bi-wing hydraulic fracturing design using a high-viscosity gel was also developed. The gel fluid was optimized using different concentrations of ingredients for the high reservoir temperature (~126°C) using source water and local ingredients. The same type of proppant used for the slickwater/liner gel treatment was used for the design of gel treatment. The proppant concentrations and amounts will be certainly different in these two types of treatments. It was found that two stages were enough to cover the whole reservoir interval (Figure 15). These two designed hydraulic fractures were then imported into the original reservoir model (right plot in Figure 14) for production prediction and comparison of the production to that predicted after slickwater linear gel stimulation.  Figure 16 shows the production prediction comparison from the two different hydraulic fracturing treatments. The red curve is the production prediction from slickwater/linear gel treatment, which is scaled down to ~2/3 of the initial prediction to account for the heterogeneity of the reservoir model due to a simplified reservoir model used for pre-stimulation condition. The blue curve is the production from conventional two-wing gel fracturing design. It is found that post-frac flow rate from slickwater stimulation is expected to be about three times the flow rate from the gel treatment in the stabilized regime (one year after stimulation). Although actual flow rates from both treatments depends on the applied drawdown, the corresponding flow rates after one year are expected to be ~55 × 10 4 m 3 /d for slickwater treatment and ~ 17 × 10 4 m 3 /d for gel treatment, respectively, with a constant drawdown of 20 MPa. Figure 16. Production prediction comparison of two different hydraulic fracturing treatments. The red curve is the production prediction from slickwater/linear gel treatment; the blue curve is the production from conventional twowing gel fracturing design.

Injectivity test and stage 1 treatment
It was decided to test the slickwater/liner gel treatment in D3 well after the study was completed. A pre-stimulation injectivity test was performed through perforations prior to Stage 1 and after the mini-frac test (Figure 17). Interestingly, the test showed the opposite behavior from what one would expect if the stimulation enhances reservoir permeability. Later-stage injectivity (during step-down) is lower than early stage injectivity (during step-up), rather than higher. Although there might be other reasons affect the test result, i.e., the un-stable injection during the whole test, it is believed the main reason was lack of access to natural fractures in the tested interval and the high closure pressure because the test was conducted in a cased and perforated hole and after a mini-frac. The Stage 1 treatment was conducted using slickwater and linear gel after the injectivity test. However, a screen out was experienced at the end of the execution and tubing leakage was discovered afterwards. Treatments in the other three zones had not occurred at the date of writing this paper. The stage 1 production test is still very promising, and it has been decided to continue slickwater/linear gel treatment in other three stages after the tubing problem is fixed.

Discussion and conclusion
In this paper we have outlined a new workflow for simulation of a complex fracture network created by stimulation using low-viscosity fluids in a fractured tight sandstone reservoir. The workflow is based on critically stressed fracture theory. This process of natural fracture stimulation is believed to be the underlying reason for the success in shale gas reservoir stimulation. The results suggested that there would be significantly higher production from this approach compared to conventional two-wing gel fracturing.
There are, however, some uncertainties in the modeling of the natural fracture stimulation for this fractured tight gas reservoir.

1.
The pressure-permeability relationship used in modeling the permeability enhancement by slickwater stimulation is taken from a shale gas field. It is unclear whether the data from the analogue field drilled through mudstones will be applicable to the modeled fractured tight sandstone reservoir. Post-stimulation production simulation, or a prestimulation injectivity test in nearby wells in open hole could help to better constrain this relationship, hence improve the accuracy of the prediction.

2.
Due to the lack of knowledge of fracture distribution between wells, the fractures interpreted from all three offset wells were used to predict the stimulation behavior of natural fractures, and it was assumed that a similar fracture distribution would be found in all formations. In reality, the fracture distribution is likely to be different, depending among other things on the lithology and structural location. For example, it is already noticed that there are fewer fractures in the lower part of the reservoir than in the upper part in the D3 well. Intervals with dense fracture networks are more likely to benefit from slickwater treatment compared to formations with no or very sparse fractures. A 3D description of the fracture distribution is always preferred.

3.
Micro-seismic imaging is not available in the study area. No wells are close enough to work as a monitoring well and surface monitoring is also impossible due to the great depth of the reservoir. The lack of microseismic data made it impossible to calibrate the prediction of the shape of SRV.
The main uncertainty in gel frac productivity estimation comes from the propped fracture conductivity estimation. This conductivity is based on proppant testing in the laboratory. The proppant inside fractures involves clogging, crashing and embedment over the production period. There is no analytical method available to model these long-term effects on propped fracture conductivity. An approximate conductivity damage factor has been used in this study to consider these effects.
Although there are still some shortcomings with the workflow, it can assist in the assessment of development concepts and the evaluation of stimulation enhancement options. The anisotropy in the slickwater treatment can be reasonably well-predicted and applied into the production simulation, which provides a more robust prediction than a simple isotropy model. The new workflow can be used in naturally fractured shale gas, tight gas/oil and CBM reservoirs.