Schedule for: 18w5085 - Hydraulic Fracturing: Modeling, Simulation, and Experiment

A buffet dinner is served daily between 5:30pm and 7:30pm in the Vistas Dining Room, the top floor of the Sally Borden Building.

Breakfast is served daily between 7 and 9am in the Vistas Dining Room, the top floor of the Sally Borden Building.

A brief introduction to BIRS with important logistical information, technology instruction, and opportunity for participants to ask questions.

In hydraulic fracture stimulation of conventional reservoirs (e.g. tight gas and deep water unconsolidated sands), the use of sophisticated design models is almost indispensable. These hydraulic fractures are typically single fracture treatments and executed from near vertical wellbores. It is well understood that the post-fracture productivity is directly linked to achieving an optimum hydraulic fracture conductivity, which is governed largely by propped fracture length and width. From an engineering and operation execution point-of-view, the goal is to pump into the fracture the desired (large) volume of proppant without encountering pre-mature ‘screen-out’. Therefore, the prediction of fracture geometry and the design of pad volume become critical for propped fracture design. Models are calibrated on-site with dedicated mini-frac tests prior to main propped fracture treatments. Two important calibration parameters are fluid efficiency (leak-off behavior) and minimum in-situ stress (stress profile/contrast). The injection pressure during fracturing, which is readily available, is a valuable source of information and often analyzed and compared with model prediction for fracture diagnostics. In practice, a wide range of models have been employed successfully. However, such considerations do not appear to be important for unconventional resources where multiple fractures are pumped from a long horizontal well. In fact, multi-fracced horizontal well technology has advanced through field trials and experimentation without much help from modelling or understanding of multiple-fracture mechanics. Perhaps one reason is the lower risk of screening-out. This could be due to 1) the extreme low permeability of unconventional shales, which renders the use of high proppant concentration unnecessary and 2), the treatment of multiple fractures in one stage of pumping allows for one or two of the fractures to screen-out without causing an unacceptable rise of pumping pressure. In fact, with the pumping of tens and even a hundred fractures in one horizontal well, the ‘system’ appears to tolerate some ‘non-performing’ fractures without impairing the ultimate production. Conventional wisdom has it that fracture length should be maximized, but in the development of onshore unconventional resources, the horizontal wells are spaced evermore closer to each other, and consequently, the fracture length may not need to be long in order to access the reserves. Operators have successfully fractured and produced from unconventional reservoirs without the use of advanced modelling technology. This begs the questions of what areas of research and model design parameters should we focus on? Can we avoid the ‘details’ while dealing with the ‘big picture’ such as fractures spacing, horizontal well length/direction, the well’s landing depth, and their impact on cost and production? Are research and model development sufficiently guided and tested by field data/observations?

Significant progress in developing ultra-tight unconventional resources has been achieved with horizontal drilling and massive hydraulic fracturing. These technologies enable economic production from sub-microDarcy rocks; however, they introduce significant uncertainty. The wells are drilled from pads and in the subsurface are spaced at about 300-400 m (1000-1300 ft). The wells are then stimulated in stages with variable number of sleeves or perforating clusters, which brings uncertainty about the injection rates into each hydraulic fracture. Poor geological characterization of potential fracture barriers and of natural fractures result in significant uncertainties about the created geometries of hydraulic fractures and distribution of proppant. The in-situ fracture conductivity is also poorly understood. This incomplete list of sources of uncertainty in hydraulic fracturing stimulation of unconventional wells explains a dilemma faced by the operators: to make the development decisions based on physics-based numerical modeling or based on field experience and statistical analysis. The latter becomes a viable alternative in a view of enormous number of producing wells and increasing amount of completion and production data. This presentation reviews outstanding challenges in modeling hydraulic fracturing that need to be addressed to make physics-based numerical modeling relevant to the development of unconventional fields. These challenges can be divided into problems that are understood but difficult to solve or implement (e.g., modeling of the wellbore hydraulics for multi-cluster treatments or numerical efficiency of solving integrated multi-well problems); and those that are not yet understood (proppant transport in realistic non-planar rough fractures in the presence of geological heterogeneities). Addressing these challenges will require advances in numerical modeling, field data acquisition and laboratory experimentation.

The main purpose of this presentation is to give a short review of the existing hydraulic fracture models in Schlumberger, and discuss their specific features, advantages, disadvantages, and modeling challenges. Such models as Planar3D, and UFM, will be presented, and main modeling challenges will be discussed, such as - Height growth modeling in heterogeneous layered media - Interaction with pre-existing natural fractures (crossing criterion) modeling. What we have modeled so far is for weak interfaces. Some of the aspects which need to be better understood: o mineralized fracs - not completely certain how they affect the crossing behavior; o if fracture propagates along the dominant weak plane – how do the smaller fracs or defect affect the fracture bifurcation; o leakoff into the interfaces (NFs). - Interaction between closely spaced hydraulic fractures/branches (stress shadow effect) and stability issues - numerical or real? How to deal with numerical instabilities? - 3D effects and CPU efficiency

Subsurface fluid injection is often followed by observations of an enlarging cloud of microseismicity. The cloud’s diffusive growth is thought to be a direct response to the diffusion of elevated pore fluid pressure reaching pre-stressed faults, triggering small instabilities; the observed high rates of this growth are interpreted to reflect a relatively high permeability of a fractured subsurface [e.g., Shapiro, GJI 1997]. We investigate an alternative mechanism for growing a microseismic cloud: the elastic transfer of stress due to slow, aseismic slip on a subset of the pre-existing faults in this damaged subsurface. We show that the growth of the slipping region of the fault may be self-similar in a diffusive manner. While this slip is driven by fluid injection, we show that, for critically stressed faults, the apparent diffusion of this slow slip may quickly exceed the poroelastically driven diffusion of the elevated pore fluid pressure. We also examine recent field injection experiments providing time series, measured at the borehole, of both fluid pressure as well as the relative displacement of a fault cross-cutting the borehole [Guglielmi et al., 2015]. We couple a hydrogeologic model for fluid flow from the borehole with a model for an expanding shear rupture of the fault. We find that such a model reproduces the observed time history, with a Bayesian inversion providing uncertainties of the model parameters for host rock stiffness and frictional strength, fault zone storage and permeability, as well as the pre-injection stress state. Remarkably, we also find that the inferred rupture front outpaces the region of significant pore pressure increase.

Lunch is served daily between 11:30am and 1:30pm in the Vistas Dining Room, the top floor of the Sally Borden Building.

Meet in the Corbett Hall Lounge for a guided tour of The Banff Centre campus.

