Multiscale Flow and Optimal Production Control Techniques in Smart Unconventional Reservoirs

Hydraulic fracturing enables the commercial development of unconventional resources in shales and tight formations. The conductivity and complexity of created fractures are critically dependent on the rheology of fracking fluid and the mechanics properties of rocks. Literatures show that both the rheology of fracturing fluid and fracture propagation dynamics are affected by the temperature of fracturing fluid. Neglecting the temperature transient behaviour may defeat the purpose of fracturing optimization during fracture initiation, propagation, and sand packing. The objective of this paper is to investigate the impact of temperature on fracturing design by studying the transient temperature behaviour across a complex wellbore using numerical modelling by coupling a finite difference heat transfer model with a dynamic fracture propagation model. The study results show that with the injection of cold fracturing fluid, hydraulic fracture propagation is decelerated, and production prediction is thus lessened compared with the case ignoring temperature effect. For multistage fractured wells, fracture geometry enlarges along the fluid flow direction in a horizontal segment. This potentially lowers the cost of hydraulic fracturing designs.


Introduction
The exploration of oil, gas, and geothermal reservoirs has been getting deeper and deeper in these decades [1,2]. Most of these deep reservoirs are identified as unconventional, including extremely tight sandstone and shale reservoirs [3], for which fracking techniques become a must to start and stimulate the production [4]. Once the well is drilled and hydraulic fracking applied, it is important to evaluate the effectiveness of hydraulic fracking in an accurate manner for further predictions. This is done by analysing fracture geometry and conductivity, which are known to be affected by a combination of pressure and temperature [5]. However, it has been commonly found in deep reservoir field practice that the predicted fractured well production exceeds true history data, implying that our knowledge of hydraulic fractures is still lacking in such cases. Studies conducted to investigate pressure impacts with isothermal assumption on these hydraulic fractures [6] may be sufficient in conventional reservoirs which are shallowly buried and have relatively smaller geothermal gradients but do not work out nicely in deep unconventional and geothermal cases which have much higher temperatures.
Research on temperature impacts of fracking is rare compared with pressure. Meyer [7] solved the coupled continuity and energy balance equation for temperature distribution on a fracture surface. Dysart and Whitsitt [8] proposed an analytical model for heat transfer between fluid in fractures and formation. Wheeler [9] related downhole fluid temperature as a function of time, space position, injection condition, and reservoir characteristics with a constant leak-off. Whitsitt and Dysart [10] studied the heat transfer during hydraulic fracturing with a linearly increasing leak-off factor. However, these studies are based on steady-state analytical solutions and do not account for field pumping practice nor formation heterogeneity. Jiang [11] proposed a model for transient heat transfer in vertical wells with porous rock media, but the study does not yield a clear relationship between fluid temperature and fracture propagation, nor did any of these studies consider nonisothermal effect on in situ stress distribution, which also has massive impact on the fracture geometry.
In this paper, a numerical solution to wellbore fluid heat transient transfer is proposed based on Jiang's transient model and used to predict the temperature behaviour during the hydraulic fracturing process. The temperature effects on fluid viscosity and thermal stresses are included into a dynamic fracture propagation model to investigate coupled temperature impacts on fracture growth by analysing the fracture geometry in deep reservoir well fracturing cases.

Materials and Methods
2.1. Transient Heat Transfer Model for the Horizontal Wellbore. Heat transfer during a hydraulic fracturing process occurs when cold fracturing fluid interacts with a hightemperature wellbore and formation through convection and conduction. To include more complexity, in this case, we define a well drilled in a deep reservoir into two separate segments: the vertical and horizontal segments. The structures of segments are shown in Figures 1 and 2, respectively.
For a completed well before fracking, we write the energy balance equation coupled with Fourier's law in the form of where ρ is the fluid density, C p is the fluid heat capacity, v r and v z are flow velocities at radial and vertical directions, respectively, k is the material heat transmissibility, and T is the temperature. And the continuity equation is We apply the spatial and time finite difference to these governing equations, and for different flow regions shown in Figure 1, the heat transfer form is derived into a tridiagonal matrix in the form of The factors A, B, C, and D for different regions are listed in Appendix A.
To connect the vertical and horizontal segments in the case of horizontal wells, a boundary condition between horizontal piece inlet temperature T in and vertical piece bottomhole temperature T bh is applied to the bottom of the vertical segment and start of the horizontal one, which is shown in We use a symmetrical initial temperature distribution around the vertical segment of a wellbore along the radial direction, while heat flux of the upper and lower sides is evaluated separately along the horizontal piece. This is calculated explicitly, and fluid temperature is updated at the start of every time step.