Strain-softening/hardening tensile laws, derived from laboratory tests on nano-cantilever beams of kerogen rich shales (KRS) were implemented as user-defined spring models in XSITE, a DEM based lattice code developed by Itasca. The macroscopic tensile strength and toughness properties of the simulated KRS materials were obtained by performing direct, self-similar notched tension tests on microscopic samples with low, medium and high kerogen contents in the lattice code. The notched tests were characterized by two dimensionless numbers, i.e., the ratio of initial crack size over sample width, and the crack resolution in the lattice (the ratio of initial crack size over the lattice resolution). With the first ratio and lattice resolution fixed, parametric studies were performed to generate log-log plots of the critical tensile stress versus the crack resolution in the lattice. The macroscopic tensile strength of the simulated KRS material was estimated from the horizontal plateau obtained at low resolution values in the plots. The toughness corresponding to the condition of LEFM (linear elastic fracture mechanics) was measured from the segment observed at high values of resolution in the logarithm plots where the slope was -½. The toughness values derived from the log-log plots were larger for higher kerogen content, but overall were quite small compared to typical values for shale. Some interesting phenomena were observed in the simulations, e.g., the extent of the process zone near crack tips appeared to remain constant in the tests at nominal low kerogen content; the value of crack resolution beyond which LEFM applies increases as the kerogen content increases (which is consistent with the lower brittleness, or higher plasticity level, observed experimentally at high kerogen content). Using the calibrated macroscopic tensile strength and toughness, fluid injection tests were simulated in meter-scale numerical models at uniform nominal kerogen content. The numerical injection results indicated that, as kerogen content increases, the cluster injection pressure and hydraulic fracture radius increase but the maximum aperture decreases. The behavior is consistent with that expected from the evolution of a penny-shape crack in the viscosity dominated regime.

Natural fractures are widely observed in unconventional oil and gas reservoirs through core samples, image logs, mineback experiments, and outcrop studies. The existence of natural fractures can strongly impact the hydraulic fracture propagation, potentially influencing the effectiveness of reservoir stimulation and hydrocarbon production. Past studies on hydraulic fracture-natural fracture interaction typically assume uniform properties on natural fractures that persist through the entire height of the reservoir/hydraulic fracture, which allow simplification of problems within a two-dimensional (2D) framework. Recent field observations, however, demonstrate that natural fractures can be partially cemented and/or with a height that is less than the reservoir/hydraulic fracture height. These spatially-varied natural fracture properties can influence the morphology of propagation patterns and fracture network differently and require a three-dimensional (3D) consideration. In this study, we present a series of analogue laboratory experiments for hydraulic fractures crossing partially-cemented and/or non-persistent natural fractures. It is observed that a strong enough region on an otherwise weak natural fracture is sufficient to promote crossing. Also, a hydraulic fracture is able to engulf limited-height weak natural fractures and continue propagation after crossing. Guided by experimental observations, an analytical criterion based on linear elastic fracture mechanics is derived to capture the dependence of crossing/no crossing behaviors on spatially-varied natural fracture properties, including the proportion of the cemented region, natural fracture height, and cementation strength. The criterion is further compared with fully-coupled 3D lattice simulations and good agreements are achieved.

Swarming morphologies, that is, those involving multiple aligned members separated by a finite spacing, emerge from systems involving interplay among three fundamental drivers: 1) Alignment: Move in the same direction as neighboring members, 2) Avoidance: Not running into other members, and 3) Attraction: Do not move too far away from other members. As with other systems resulting in swarm-like morphologies, simulation of multiple hydraulic fractures requires a model accounting for the interplay of these fundamental drivers. Specifically, alignment corresponds to the control of the ambient stress field on hydraulic fracture orientation that leads to predominance of certain strike directions. Avoidance drives hydraulic fractures to separate from one another and/or suppress one another’s growth due to the energetic consequences of propagation within the region of elevated compressive stresses surrounding each fracture. Finally, attraction arises from the reduction of viscous energy dissipation associated with splitting the injected fluid among many growing hydraulic fractures rather than just one dominant fracture. When combined, theoretically predicted alignment and emergent spacing in hydraulic fracture swarms matches match with field observations for both hydraulic fractures and naturally occurring dyke swarm analogues. Unfortunately, however, some of the most tempting simplifications, such as neglecting fluid flow or using a two-dimensional modeling domain, result in omitting or fundamentally altering the energetics associated with one or more of the three drivers of hydraulic fracture swarms. As a result, certain simplifications result in a complete loss of model fidelity. On the other hand, reasonably accurate simulations can be obtained from heavily simplified models as long as they preserve the three basic drivers and first principles such as volume balance.

Stress shadowing, a well-known effect that occurs in multi stage hydraulic fracture operations when the hydraulic fractures are placed close to each other, is an important challenge to obtain the highest estimated ultimate recovery (EUR) from a horizontal well. In industry, the most common practices of wellbore completion include three to five perforation clusters (i.e. entry points from the cased wellbore to the formation) per stage. Ideally, each cluster takes the same fluid volume during the hydraulic fracturing operations, leading to uniform stimulation of the reservoir. However, because of stress shadowing, some of the clusters tend to dominate others resulting in an unequal growth of the hydraulic fractures. Motivated by a need to validate and benchmark models used to select perforation spacing, fluid viscosity, injection rate, and so forth that will minimize the negative impacts of stress shadowing, our research focuses on laboratory experiments on the behavior of multiple, simultaneously growing hydraulic fractures. The experimental results show the impact of fracture spacing, fracture height, and number of fractures on multiple fracture growth. We demonstrate qualitative similarity in many respects to existing numerical simulations. However, we also find that certain predicted geometries are apparently less stable than others are when subjected to natural perturbations associated with laboratory experiments.

Numerical modeling is one of the tools that can be used for designing an optimal hydraulic fracturing treatment. One approach is to solve a fully 3D problem of fracture propagation numerically. However, numerical solution of the latter problem is computationally expensive and may preclude one to use it for problems involving optimization or sensitivity analysis. At the same time, there are more specialized models that typically rely on a series of assumptions, but are substantially faster to run. For instance, such models include plane strain, radial, Perkins-Kern-Nordgren (PKN), and pseudo-3D (P3D) hydraulic fractures. The primary aim of this talk is to present a numerical model that extends the aforementioned specialized models for a single fracture into the hydraulic fracturing simulator for multiple cracks. This is done by developing an algorithm for a single fracture, which is then extended to multiple cracks. The numerical algorithm utilizes a fixed mesh approach, in which fracture grows by extension of the tip elements that are eventually split into two parts. Tip element extension utilizes the universal asymptotic solution that originates from the problem of a semi-infinite crack, which includes the effects of toughness, fluid viscosity, and leak-off. The algorithm has been tested against the solution for a plane strain hydraulic fracture in different regimes. In addition, the approach has been extended and tested against enhanced PKN and enhanced P3D models. One of the advantages of the developed model is the fixed mesh methodology, which enabled us to extend the model to multiple fractures that can change their direction of propagation with time. Extension to multiple fractures poses an additional challenging problem of solving for the elastic interaction between the cracks. To address this problem, we use Displacement Discontinuity Method, which has been modified by using elliptical fracture elements. To check accuracy of the developed simulator, its predictions are compared to the reference solution, that is computed using Implicit Level Set Algorithm.

In engineering design for multi-stage HF treatments of horizontal well stimulation, it is ideal to promote simultaneous growth of all fractures in each stage in order to reduce the number of non-producing perforation clusters. While increased attention has been given to studies of multiple HF growth, time dependence is not typically considered as a factor affecting the HF initiation and following growth. A combined experimental and modeling study is carried out to explore the occurrence of the time-dependent initiation of single/multiple hydraulic fracture(s) and their subsequent propagation. By showing the existence of HF initiation at wellbore pressures that are insufficient to induce instantaneous initiation, and explaining that its underlying mechanism is due to the stable growth of the hydraulic fracture under subcritical conditions, this research leads to new insights for promoting more evenly growth of multiple hydraulic fractures in multi-stage HF treatments. Furthermore, our experimental results indicate that the time delay associated with hydraulic fracture initiation can be affected by various factors, such as the fluid viscosity and acidity, and the confining stresses, thereby leading to the practically-relevant outcome that fluid(s) can be chosen in order to promote initiation and growth of multiple hydraulic fractures and/or single hydraulic fractures under conditions where the required wellbore pressure for instantaneous initiation cannot be reached.