Viscosity Model.
During the hydraulic fracturing process, a desirable viscosity is required for proppant transportation along a fracture. Fluid viscosity is also reported to have significant impact on fracture propagation [12]. Since the fracturing fluid viscosity is temperature sensitive, we use a derived model from Kestin et al. [13] for brine viscosity change with temperature. This viscosity is calculated explicitly with respect to fluid temperature and coupled into the fracture propagation model.
where μ 0 is the viscosity at reference pressure, β is the correlation factor related to temperature and salinity, and p, T, and m are the fluid pressure, temperature, and density, respectively.

Thermal Stress Model.
Stress state changes when nearwellbore field temperature distribution varies, which significantly affects hydraulic fracture propagation. We use a revised model from Gogoi [14] for estimating the thermal stresses around a wellbore in a radial coordinate system and proposed the following model for the in situ stress tensor with respect to transient temperature.

Geofluids
where σ r,th , σ θ,th , and σ z,th are thermal stress components in cylindrical coordinates. T m is the wellbore temperature, and T f ,u is the undisturbed formation temperature. r D is the dimensionless radial distance measured in the unit of wellbore radius, and R D is the dimensionless radius of thermal influence. α l is the thermal linear expansion coefficient, and E is rock Young's modulus. Gogoi's model is calculated explicitly with respect to formation temperature and coupled with the fracture propagation model using Kirsch's solution [15].

Fracture Propagation Model.
We use a fracture geometry model to couple the temperature-influenced parameter for analysing the impact on fracture geometry and production forecast. A few assumptions are made before setting up the model: (1) Constant fluid leak-off coefficient The models of fracture width and length are shown as follows: where t D , L D , and w D are the dimensionless propagation time, propagation fracture half-length, and fracture width, respectively. μ is the fluid viscosity, Q is the fluid injection rate, ν is the Poisson ratio of rock, and G is the bulk modulus.
To account for the change in the stress state, more implementations are made to further calibrate the fracture length and width.
The fracture half-length is a function of pore pressure difference ΔP, shown in The correlation for fracture width is derived and shown in Appendix B. Fracture porosity ∅ F and permeability k F are measured using A workflow of step-by-step modelling and implementation of a hydraulic fracture is shown as follows: (1) Model the temperature behaviour at initial reservoir condition for vertical and horizontal segments The transient temperature profile of formation in the radial direction at the bottomhole is simulated using the heat transfer model, and results are shown in Figure 3. We notice a significant temperature drop around the wellbore during the 300-minute fracking process. The wellbore is gradually cooled down when cold fracking fluid is being injected into it, causing a temperature drop to about 50 degrees Celsius below geothermal temperature. The surrounding formation around the wellbore is also cooled down by fracking fluids. The nearer to the wellbore, the lower the formation temperature becomes.
Brine viscosity is also affected by this process. While the wellbore and surrounding formation are being cooled, the heat is transferred into fracking fluids, heating them up. We Thermal-induced stress during pumping is calculated and superposed into the original stress state profiles, shown in Figure 4. We notice that the induced stress is highest along    4 Geofluids the wellbore (which has the lowest temperature) and decreases along the radial direction. The induced stress inside the formation also changes with time. During the fracking process, the induced stress gradually rises up and finally stays at a steady position at the end of pumping. Coupling all the temperature impacts described in the above sections, we modelled the dynamic growth of a hydraulic fracture with respect to transient heat behaviour. A comparison between the isothermal fracture half-length and the transient heat model fracture half-length is shown in Figure 5. We notice that there has been a significant decrease in fracture half-length, from 162 m to 130 m. The difference, while this may not be this significant in conventional reservoirs with lower reservoir temperature and   On the other hand, the impact on fracture width is much smaller. We only noticed a change of about 0.016 mm in fracture width, which is negligible compared with the total fracture width. This may result from the shape of hydraulic fractures, where the fracture width itself is already small enough.
Simulation has been made by using the different fracture geometries with and without transient heat transfer. We compare the gas production with true field data, and the results are shown in Figure 6.
In this specific case, we can easily infer that the isothermal fracture geometry failed to predict the production at a good accuracy in deep reservoirs. It tends to overestimate the fracture length and thus overestimate the total production. On the other hand, the transient heat model gives a better estimation of production trends, inferring that the method we use is validated.

Horizontal Well Study.
Horizontal wells tend to have a much longer horizontal segment compared with the vertical piece. When we have a long horizontal well, fracking fluids have longer travel time inside the wellbore and cool the formation in a slightly different manner. In this case, the formation temperatures at the inlet of the horizontal segment and the end of the horizontal wellbore are different, causing difference in multistage hydraulic fracturing cases.
We added a horizontal wellbore of 2000 meters to the vertical case and apply a 5-stage hydraulic fracturing process. Stage 1 is at the ending of the horizontal segment while stage 5 is at the inlet position. The simulated temperature profile for the formation near the wellbore along the horizontal segment is shown in Figure 7.
We notice that the formation is being cooled down during pumping, and the end of the horizontal wellbore tends to have a much higher temperature than the inlet. Once we take this into consideration, we can model the multistage averaged fracture half-length at different locations, shown in Figure 8.    6 Geofluids We notice that in this case, the hydraulic fractures at different stages (and thus different distances) have shown their relationship with temperature. The 1 st -stage fractures at the ending of the horizontal wellbore have the longest half-length, and the 5 th -stage fractures at the inlet of the horizontal segment are the shortest. This trend cannot be found if we use isothermal fracture conditions, for which this difference around 15 m will also have impact on production predictions.

Conclusions
Based on the previous analysis and discussions, the following conclusions are drawn: (1) Temperature has shown its impact on different mechanisms during the hydraulic fracking process, including changing the fracking fluid viscosity and formation stress distributions. These changes affect the half-length of hydraulic fractures and therefore have an impact on production predictions. Current fracture models with isothermal assumption may be valid in conventional reservoirs with low reservoir temperature and small geothermal gradients but become inaccurate in deep reservoirs as the temperature impacts become nonnegligible (2) By coupling the mechanisms into numerical simulations, considering these temperature impacts tends to reduce simulated fracture half-length compared with isothermal fracture conditions in deep reservoirs and thus reduce production prediction. This temperature impact is significant and cannot be ignored. On the other hand, the change of fracture width is much smaller and still negligible (3) For multistage hydraulic fractures, the temperature behaviour along the horizontal wellbore tends to reduce the hydraulic fracture half-length with respect to different stage locations. As multistage fracturing usually starts from the ending of horizontal segments and moves along to the inlet, hydraulic fracture halflength also tends to decrease along this process. As the horizontal wellbore is usually much longer than the vertical piece, this impact on fracture length cannot be neglected

Appendix A. Heat Transfer Model Coefficients
The A, B, C, and D factors used in the transient heat model for different regions are listed here. The subscript j is used to infer that these derivations will be used into a numerical system and stands for grid numbers.
For an open-hole completed production well which is meant to be fractured, there are 3 different regions to be analysed: pipeflow, casing and cementing, and formation: For the pipeflow region (subscript 1), ðA:1Þ For the casing and cementing region (subscript 2),  where v stands for fluid velocity, Δt stands for fluid flow time, k is the heat transmissibility, h is Newton's surface heat transfer coefficient, r in , r out , and r ann stand for radius of inlet, outlet, and annulus regions, k eff is the effective rock heat transmissibility in porous media, ρ is the density, and C p is the heat capacity. where y is the superposition parameter in space. Table 1 lists the geological and geomechanical properties of reservoir formation as well as fracking schedule information.  Table 2 lists the thermal properties of formation, completion materials, and fracking fluids.

Data Availability
The data used to support the findings of this study are included in the article.