A buffet dinner is served daily between 5:30pm and 7:30pm in the Vistas Dining Room, the top floor of the Sally Borden Building.

A computational method is presented for the simulation of fluid-driven fracture propagation in a planar crack using unstructured triangular mesh elements. The crack opening elastic interactions between each element are determined using the displacement discontinuity boundary element method. A moving mesh, incorporating appropriate fracture tip asymptotic representations that depend on the fluid viscosity and the fracture toughness, is advanced and periodically regenerated at the crack front. The method includes special logic to represent intermediate viscosity-toughness crack tip asymptotic behaviour. The approach is validated using the simple geometry of a penny-shaped hydraulic fracture and is applied as well to the case of fracture propagation in a discontinuous stress field.

Hydraulic fracturing is a technique for stimulating oil and gas wells, in which a viscous fluid is injected deep into a rock formation to produce high conductivity channels that facilitate flow of hydrocarbons back to the surface. Even for simple fracture geometries, such as plane strain or axisymmetric fractures, the solution features an interesting behavior due to interplay of physical mechanisms associated with the fluid viscosity, fracture toughness, and fluid leak-off into the formation. In particular, it is known that there are four types of self-similar solutions that correspond to the so-called regimes of propagation. The latter solutions occur for some limiting parameters that correspond to domination of one physical process, such as viscosity or toughness. The global solution, on the other hand, gradually transitions from one regime (or self-similar solution) to another in time. The “structure” of the global solution in the parametric space is investigated for plane strain and radially symmetric fractures. That is, location of the solution relative to the limiting cases is obtained for any problem parameters. Developments are extended to the case of a planar fracture driven by a power-law fluid in an anisotropic (but homogeneous) rock formation. Propagation of multiple closely spaced hydraulic fractures with limited entry design is studied with respect to the regime of propagation. It is found that the fracture shapes evolve from “pancakes” to a “flower” during the transition from the viscosity to the toughness dominated regimes. In the former case, all fractures are mostly radially symmetric and have approximately the same size. At the same time, the fracture “flower” is formed when each fracture has a shape of a petal. There is almost no overlap between the fractures if one observes them from the side. This enables one to influence fracture geometry of multiple fractures in field applications by controlling the regime of propagation.

The geometry and propagation of fluid-driven fractures is determined by a competition between the flow of viscous fluid, the elastic deformation of the solid, and the energy required to create new surfaces through fracturing. To date, much research has focused on the formation of idealised penny-shaped cracks in elastic media [1]. However, the dynamics of fluid-driven fracturing of thin adhered elastica remain unexplored and unobserved, and provide an experimentally accessible and theoretically simpler setting in which to assess the underlying physical processes. We present a theoretical and experimental approach to model a ‘fracture’ produced when fluid is injected from a point source between a solid horizontal plane and an elastic sheet, which is adhered to the plane. Divergence of viscous stresses necessitates the formation of a vapour tip between the fluid front and fracture front. This results in two dynamical regimes of spreading: viscosity dominant spreading controlled by the flow of viscous fluid into the vapour tip, and adhesion dominant spreading controlled by the energy required to fracture the two layers. Constant flux experiments using clear elastic sheets (PDMS) enable new, direct measurements of the vapour tip and confirm the existence of spreading regimes controlled by viscosity and adhesion. We extend this work to consider the possibility of turbulent flow within the body of the fracture and assess the scale of the laminar tip at the fracture front. This analysis identifies the transition from turbulent to laminar control of the spreading, or equivalently the transition from bulk to tip control. These processes primarily feature industrially in the hydraulic fracturing of shale [2], but are also commonplace in nature, from magmatic intrusions in the Earth’s crust [3, 4], to the propagation of cracks at the base of glaciers [5]. [1] D. I. Garagash and E. Detournay, “The Tip region of a Fluid-Driven fracture in an Elastic Medium,” J. Appl. Mech. 67, 183-192 (1999) [2] E. Detournay, “Mechanics of Hydraulic Fractures,” Annu. Rev. Fluid Mech. 48, 311-339 (2016) [3] C. Michaut, “Dynamics of Magmatic Intrusions in the Upper Crust: Theory and Applications to Laccoliths on Earth and the Moon,” J. Geophys. Res. 116, 1-19 (2011) [4] A. M. Rubin, “Propagation of Magma Filled Cracks,” Annu. Rev. Earth Planet. Sci. 23, 287-336 (1995) [5] V. C. Tsai and J. R. Rice, “A Model for Turbulent Hydraulic Fracture and Application to Crack Propagation at Glacier Beds,” J. Geophys. Res. Earth Surf. 115, 1-18 (2010)

In this talk, we will review the implications of the use of high rate slickwater hydraulic fracture treatments as often performed in unconventional gas reservoirs. In particular, due to the very large injection rate used (with rates up to 25 Barrels per minutes in a multistage context when not all the fractures within a stage propagate), the assumption of laminar flow in the fracture may be challenged - at least in the near-wellbore. This is particularly striking if one takes the properties of water to compute a fracture inlet Reynolds number: e.g. thus obtaining inlet Reynolds number up to 4; 000 for a PKN fracture geometry [6]. However, in practice, so-called “friction reducers” are always added in small quantities to the injected water in order to reduce the pressure drop in the wellbore (where the flow is turbulent) and thus minimize the pumping energy required on site (i.e. the numbers of pumping trucks). These friction reducers are high molecular weight (macro molecules) polymers (typically polyacrylamide-based), whose micellar structures completely change the transition from laminar to turbulent flow [5] - making the water “slick”. The effect of the addition of these polymers do saturate at a finite concentration, where the so-called maximum drag reduction asymptote is reached. Such a “saturating” concentration is actually quite low such that it is always targeted in engineering practice. We will present limiting solutions for hydraulic fracture growth in the case of a turbulent maximum drag reduction flow for both the height contained (PKN) and radial fracture geometries in the zero toughness limit. In particular, we will show that most turbulent flow regimes can be recasted as modified power-law fluids as first discussed in [4] for the turbulent rough Gauckler-Manning-Strickler regime. This allows to partly adapt a number of existing solutions for hydraulic fracture growth [3, 1, 2]. We will discuss our results in light of typical operational parameters, highlighting the importance of the drag reduction of slickwater at large Reynolds number. References [1] Adachi, J. I. and Detournay, E. [2002], ‘Self-similar solution of a plane-strain fracture driven by a powerlaw fluid’, International Journal for Numerical and Analytical Methods in Geomechanics 26(6), 579–604. [2] Madyarova, M. and Detournay, E. [2004], Radial Fracture driven by a Power-law Fluid in a Permeable Elastic Rock, Technical report, Schlumberger. Report of UMN to Modeling & Mechanics Group, EAD, Schlumberger. [3] Savitski, A. and Detournay, E. [2002], ‘Propagation of a penny-shaped fluid-driven fracture in an impermeable rock: asymptotic solutions’, International Journal of Solids and Structures 39(26), 6311–6337. [4] Tsai, V. and Rice, J. R. [2010], ‘A model for turbulent hydraulic fracture and application to crack propagation at glacier beds’, J. Geoph. Res. - Earth Surface. [5] Virk, P. S. [1975], ‘Drag reduction fundamentals’, AIChE Journal 21(4), 625–656. [6] Zia, H. and Lecampion, B. [2017], ‘Propagation of a height contained hydraulic fracture in turbulent flow regimes’, International Journal of Solids and Structures 110-111, 265–278.

Meet in foyer of TCPL to participate in the BIRS group photo. The photograph will be taken outdoors, so dress appropriately for the weather. Please don't be late, or you might not be in the official group photo!

Fluid-driven fracture presents an interesting case of crack elasticity and fracture propagation nonlinearly coupled to fluid flow. With the exceptions of a few numerical studies, previous hydraulic fracture modeling efforts have been based on the premise of Linear Elastic Fracture Mechanics (LEFM): specifically, that the damage (aka cohesive) zone associated with the rock breakage near the advancing fracture front is lumped into a singular point, under the tacit assumption that the extent of the cohesive zone is small compared to lengthscales of other physical processes relevant in the HF propagation. The latter include the dissipation in the viscous fluid flow in the fracture channel, of which the fluid lag - a region adjacent to the fracture tip filled with fracturing fluid volatiles and/or infiltrated formation pore fluid - is the extreme manifestation. In this work, we address the validity of the LEFM approach in hydraulic fracturing by considering the solution in the near tip region of a cohesive fracture driven by Newtonian fluid in an impermeable linearelastic rock. First, we show that the solution in general possesses an intricate structure supported by a number of nested lengthscales (a general sentiment for HF), on which different dissipation processes are realized (or are dominant). The latter processes can be cataloged as (1) dissipation in the fracture cohesive zone, “c”, parameterized by the fracture energy Gc (cohesive energy release per unit fracture advance), the peak cohesive stress $\sigma_c$ destroyed by fracturing, and corresponding fracture aperture scale $w_c = G_c/\sigma_c$; (2) the LEFM “reduction” of the cohesive zone process, “$k$”, quantified by $G_c$, but with the cohesive zone replaced by a singularity ($\sigma_c \rightarrow \infty$ and $w_c \rightarrow 0$); (3) viscous fluid dissipation associated with the fluid lag region, “o”, parametrized by an equivalent fracture energy $G_o = \sigma_o w_o$, where $\sigma_o$ is the ${ in \ situ}$ confining stress (signifying the fracturing fluid pressure drop in the lag region from value $\sim\sigma_o$ to near zero) and $w_o$ is the corresponding fracture aperture scale given previously by Garagash and Detournay (2000); and (4) the viscous dissipation along the rest of the fracture (away from the fluid lag), “m” (Desroches et al, 1994). Furthermore, each of the above limiting processes corresponds to distinct solution asymptotes. The HF tip solution structure is bookended by the “c” or “o” (solid or fluid process zones) asymptotes near the tip and by the “m” asymptote away from the tip, while the LEFM “k” asymptote may emerge within the transitional region, as an intermediate asymptote, depending on values of the two governing parameters: the cohesive-to-lag fracture energy ratio $G_c=G_o$ and the cohesive-to-in-situ stress ratio $\sigma_c=\sigma_o$. For typical sets of parameters representative of both field hydraulic fractures and their lab siblings, $G_c/G_o$ is either $\sim 1$ (low-viscosity frac. fluid) or $<<1$ (high viscosity frac. fluid). Under the above conditions, $\sigma_c/\sigma_o>> 1$ is shown to be required for appearance of the LEFM intermediate asymptote near the HF tip. Since $\sigma_c \sim 1$ MPa for most rocks, it can be easily recognized that the latter condition for the relevance of the LEFM to hydraulic fracturing is mostly realized in laboratory experiments conducted under reduced levels of confining stress, and would almost never occur in the field (with the exception of very-near-surface fracturing and/or highly overpressured permeable formations).

Hydraulic fracturing is a process in which fractures are generated in the rock by injection of highly pressurized fluid. Hydraulic fracturing technique together with horizontal drilling allowed to effectively increase oil and gas recovery from low permeable shale formations. The global behavior of a hydraulic fracture is strongly influenced by the processes occurring near the fracture tip, which are related to rock toughness, fluid viscosity and leak-off. The near tip region is modeled as a semi-infinite fracture. The governing equations include elasticity equation, lubrication equation and a propagation criterion. Analytical solutions can be found only for the particular cases of toughness, viscosity and leak-off dominated regimes of propagation. To find a general solution, we employ non-singular formulation to solve the problem numerically. To study the effect of fluid yield stress, the problem of a semi-infinite fracture driven by Herschel-Bulkley fluid is investigated. Numerical results demonstrate that the yield stress influences fracture width solution at larger distances from the tip. At the same time, the solution follows the behavior of a power-law fluid ahead of this zone. Analytical solution for the yield stress dominated regime is obtained and boundaries of its applicability are found. The near tip behavior of a hydraulic fracture can also be strongly affected by proppant - granular material which is mixed with fracturing fluid to prevent fracture from closing after the pressure is removed. Proppant can accumulate near the fracture tip due to settling, bridging, and/or dehydration of the slurry. To investigate the effect of proppant, the problem of a semi-infinite fracture with a localized proppant plug near the tip is analyzed for the case of Newtonian fluid. Fluid filtration through the proppant plug is modeled according to Darcy’s law. Boundaries of the proppant plug are determined by a particle-size dependent bridging criterion and total volume of particles. Proppant causes a noticeable pressure drop over the plug, which in turn leads to fracture widening behind proppant. The effect of proppant can be equivalently represented by a stress barrier solution without proppant. Expressions for magnitude and location of the stress jump are explicitly calculated. Results indicate that such a representation leads to a solution that agrees reasonably well with the numerical solution with proppant.

Although a large number of fluids used in hydraulic fracturing practice exhibit a shear thinning behaviour, little is known on the impact of such a complex fluid rheology on the propagation of a hydraulic fracture. We focus our investigation on the configuration of a semi-infinite hydraulic fracture propagating at a constant velocity in an impermeable linearly elastic material. We allow for the occurrence of a region without fluid of a-priori unknown length at the fracture tip. We use the Carreau rheological model in order to properly account for the shear thinning of fracturing fluid between the low and large shear rates Newtonian limits. We solve this problem numerically combining a Gauss-Chebyshev method for the discretization of the elasticity equation, the quasi-static fracture propagation condition and a finite difference scheme for the width-averaged lubrication flow. This yields in a system of non-linear equations for the fluid pressure in the filled region of the fracture and the extent of the fluid lag region near the fracture tip. We show that for a Carreau rheology, the solution depends on four dimensionless parameters: a dimensionless toughness (function of the fracture velocity, confining stress, material and fluid parameters), a dimensionless transition shear stress (related to both fluid and material behaviour), the fluid shear thinning index and the amplitude of the shear thinning behaviour of the fluid (captured by the ratio of the high and low shear rate viscosities). The solution exhibits a complex structure with up to four distinct asymptotic regions as one moves away from the fracture tip: a region governed by the classical linear elastic fracture mechanics behaviour near the tip, a high shear rate viscosity asymptotic and power-law asymptotic region in the intermediate field and a low shear rate viscosity asymptotic far away from the fracture tip. The occurrence and order of magnitude of the extent of these different viscous asymptotic regions are obtained analytically. Our results also quantify how shear thinning drastically reduces the size of the fluid lag compared to a Newtonian fluid. We also investigate the response obtained with simpler rheological models (powerlaw, Ellis). In most cases, the power-law model does not accurately match the predictions obtained with a Carreau rheology. In the zero lag limit, the Ellis model properly reproduces the results of a Carreau rheology, albeit only for a dimensionless transition shear stress below a critical dimensionless transition shear stress whose expression is given analytically as function of the shear thinning index and magnitude.

This talk investigates the tip region of a hydraulic fracture propagating near a free surface via the related problem of the steady fluid-driven peeling of a thin elastic layer from a rigid substrate. The solution of this problem requires accounting for the existence of a fluid lag, as the pressure singularity that would otherwise exist at the crack tip is incompatible with the underlying linear beam theory governing the deflection of the thin layer. These considerations lead to the formulation of a nonlinear traveling wave problem with a free boundary, which is solved numerically. The scaled solution depends only on one number K, which has the meaning of a dimensionless toughness. The asymptotic viscosity- and toughness-dominated regimes, respectively, corresponding to small and large K, represent the end members of a family of solutions. It is shown that the far-field curvature can be interpreted as an apparent toughness, which is a universal function of K. In the viscosity regime, the apparent toughness does not depend on K, while in the toughness regime, it is equal to K. By noting that the apparent toughness represents an intermediate asymptote for the layer curvature under certain conditions, the obtention of time-dependent solutions for propagating near-surface hydraulic fractures can be greatly simplified. Indeed, any such solutions can be constructed by a matched asymptotics approach, with the outer solution corresponding to a uniformly pressurized fracture and the inner solution to the tip solution derived in this talk.

We discuss the Hydraulic Fracture (HF) model introduced in [1] accounting for the hydraulically induced shear stresses at the crack faces. The model utilizes a general form of the boundary integral operator alongside a revised fracture propagation condition based on the critical value of the energy release rate. The tip asymptotics of the revised model is always consistent with that of the Linear Elastic Fracture Mechanics. We have found that energy release rate criterion takes a more general form in this case and, in fact, plays the role of a natural regulariser in the numerical simulations. As a result, the hydraulically induced tangential tractions may play a significant role in the small toughness and viscosity dominated regimes of crack propagation, while for other regimes the reported results are close to those obtained in the classic model. We have also found that, in case of small toughness or viscosity dominated regimes, the crack redirection angle may change rather significantly for a mixed mode loading [2]. Certain aspects of the recent discussion on the topic [3-5] will be presented and commented. The potential of the revised formulation in tackling some challenges of HF modelling will be demonstrated. References [1] Wrobel, M., Mishuris, G., & Piccolroaz, A. (2017). Energy release rate in hydraulic fracture: Can we neglect an impact of the hydraulically induced shear stress? International Journal of Engineering Science, 111, 28–51. [2] Perkowska, M., Piccolroaz, A., Wrobel, M., & Mishuris, G. (2017). Redirection of a crack driven by viscous fluid. International Journal of Engineering Science, 121, 182–193. [3] Linkov, A. M. (2017). On influence of shear traction on hydraulic fracture propagation. Material Physics and Mechanics, 32, 272–277. [4] Linkov, A. M. (2018). Response to the paper by M. Wrobel, G. Mishuris, A. Piccolroaz “Energy release rate in hydraulic fracture: Can we neglect an impact of the hydraulically induced shear stress?” International Journal of Engineering Science, 127, 217–219. [5] Wrobel, M., Mishuris, G., & Piccolroaz, A. On the impact of tangential traction on the crack surfaces induced by fluid in hydraulic fracture: Response to the letter of A.M. Linkov. Int. J. Eng. Sci. (2018) 127, 217–219

Hydraulic fractures occur miles underground, below complex, layered, heterogeneous rocks, making direct measurements of their dynamics or structure extremely challenging As such, these fractures are typically studied in the lab within blocks of classically brittle materials like glass, PMMA, or rocks that are hydraulically broken with air or fluid (Bunger (2008), Alpern (2012)). Developments in polymer science have shown that heavily cross-linked hydrogels behave nearly identically both qualitatively and quantitatively to these same brittle materials and thus are another good material in which one can study hydraulic fractures (Livne et al (2004)). We have developed a system to study hydraulic fractures within these hydrogels, which have the benefits of highly tunable material properties, being optically clear, and fracture speeds and breakdown pressures 2-3 orders of magnitude lower than PMMA. Using a combination of fast camera photography and laser sheet microscopy, we can study the three dimensional morphology and dynamics of hydraulic fractures at extremely high spatiotemporal fidelity. While the fractures in the gels show excellent agreement with the tip asymptotics outlined in Rice (1968) and Spence and Sharp (1985). However, we also observe instabilities in the propagating fracture front that generate small steps, which leave behind “step lines” that segment an otherwise smooth fracture surface. We show that the density of these lines are the result of increasing mechanical heterogeneity, which we can control in our system, and that at high density, the lines interact resulting in a very rough and uneven fracture surface. This has important practical applications as roughness can be a dominant effect in hydraulic fracture propagation, as well as acting as a nucleation point for the clogging of proppants.

Well pads with multiple horizontal wells are widely used to develop the unconventional tight and shale oil reservoirs, which is driven by both the economic and environmental considerations. Thus, determining the optimal well spacing, placement configuration, and stimulation design is critical to optimizing hydrocarbon production from multi-well pads in unconventional. Recent research efforts have been devoted to maximizing the oil/gas production or the NPV in unconventional reservoirs by utilizing analytical models and reservoir simulations. However, owing to the complexity and computationally expensive simulation of the field-scale problem, the optimization process is mostly restricted to the parametric-sensitivity analysis, where a single variable is varied while others are fixed as constant values. In this work, a new Generalized Differential Evolution (GDE) algorithm has been developed and successfully applied to optimize the well placement as well as fracture parameters of a multi-well pad under constraints. A new well completion economic model based on the field dataset is developed and incorporated into the optimization framework, allowing us to find a practical optimum scenario for the multi-well pad development. A field case in Cardium tight oil reservoir is finally used to demonstrate the successful application of the newly developed optimization framework. It is shown from optimum solutions that the well spacing between 230 and 280 m is considered to be the optimum range for the multi-well pad development in Cardium tight formation. The optimum fracture half-length ranges from 82 to 97 m, and the optimum value of fracture conductivity is between 220 and 240 md⸱m. Under an optimal design of well placements and fracture parameters, the proppant pumped per stage ranges from 15 to 20 tonnages and fracturing fluid injection volume is between 100 and 130 m3 per stage. In summary, the relationship between the overall NPV and total fracture volume is complicated and it is of practical importance to optimize the total fracture volume and strike a balance between the oil production and stimulation cost in order to achieve a higher NPV.

We describe methodologies and robust flow and mechanics algorithms for modeling diffusive fracture network representations in tight formations. These include a priori and a posteriori error estimates for modeling Biot systems, generating natural fracture networks and applying phase field for stimulation .

The computational modeling of the formation and growth of the pressurized and fluid filled fractures in poroelastic media is difficult with complex fracture topologies. Here we study the fracture propagation by approximating lower-dimensional fracture surface using the phase field function. The major advantages of using phase-field modeling for crack propagation are i) it is a fixed-topology approach in which remeshing is avoided, ii) crack propagation and joining path are automatically determined based on energy minimization, and iii) joining and branching of multiple cracks also do not require any additional techniques. Recently, the phase field approach has been widely employed to different applications and developed for various softwares. The two-field displacement phase-field system solves fully-coupled constrained minimization problem due to the crack irreversibility. Here, this constrained optimization problem is handled by using active set strategy. The pressure is obtained by using a diffraction equation where the phase-field variable serves as an indicator function that distinguishes between the fracture and the reservoir. Then the above system is coupled via a fixed-stress iteration. In addition, we couple with transport system for proppant filled fracture by using a power-law fluid system.The numerical discretization in space is based on Galerkin finite elements for displacements and phase-field, and an Enriched Galerkin method is applied for the pressure equation and transport equation in order to obtain local mass conservation. Nonlinear equations are treated with Newton’s method. Predictor-corrector dynamic mesh refinement allows to capture more accurate interface of the fractures with reasonable number for degrees of freedom. In addition, we will discuss how to couple these phase field model to multi scale and optimization problems.

Since their inception in the mid-90's, variational phase-field models of fracture [1] have steadily gained popularity. One of their strengths, the ability to handle complex topologies with unknown crack paths and the interaction between multiple cracks, which is a fundamental requirement for the numerical simulation of hydraulic fracturing in complex situations. Following the technique of [2], crack propagation subject to given pressure $p$ acting along the fracture surfaces of a brittle material occupying a domain $\Omega$ is computed as the minimizer of the energy functional $$\mathcal{E}_\ell (u,\alpha)= \int_{\Omega} \frac{1}{2} (1-\alpha)^2 \mathtt{A} e(u):e(u) \mathrm{d}x - p \int_{\Omega} \alpha \div(u) \mathrm{d}x + \frac{3G_c}{8} \int_{\Omega} \frac{\alpha}{\ell} + \ell | \nabla \alpha |^2 \mathrm{d}x$$ where $u$ denotes the displacement field, $\alpha$ the phase field representing the fracture geometry and $\mathtt{A}$ its Hooke's law. As in [2], in the case of a single pre-existing line or penny-shaped crack in an infinite medium the pressure and volume of fluid recovered from minimizers of (1) can be compared with solutions from the literature [3]. This formalism can also be used to address the issue of hydraulic stimulation of multiple cracks. Symmetry arguments are routinely used to suggest that the propagation of an infinite array of cracks of equal length in an infinite reservoir is possible. Yet a simple stability analysis reveals that this is not the case and that loss of symmetry is always energetically favored [4]. References [1] Bourdin, B., Francfort, G., and Marigo, J.-J., Numerical experiments in revisited brittle fracture. J. Mech. Phys. Solids, (2000), 48(4) 797-826 [2] Bourdin, B., Chukwudozie, C., and Yoshioka, K. (2012). A variational approach to the numerical simulation of hydraulic fracturing. In Proceedings of the 2012 SPE Annual Technical Conference and Exhibition, volume SPE 159154. [3] Sneddon, I. and Lowengrub, M. (1969). Crack problems in the classical theory of elasticity. The SIAM series in Applied Mathematics. John Wiley & Sons. [4] Tanné, E. (2017). Variational phase-field models from brittle to ductile fracture: nucleation and propagation. PhD thesis, Université Paris-Saclay, Ecole Polytechnique.

Since their inception in late 90's, the phase-field models of fracture simulation have steadily gained popularity. One of the appeals is its ability to handle complex topologies with unknown crack paths in relatively coarse meshes as well as multiple-crack interaction, which is a fundamental requirement for the numerical simulation of hydraulic fracturing in complex situations and is technically more difficult to achieve with many other methods. In this talk, we will first describe the construction of a phase-field based coupled hydromechanical reservoir simulator. We will then revisit the problem of a single hydraulic fracture propagating in an infinite impermeable medium in order to validate the computation of fracture width and frac pressure from the phase-field model. Finally we will show how a phase-field description of a system of cracks can be leveraged to model flow in a fractured porous medium and describe the coupling of the flow and mechanics problems, and illustrate the properties of this model through various numerical simulations.

The excursion to the Sulfur Mountain Peak offers an unparalleled 360 degree view of the Canadian Rockies and Bow River Valley in addition to the along-the-ridge broad-walk to the Cosmic Ray Station and the Mountain peak. The Mountain can be ascended by the Banff Gondola (~15 min), or by a hiking trail (1.5-2 hours). https://www.banffjaspercollection.com/attractions/banff-gondola/ The hiking switch-back trail is 5.5 km long and has a ~650 meter elevation gain, and is rated as moderate. The mountain Top offers an observation deck, restaurant/wine-bar, gift shop. If time allows, the Banff Upper Hot Springs are located few minutes walk from the Gondola/mountain-trail base, and provides a way to relax in the naturally hot spring-water and open pool after the mountaineering exercise. http://www.hotsprings.ca/banff-upper-hot-springs (bring you swim suits).

The eXtended Finite Element Method (XFEM) has developed to a standard tool in fracture mechanics. The method enriches the approximation space such that inner-element cracks are considered without loss of accuracy. The hybrid explicit-implicit XFEM uses both, an explicit surface mesh and the implicit level-set method for the description of the crack geometry. The implicit description is needed for the integration in cut elements and it defines where to enrich and how. The explicit description facilitates the (non-planar) crack propagation and provides the basis for solving general transport models on the surface mesh in order to consider for the fluid in Hydraulic Fracturing. One may span the full range from the Reynolds equation, general scalar advection-diffusion equations up to Stokes and Navier-Stokes equations. Because the crack surfaces in hydraulic fracturing may be non-planar, these models have to be extended from the flat case to the situation on curved manifolds. Tangential differential calculus and surface operators play an important role and approximations based on finite elements have to be provided. Future research will show which of these transport models are necessary and sufficient in the field of hydraulic fracturing.

Hydraulic fracturing (HF) is a coupled process involving the simultaneous consideration of both solid deformation and fracturing of the rock mass and the flow of the fracturing fluid. While fully coupled and simple iterative schemes have been shown in the literature, relatively little focus has been placed on the effectiveness and stability of the iterative schemes. This is in contrast to the porous media simulation literature. In the context of the Finite Element Method (FEM) and the eXtended FEM, the most commonly adopted iterative scheme for HF simulation is analogous to the drained split, which in the context of porous media simulation has been shown to be unstable. In this presentation we contrast and compare the HF drained split with a new HF split analogous to the stable undrained split developed for porous media. Through two-dimensional examples of non-planar hydraulic fracture propagation and a benchmark comparison with the KGD model, it will be shown that the new HF undrained split has superior performance in terms of accuracy and load step size, leading to increased computational efficiency.

We describe our approach to simulating curvilinear brittle fractures. Key to our approach is the ability to compute the values of the stress intensity factors around the crack tip with high order of accuracy, in practice fourth order. The practical consequences of this feature are that (a) converged crack paths can be obtained with relatively coarse meshes and, more importantly, (b) it is not necessary to refine the mesh around the tip at each crack step, except perhaps, around high curvature regions of the crack paths. The ability to compute accurate stress intensity factors relies on two novel developments in my group: (a) the use of Universal Meshes to deform an underlying mesh so that it precisely matches the geometry of the fracture as it evolves, and (b) the computation of high-order solutions to elasticity problems with singularities. We will briefly illustrate this method through the simulation of the propagation of thermally driven cracks. We will then focus on the application of the Universal Meshes to simulate hydraulic fractures with lag in two-dimensions. The key advantage here is that the Universal Mesh provides a mesh of good quality on the crack faces for the computation of the lubrication equations. We are still working on the method for the computation of high-order solutions in the presence of the fluid, so we only obtain first order convergence of the stress intensity factors in this case.

Hydraulic fracturing in naturally fractured rocks often leads to the creation of a stimulated zone of enhanced permeability in which the target formation is irreversibly deformed through shear dilation of natural fractures, plastic deformation, and induced bulk damage. The current dominant modeling approach - explicitly accounting for each fracture with microscale resolution of the fracture network (e.g., discrete fracture network or distinct element method) is a computationally expensive and complex task. There also remain large uncertainties with respect to natural fracture distribution and reservoir parameters. Addressing these issues leads to identification of the need to develop up-scaled continuum model that is able to, in an average sense, capture the irreversible behavior of naturally fractured rock masses. We present a novel mathematical approach with the goal of simulating the evolution of the Stimulated Rock Volume (SRV) in a 2D/3D geomechanical model. This is achieved by introducing a homogenized non-local poro-elastic-plastic continuum zone for the stimulated region, described by an internal characteristic length scale. The up-scaled mechanism of fracturing and deformation is described by a non-local Drucker-Prager model coupled to a Biot poroelastic medium, and implemented within a standard Galerkin Finite Element Method framework. We first quantify the evolution of the SRV and pressure change in the reservoir for a typical example of hydraulic fracture stimulation in a tight formation. After the creation of a sufficiently large SRV, the well is shut-in for an extended period of time and the wellbore pressure is allowed to fall-off. The analysis of post-shut pressure curves confirms the existence of the well-known flow regimes- storage and bi-linear flow- characteristic of the simple bi-wing hydraulic fractures in homogenous rocks. Using the existing analytical solution for the finite conductive fracture, the flow capacity of the stimulated zone is calculated and correlated to the size of the stimulated zone through the non-local length scale. The performance of the developed methodology is tested by considering examples of 2D and 3D SRV calculations. For each example, the stimulated zone and the fluid pressure in relation to the local in-situ stress field are quantified. The influence of reservoir complexities, such as sedimentary layering, complex initial in-situ stress field, and wellbore effects on the evolution of the SRV and fluid pressure will be discussed.

The inversion of remote tilt measurements to determine the geometry of an evolving hydraulic fracture (HF) is a classically ill-posed problem due to the fact that the elliptic PDE that governs the behaviour of the displacement gradient field rapidly smooths the geometric details of the fracture with distance. Indeed, it is possible to show that it is only possible to obtain reasonable estimates of the first few moments of the crack opening displacement field from such measurements. On the other hand, numerical and analytic models of evolving HF cannot be expected to provide a completely accurate prediction of evolving HF geometries taking place in the field due to the large number of uncertainties in the data and the un-modeled dynamics due to physical processes that have, of necessity, been ignored. Our approach is to feed the time series of tilt data as input to the implicit level set algorithm (ILSA) model for an evolving HF via the Extended Kalman Filter (EKF). Form an inversion point of view, the dynamics from the coupled ILSA model enables the tilt data snapshots in the time series to be connected, where in previous inversion algorithms these data were regarded as independent. We illustrate the EKF-ILSA algorithm using numerical experiments for planar HF propagating in 3D elastic media. By varying the confining stress field, synthetic tiltmeter data are generated that result in substantial changes to the geometry of the evolving HF. The ILSA model is assumed to have no knowledge of this confining stress field except for feedback from the tiltmeters via the Extended Kalman Filter. Indeed, without this feedback the ILSA HF model would propagate with radial symmetry. We compare the EKF-ILSA estimates of the fracture geometry and width with those of the HF used to generate the synthetic data with and without Gaussian noise. We also present results in which the algorithm is tested on real field data from a mining situation in which HF have been deliberately generated to enhance caving in longwall coal mining. The model is able to detect asymmetry in the growth of the HF, which is corroborated by measured intersection times of the HF with monitoring boreholes.

Models of hydraulic fractures in conventional reservoirs assume that Carter’s leak-off law —the leak-off rate is proportional to the inverse of the square root of the time elapsed since exposure to fracturing fluid — is applicable. The validity of Carter’s leak-off law stems from the cake-building properties of the fracturing fluid. In some situations where water is essentially the fracturing fluid, Carter’s leak-off law can also be justified (through an reinterpretation of the leak-off coefficient) as an early-time solution of the diffusion equation. However, in water flooding operations of very permeable reservoirs, the fracture propagates in a region where the pore pressure perturbations caused by injection of water is quasi-stationary. The talk will present the construction of a new class of solutions for hydraulic fractures propagating under these asymptotic conditions. We will first present a KGD-type model of an hydraulic fracture created by injecting fluid in weak, poorly consolidated rocks. By further assuming a “small” or negligible toughness (with the consequence that the crack aperture is “small”), we prove that the system is characterized by two asymptotic fracture propagation regimes: rock-flow dominated at small time and fracture-flow dominated at large time. The timescale that legislates the transition between the small and large time asymptotic regimes is shown to be a strongly nonlinear function of a dimensionless injection rate. The rock-flow dominated regime is characterized by an increasing injection pressure while the fracture-flow dominated regime is associated with an injection pressure decreasing with time. The peak injection pressure takes place during the transition between the two regimes. The KGD model collapses, however, when the total crack length (two wings) becomes larger than the thickness of the reservoir layer, assumed to be bounded by impermeable strata. The changing geometric ratio of the constant crack height over its length affects the fracture compliance, i.e. the relationship between fracture aperture and net pressure. As this ratio decreases, the non-local elastic interaction characteristic of the KGD model progressively vanishes and for ratio approximately larger than 5 the compliance becomes essentially local as in the PKN model. It will be shown that evolution of the fracture from a KGD to a PKN mode causes an unexpected reversal of the regime of propagation, with the rock-flow dominated regime in the PKN geometry becoming the long-term solution.

The model and algorithm for the numerical solution of the three-dimensional problem of hydraulic fracture initiation and further propagation will be presented. The model is fully coupled and takes into account three important processes: elastic deformation of the rock, fluid flow in the fracture, and its further propagation in the rock. The mathematical model consists of three groups of equations. Each of them responses for one process defined above. The elasticity equations are solved by the dual boundary element method (DBEM), the lubrication equations by the finite element method (FEM) improved by simple conservative correction. This correction allows us to preserve the total volume of injected fluid on the discrete level. The fracture propagation criterion gives the system of non-linear equations, which is solved by special modification of relaxation method. In the early stages of propagation we need to explicitly consider the fluid lag, which in general varies along fracture front. It essentially increases needed computational resources. One of the ways to overcome this challenge is to use any approximation of behavior of variables near the fracture front (tip). Are the already developed asymptotic solution applicable here? The results obtained by the model include the initiation pressure for the real configuration of perforated well, the shape of the fracture, its position and orientation, as well as the possibility of reorientation and the size of the domain where it reorients. The cementing and casing of the well can be taken into account. From the point of oilfield engineer's view, the model can be useful in the understanding of the early stage of hydraulic fracturing when there are many stop cases that sometimes lead to an unsuccessful hydraulic fracturing.

In the talk, we present our recent developments of the modelling of a planar hydraulic fracture in an inhomogeneous reservoir. The hierarchy of models includes the Enhanced Pseudo 3D (EP3D) model [1] coupled with the proppant transport, the Planar 3D model under the modification of the Implicit Level Set Algorithm (ILSA) [2], and the Planar 3D Biot model that accounts for the effect of poroelasticity [3]. In the EP3D model, the three-dimensional fracture is modeled in terms of quantities that are averaged along the fracture’s height dimension. This model is computationally fast, but is limited to a certain shape of the fracture, and takes into account the layered structure of the reservoir only in terms of the confining stress difference. This model is used for the development of the coupling procedure with the proppant transport module. For the latter we use the one-speed transport model with the effective viscosity varying with particle concentration. We demonstrate effects of the Saffman-Taylor instability development, proppant bridging, breakage of the proppant plug due to the displacement instability. The Planar 3D model describes the fracture development in a layered reservoir, where the effects of the poroelasticity are neglected. The model is implemented using a modification of the ILSA approach [2]. In particular, the model accounts for the inhomogeneity of the elastic properties of the reservoir (only the layered structure of the reservoir is allowed), and is able to simulate cases of fracture local closure due to the fluid re-distribution and/or leak-off. The most advanced Planar 3D Biot model describes the fully coupled interaction of the stresses with the fluid filtration. The numerical model is implemented using the Finite Element Method and allows us to account for an arbitrary inhomogeneity of all physical characteristics of the reservoir. In particular, we show that in the case of the layered structure of the formation, where the layers differ only by permeability, the fracture propagation can demonstrate counter-intuitive non-monotonical behavior. Both, Planar 3D and Planar 3D Biot models, are thoroughly verified by comparison with fast approximate solution for the radial fracture [4], and matched with the existing experimental data [5]. Finally, we present results of the multi-parameter and multi-objective optimization for the Net Present Value and the Fracture Production depending on the applied flow rate, volume of fluid and proppant, and other characteristics of the fracturing process. The optimization is based on the fast algorithms for estimation of fracture’s characteristics [6], on the module for computation of the production of the multiply-fractured wellbore, and on application of genetic algorithms [7] for construction of the Pareto front in the space of objective functions. The work was supported by the Ministry of Science and Education of Russian Federation (grant 2016-220-05-2642). Literature 1. E.V. Dontsov, A. P. Peirce. (2015) An enhanced Pseudo-3D model for hydraulic fracturing accounting for viscous height growth, non-local elasticity, and lateral toughness, Eng. Fract. Mech., V. 142, p 116-139 2. E.V. Dontsov, A.P. Peirce. (2017) A multiscale Implicit Level Set Algorithm (ILSA) to model hydraulic fracture propagation incorporating combined viscous, toughness, and leak-off asymptotics. Comput. Methods Appl. Mech. Engrg. V. 313. P. 53–84. 3. A.N. Baykin, S.V. Golovin. (2016) Modelling of hydraulic fracture development in inhomogeneous poroelastic medium. J. Phys.: Conf. Ser., V. 722. 012003 4. Dontsov, E.V. (2016) An approximate solution for a penny-shaped hydraulic fracture that accounts for fracture toughness, fluid viscosity, and leak-off. R.Soc. Open Sci. V. 3. P. 160737. 5. R. Wu, A.P. Bunger, R.G. Jeffrey, E. Siebrits. (2008) A comparison of numerical and experimental results of hydraulic fracture growth into a zone of lower confining stress. ARMA 08-267 6. Dontsov, E. V. (2016) An approximate solution for a penny-shaped hydraulic fracture that accounts for fracture toughness, fluid viscosity and leak-off. Royal Society open science V.3 N.12 P. 160737. 7. Deb, K. (2001) Multi-objective optimization using evolutionary algorithms. John Wiley & Sons.

A fast multipole method (FMM) is used to decrease the computational time of a fully coupled poroelastic hydraulic fracture model with a controllable effect on its accuracy. The hydraulic fracture model is based on the fully poroelastic formulation of the displacement discontinuity method (DDM) which is a special formulation of boundary element method (BEM). DDM is a powerful and efficient method for problems involving fractures. However, this method becomes slow as the number of spatial elements increases, or necessary details such as poroelasticity, that makes the solution history-dependent, are added to the model. FMM is a technique to expedite matrix-vector multiplications within a controllable error without forming the matrix explicitly. Several examples are provided to show the efficiency of the proposed approach in problems with large degrees of freedom (in time and space). Examples include hydraulic fracturing of a horizontal well and randomly distributed pressurized fractures at different orientations with respect to horizontal stresses. The results are compared to the conventional DDM in terms of computational processing time and accuracy. It is demonstrated that FMM may decrease the computation time by up to 70 times with a negligible error. The solution of tip displacements using both methods are then used to compare the computation of stress intensity factors (SIF) in mode I and II, which are needed for fracture propagation. The error of SIF calculation using the proposed modification was also found to be negligible. Consequently, this method will not affect the estimation of the fracture propagation direction. Accordingly, the proposed algorithm may be used for fracture propagation studies while substantially reducing the processing time.

The aim of this work is to analyse the nonlinear response of a porous thick-walled hollow sphere subjected to inner fluid pressure. In particular, we aim to find out how Biot coefficient and Poisson's ratio of the material influence the fracture process. Spherical symmetry is assumed so that the hydromechanical fracture problem is expressed by means of an ordinary differential equation of the radial displacement with respect to the radius of the sphere. For the elastic response, we derive the analytical hydro-mechanical solution based on the work in [2], which is a hydro-mechanical extension of the well known mechanical case described for instance in [1]. It is assumed that the material is fully saturated and that the application of the fluid pressure changes so slowly that a steady state always exists. For the extension to fracture, the crack openings are smeared out into an inelastic strain, which is used in a one-dimensional damage model for the stiffness reduction. The resulting nonlinear ordinary differential equation, which is presented here for the general case of non-zero Poisson's ratio, is solved numerically by means of a finite difference scheme. A sensitivity study reveals that both Biot coefficient and Poisson's ratio have a very strong influence on the hydromechanical response of the thick-walled hollow sphere. References [1] S. P. Timoshenko, J. N. Goodier, Theory of elasticity. McGraw-Hill, 1987. [2] P. Grassl and C. Fahy and D. Gallipoli and S. J. Wheeler, On a 2D hydro-mechanical lattice approach for modelling hydraulic fracture. Journal of the Mechanics and Physics of Solids, Vol. 75, pp. 104{ 118, 2015. 1

Fracture surface roughness can play a significant role in enhanced oil recovery, due to its impact on the fluid dynamics within the thin fracture. In theory, fracture surfaces should be smooth; however, undulations and crinkles of the fracture front resulting from heterogeneities and dynamic instability produce a complex fracture surface. These effects are fast, multi scale and generically three dimensional, and as such, are intractable both experimentally and computationally. Experimentally, these difficulties are mediated by studying fracture dynamics in brittle hydrogels. In these transparent materials, the fracture dynamics are slow and can be visualized. We combine high speed photography and laser sheet microscopy and directly observe in full 3D how roughness is dynamically generated by the fracture front. Specifically, I will discuss the behavior of small step-like discontinuities that form, propagate, and interact with each other along the advancing front. The interaction of these step lines can create significant relief along the fracture front and can even result in fracture propagation on separate parallel planes.

5-day workshop participants are welcome to use BIRS facilities (BIRS Coffee Lounge, TCPL and Reading Room) until 3 pm on Friday, although participants are still required to checkout of the guest rooms by 12 noon.