Abstract:
Embodiments herein relate to a method for hydraulic fracturing a subterranean formation traversed by a wellbore including characterizing the formation using measured properties of the formation, including mechanical properties of geological interfaces, identifying a formation fracture height wherein the identifying comprises calculating a contact of a hydraulic fracture surface with geological interfaces, and fracturing the formation wherein a fluid viscosity or a fluid flow rate or both are selected using the calculating. Embodiments herein also relate to a method for hydraulic fracturing a subterranean formation traversed by a wellbore including measuring the formation comprising mechanical properties of geological interfaces, characterizing the formation using the measurements, calculating a formation fracture height using the formation characterization, calculating an optimum fracture height using the measurements, and comparing the optimum fracture height to the formation fracture height.

Description:
RELATED APPLICATION INFORMATION 
       [0001]    This application claims the benefit of U.S. Provisional Application No. 62/008,082 filed Jun. 5, 2014, which is incorporated herein in its entirety. 
     
    
     FIELD 
       [0002]    This relates to the field of geomechanics and hydraulic fracture mechanics. This relates to oil-and-gas reservoir stimulation, performed by hydraulic fracturing of rock from the wellbore, including providing a technique to predict hydraulic fracture height growth in the rock affected by pre-existing weak mechanical horizontal interfaces such as bedding planes, lamination interfaces, slickensides, and others. 
       BACKGROUND 
       [0003]    For context, we demonstrate the results of two fracture propagation modeling cases with different structure of rock interfaces with respect to the horizontal wellbore. In both examples, one hydraulic fracture is initiated at the horizontal wellbore and propagates in vertical and horizontal directions. The rock properties and in-situ stresses are the same in different layers dividing by the prescribed interfaces for both presented examples. The interfaces are cohesionless but frictional planes of weakness. 
         [0000]    Case of Symmetrical Interfaces with Respect to Wellbore 
         [0004]    In the first example, horizontal interfaces are located symmetrically with respect to the horizontal wellbore. Hydraulic fracture initiated and propagates across these interfaces as well as along them in the horizontal direction, as shown in the  FIG. 1 .  FIG. 1  shows a hydraulic fracture propagating from the horizontal wellbore in the case of symmetrical placement of horizontal interfaces with respect to wellbore. 
         [0005]    Propagating of both vertical tips of hydraulic fracture across the interfaces is relatively slow because of continuous stops at each if these interfaces. At the same time, lateral tips of the hydraulic fracture propagate without interaction with interfaces (parallel to them). As a result, the length of hydraulic fracture appears to be much longer than its height ( FIG. 2 ). 
         [0006]      FIG. 2  shows an upper, lower, and lateral fracture tip propagation with time of fluid injection (upper graph), and corresponding pressure response at the fracture inlet (lower graph) for symmetrical placement of the interfaces. 
         [0000]    Case of Asymmetrical Interfaces with Respect to Wellbore 
         [0007]    In the second modeling case, cohesionless horizontal interfaces are positioned asymmetrically with respect to the wellbore. Number of interfaces below the wellbore is less than that above the wellbore (see  FIG. 3 ). The pumping schedule, the spacing between the interfaces, and all other parameters of the rock and fracture remain the same, as in the first example.  FIG. 3  shows hydraulic fracture propagating from the horizontal wellbore in the case of asymmetrical placement of horizontal interfaces with respect to wellbore. 
         [0008]    Modeling shows that in this case after crossing two interfaces below the wellbore, the hydraulic fracture will be completely stopped at one of the upper interfaces while freely propagates downward ( FIG. 4 ).  FIG. 4  illustrates an upper, lower, and lateral fracture tip propagation with time of fluid injection (upper graph), and corresponding pressure response at the fracture inlet (lower graph) for asymmetrical placement of the interfaces. 
         [0009]    These two examples indicate that the preliminary measurement of the weakness planes in rock and adequate modeling of fracture propagation in a layered formation are needed to identify fracture height containment in a layered rock adequately. And oppositely, missing the information about the heterogeneous profile of the rock strength in the vertical direction and prominent interfaces can result in wrong results in prediction of the fracture height containment conditioned by interaction of the hydraulic fracture with weakness planes. 
         [0010]    Hydraulic fracturing used for the purpose of reservoir stimulation typically aims at propagating sufficiently long fractures in a reservoir. The fracture length can be as large as several hundred meters in horizontal direction. With such fracture extent the layered rock structure reveals severe heterogeneity vertically. Depending of the rock type, sedimentary laminations or beddings can have thickness in the range of millimeters to meters. Unequal variation of rock properties in vertical and horizontal directions results in noticeable restriction of the fracture height growth with respect to lateral fracture propagation. Since the beginning of fracturing era attention to the hydraulic fracture height containment was always recognized. 
         [0011]    Subsurface three-dimensional propagation of hydraulic fractures (hereafter HF) typically implies simultaneous fracture growth in horizontal and vertical directions. Typical horizontal HF extent during field treatments varies from tens to hundreds meters along the intended formation layer. As opposed to that, vertical fracture extent appears much shorter in size because of large contrast of rock properties and tectonic stresses, as well as pre-existing horizontal bedding and lamination interfaces. There are several recognized mechanisms controlling the vertical HF growth (upward or downward) in geologic formations: (1) minimum horizontal stress variation as a function of depth (hereafter called “stress contrast” or “mechanism 1”), (2) elastic moduli contrast between adjacent and different lithological layers (hereafter called “elasticity contrast” or “mechanism 2”), and (3) weak mechanical interface between similar or different lithological layers (hereafter called “weak interface” or “mechanism 3”). A “weak mechanical interface” or “weak interface” or “plane of weakness” refers to any mechanical discontinuity that has low bonding strength (shear, tensile, stress intensity, friction) with respect to the strength of the rock matrix. A weak interface represents a potential barrier for fracture propagation as follows: when the HF reaches the weak interface, it creates a slip zone near the contact as shown by both analytical and numerical studies. Slip near the contact zone can arrest fracture propagation and lead to extensive fluid infiltration or even hydraulic opening of the interface by forming so called T-shape fractures. Such T-shape fractures have been repeatedly observed in various mineback observations in coal bed formations. 
         [0012]    Nowadays, the “stress contrast” mechanism is the main used in most HF modeling codes to control vertical height growth, both for pseudo3D and planar3D models. The “elastic contrast” mechanism is usually not explicitly modeled in most HF modeling codes, but is in some way addressed by the “stress contrast” mechanism as vertical stress profile of minimum horizontal stress are often derived from a calibrated poroelastic model and overburden stress profile (isotropic and transverse isotropy can be treated) that depends on the elasticity of the formation. The “weak interface” mechanism has drawn less attention in the hydraulic fracturing community up to date, though it has been well recognized from field fracturing jobs and discussed in literature as far back as the 1980s. This lack of interest may have been caused by the lack of characterization of the location of the weak interfaces in deep formations and/or the lack of measurements of their mechanical properties (shear and tensile strength, fracture toughness, friction coefficient and permeability). At the same time the “weak interface” mechanism is one of the only of the above mechanisms that can completely stop the HF from further propagating upward or downward in formations. The main reasons for fracture tip termination at weak interfaces are the interface slippage, pressurization by penetrated fracturing fluid, or even mechanical opening of the interface. In contrast, the first two mechanisms may only temporarily stop the HF until the net pressure is increased in the HF up to a threshold level that will allow the HF to further propagate. The “weak interface” containment mechanism may be more important than “stress” or “elastic contrast” mechanisms and may be the reason why HF are often well contained in vertical extent despite apparent absence of any observed “stress” or “elastic contrast.” In any event, more effective methods for formation characterization, existing fracture influence on fracture development, and characterization of fracture generation are needed. 
     
    
     
       FIGURES 
         [0013]      FIG. 1  shows a hydraulic fracture propagating from the horizontal wellbore in the case of symmetrical placement of horizontal interfaces with respect to wellbore. 
           [0014]      FIG. 2 . Upper, lower and lateral fracture tip propagation with time of fluid injection (upper graph), and corresponding pressure response at the fracture inlet (lower graph) for symmetrical placement of the interfaces. 
           [0015]      FIG. 3 . Hydraulic fracture propagating from the horizontal wellbore in the case of asymmetrical placement of horizontal interfaces with respect to wellbore. 
           [0016]      FIG. 4  includes upper, lower and lateral fracture tip propagation with time of fluid injection (upper graph), and corresponding pressure response at the fracture inlet (lower graph) for asymmetrical placement of the interfaces. 
           [0017]      FIG. 5  is a schematic drawing of a vertical hydraulic fracture (HF) growth in a subterranean layered rock with horizontal interfaces. 
           [0018]      FIG. 6  is a flow chart listing the information that may be used for an embodiment herein. 
           [0019]      FIG. 7  provides examples of stages for 3D frac propagation across weak planes. 
           [0020]      FIG. 8  is a flow chart of methods for an embodiment. 
           [0021]      FIG. 9  is a flow chart of a component of a method for an embodiment. 
           [0022]      FIG. 10  depicts an embodiment of an algorithm of the HF simulator ( 200 ) workflow from the beginning of the fracturing job t0 up to the end T. 
           [0023]      FIG. 11  illustrates a horizontal interface crossed by the vertical hydraulic fracture (top), and schematic distribution of the percolated fluid pressure along the interface (bottom). 
           [0024]      FIG. 12  provides a profile of fluid pressure along the interface for the “in-slip” (top) and “out-of-slip” (bottom) regimes of percolation. 
           [0025]      FIG. 13  is a series of schematic diagrams to show a hydraulic fracture propagating upward and downward in plane-strain geometry (vertical cross-section). 
           [0026]      FIG. 14  is a plot that shows the injected, fracture and leaked-off fluid volumes (top), net pressure (middle), and hydraulic fracture halfheight (bottom) during the whole cycle of fluid injection into the fracture. 
           [0027]      FIG. 15  is a two-sided contact of a vertically growing fracture and weak horizontal interfaces (left), interface activation, and fracture tip blunting as a result of the contact with the interfaces (right) 
           [0028]      FIG. 16  provides profiles of the vertical fracture opening (left) at the contact with two cohesionless interfaces and normalized fracture volume versus stress ratio (right). 
           [0029]      FIG. 17  includes the maximum tensile stress component generated on the opposite side of the cohesionless (left) and cohesional interface with κ IIC =1 (right). 
           [0030]      FIG. 18  shows fracture tip propagation (top) and inlet pressure decline (bottom) in the case of an elliptical fracture with Newtonian fluid with viscosity of 1 cP (left) and 10000 cP (right), respectively 
           [0031]      FIG. 19  is a flow chart of a component of a method for an embodiment (solver for hydraulic fracture tip propagation in the absence of interfaces). 
           [0032]      FIG. 20  is a flow chart of a component of a method for an embodiment (sub-component of the above: a coupled solid-fluid solver for hydraulic fracture with given fracture tip position). 
           [0033]      FIG. 21  is a flow chart of outputs of an embodiment of a method. 
       
    
    
     SUMMARY 
       [0034]    Embodiments herein relate to a method for hydraulic fracturing a subterranean formation traversed by a wellbore including characterizing the formation using measured properties of the formation, including mechanical properties of geological interfaces, identifying a formation fracture height wherein the identifying comprises calculating a contact of a hydraulic fracture surface with geological interfaces, and fracturing the formation wherein a fluid viscosity or a fluid flow rate or both are selected using the calculating. Embodiments herein also relate to a method for hydraulic fracturing a subterranean formation traversed by a wellbore including measuring the formation comprising mechanical properties of geological interfaces, characterizing the formation using the measurements, calculating a formation fracture height using the formation characterization, calculating an optimum fracture height using the measurements, and comparing the optimum fracture height to the formation fracture height. 
       DETAILED DESCRIPTION 
       [0035]    Herein, we provide an approach to predict hydraulic fracture height growth in rocks having laminated structure. This method includes (i) a preliminary vertical characterization of the bulk rock mechanical properties, the mechanical discontinuities and in-situ stresses, and (ii) running the computational model of 3D or pseudo-3D hydraulic fracture propagation in the given layered rock formation and taking into account the interaction with the given weak mechanical and/or permeable horizontal interfaces. Methods herein for rock characterization and advanced fracture simulation produce a more accurate prediction of a fracture height growth, fracturing fluid leak-off along weak interfaces, forming T-shaped fracture contacts with horizontal interfaces, and switching from vertical orientation of the fracture to a horizontal one. 
         [0036]    3 mechanisms that control height growth are described in more detail below.
       1. Mechanism 1 (conventional): minimum horizontal stress variation as a function of depth called “stress contrast”   2. Mechanism 2 (conventional): elastic moduli contrast between adjacent and different lithological layers called “elasticity contrast”   3. Mechanism 3 (most important, it is the novelty of this application): weak mechanical interface between similar or different lithological layers called “weak interface”   a. Sub-mechanism 3a: elastic interaction, crossing criterion and re-initiation past-interface   b. Sub-mechanism 3b: enhanced leak-off of the fracturing fluid into the interface       
 
         [0042]    Characterization of Vertical Rock Texture 
         [0043]    In order to make the prediction of fracture height growth precise, information about rock properties, its mechanical discontinuities, and in-situ stresses is required. Information about rock comprises the detailed vertical distribution of mechanical properties of the rock mass, including variation of rock strength, in terms of, for example, tensile strength, compressive strength (e.g. uniaxial confined strength or UCS) and fracture toughness, which should provide information about placement of weakness planes in rock with elastic properties (e.g. Young modulus and Poisson&#39;s ratio). Measurement of rock stresses should bring information about the vertical stress and the minimum horizontal stress in the normal stress conditions, where vertical stress component is the largest compressive stress component (or strike-slip conditions where the vertical stress is the intermediate compressive stress component). 
         [0044]    There are available rock property characterization tools that can be used for mechanical rock property measurement. These are Sonic Scanner, and image logs (e.g. REW: FMI, UBI; OBMI; e.g. LWD: MicroScope, geoVISION, EcoScope, PathFinder Density Imager), which can give information about elastic properties and locations of pre-existing interfaces. If coring is available, in the lab test one can perform heterogeneous rock analysis (HRA) on cores extracted from this rock mass, and scratch test, which provides information about statistical distribution of weakness planes on a core scale and their properties (tensile and compressive strength, fracture toughness). 
         [0045]    In summary, the input properties to be characterized are: 
         [0046]    Density (i.e. inverse of spacing) and orientation (mainly horizontal) of weak interfaces as a function of depth 
         [0047]    Mechanical and hydraulic properties of the weak interfaces (respectively, friction, cohesion, tensile strength, and toughness, and permeability and filling) 
         [0048]    Vertical stress (Sv) as a function of depth 
         [0049]    Minimum horizontal stress (Sh) as a function of depth 
         [0050]    Elasticity of bulk rock (e.g. Young Moduli and Poisson Ratio) as a function of depth 
         [0051]    Chart 1 provides an inventory of data sources and model parameters for a given type of rock and reservoir. SONICSCANNER™ and ISOLATION SCANNER™ tools are commercially available from Schlumberger Technology Corporation of Sugar Land, Tex. 
         [0000]    
       
         
               
               
               
             
           
               
                   
               
               
                   
                 Refers 
                   
               
               
                 Model parameter 
                 to 
                 Potential data source 
               
               
                   
               
             
             
               
                 Vertical profile of rock 
                   
                 High Res petrophysics, 
               
               
                 layers and interfaces 
                   
                 image and sonic logs 
               
               
                 E′ - Young modulus 
                 LAYER 
                 Sonic Scanner or Isolation 
               
               
                 (plain-strain) 
                   
                 Scanner logs 
               
               
                 K IC  - fracture toughness 
                   
                 HRA including high res 
               
               
                   
                   
                 sonic and lab toughness 
               
               
                 T 0  - tensile strength 
                   
                 HRA including high res 
               
               
                   
                   
                 sonic and lab tensile strength 
               
               
                 σ h  - min horizontal stress 
                   
                 Calibrated MEM (sonic, 
               
               
                   
                   
                 MDT stress) 
               
               
                 σ v  - overburden stress 
                   
                 Density logs 
               
               
                 p p  - pore pressure 
                   
                 Known from local field 
               
               
                   
                   
                 knowledge or measurements 
               
               
                 λ - coefficient of friction 
                 INTERFACE 
                 Lab measurements on cores 
               
               
                   
                   
                 or correlation to sonic 
               
               
                 K IIC  - fracture toughness 
                   
                 Lab measurements on cores 
               
               
                 (Mode II) 
                   
                 or correlation to sonic 
               
               
                 w int  κ i  - conductivity in 
                   
                 Lab measurements on cores 
               
               
                 intact zone 
               
               
                 w int  κ s  - conductivity in 
                   
                 Lab measurements on cores 
               
               
                 activated zone 
               
               
                   
               
             
          
         
       
     
         [0052]      FIG. 5  is a schematic drawing of a vertical hydraulic fracture (HF) growth in a subterranean layered rock. The HF propagates vertically (in the slide plane) and laterally (across the slide plane) by pumping of a fracturing fluid (in gray) from the well. Vertical propagation takes place upward and downward and characterized by the coordinates b 1  and b 2  respectively. The height growth in both sides is affected by the mechanical properties of the rock layers where the fracture tips are (e.g. fracture toughness), confining rock stresses, and hydromechanical properties of the interfaces between the adjoining layers (e.g. friction coefficient, fracture toughness, hydraulic conductivity). The HF propagation is associated with the leak-off of a fracturing fluid from the HF along the hydraulically conductive interfaces. 
         [0053]      FIG. 6  gives detailed overview of the families of input parameters and the names of every parameter in the family required for the HF simulator. 
         [0054]    Next, a discussion of a framework is needed. There are three main mechanisms related to the limitation of the HF growth in height: (i) the contrasts of the rock stresses and strengths between the adjoining rock layers (“mechanism 1” as introduced above), ( 201 ), (ii) the enhanced leak-off of the fracturing fluid into the bedding planes, presented here by the physical model ILeak ( 202 ) (sub-mechanism of “mechanism 3” as introduced above), and (iii) the elastic interaction with weakly cohesive slipping interfaces, presented here by the physical model FracT ( 203 ) (sub-mechanism of “mechanism 3” as introduced above). 
         [0055]      FIG. 7  presents an example of sequential HF growth in height affected by the interaction with weakly cohesive and conductive interfaces. The uniform HF growth is temporarily arrested by direct contact of the fracture tips with the upper and lower interfaces, meanwhile continuing its propagation laterally. After some delay of the HF tips at the interfaces, the HF reinitiate its vertical growth across them. The stages follow. 
         [0056]    Radial fracture: equal propagation in all directions 
         [0057]    Tips reach interface 
         [0058]    Vertical tips are temporarily arrested, horizontal tips continue to grow 
         [0059]    Fracture breaks interface and propagate vertically 
         [0060]      FIG. 8  demonstrates the HF height growth design workflow at a high level. It includes the input of the pre-given measured or estimated rock and interface properties on the one hand, and the input of the controlling parameters of the HF pumping schedule, on the other hand. They feed the model of the HF growth simulation ( 000 ), which is explained below. The results of the simulation go to the comparing module to find out the deviation of the simulated fracture height with respect to the optimum one. Depending on the tolerance of the fracture height growth obtained in the simulation, it either adjusts the fluid pumping parameters for the next cycle of the HF simulation, or outputs the used pumping parameters, which produce the optimum HF height in the given rock. 
         [0061]    Next, we discuss modeling of fracture propagation in a vertically heterogeneous layered medium. The implying fracture model has to provide a solution for the coupled system of equations for the mechanical response of the rock surrounding the fracture and viscous fluid flow injected into the fracture. It should be assumed that the finite strength of the rock and continuing fluid flow into the fracture will result in the propagation of the fracture tips (a contour in 3D geometry) and the injected fluid within the rock mass. Used equations describing the mechanics of both rock solid response and fluid flow within the fracture must be principally three-dimensional in order to account for the fracture growth in horizontal and vertical directions. Coupling of fracture propagation in both directions with the injected fluid volume will allow assessing fracture height containment in rock for the industrial volumes of injected fluid. 
         [0062]    Fracture model must take into account not only different stress and rock properties in different rock layers, but also interaction of the fracture tips with planes of weakness, such as bedding planes and lamination interfaces. It should be assumed that mechanical interaction of the hydraulic fracture with these interfaces can inevitably lead to creating zones of enhanced hydraulic permeability along these interfaces and significant fracturing fluid leak-off. Effect of weakness planes and enhanced interface permeability should be the key components of the intended computational model of fracture propagation in layered formations. 
         [0063]    Herein, we develop an extensive analytical model of hydraulic fracture interaction, crossing and subsequent growth across weak horizontal interfaces in the limiting case of low-viscous fluid friction (toughness-dominated regime). The latter is justified provided that the vertical fracture tip propagation velocity is reduced. We evaluate modified mechanical characteristics of a fracture such as net pressure, opening (width) and slippage zone extent when the fracture is deflected by an interface. Evaluation of the condition for crossing of the interface gives rise to finding out time delay of fracture termination at the interface. Overall picture of the intermittent character of fracture growth through a series of weakness planes is further used in the fluid-coupled description of fracture propagation in height in both plain-strain and three-dimensional elliptical fracture geometries. 
         [0064]    Construction of effective fracture propagation model in a finely laminated medium leads to the model of anisotropic medium with different fracture toughness in the vertical and horizontal directions. We estimate the aspect ratio of the length and height of the elliptical fracture in such medium for the given frictional and cohesional properties of the interfaces. The other mechanisms of fracture containment caused by the stress and rock property contrasts between layers can be applied on top of this model to use it in the modern fracture simulation tools. 
         [0065]      FIG. 9  explains the conceptual structure of the HF simulator ( 000 ). It consists of the input ( 100 ), explained in more detail above, simulation engine ( 200 ), and output ( 300 ). The simulation engine and output are explained in more detail below.  FIG. 10  depicts an embodiment of an algorithm of the HF simulator ( 200 ) workflow from the beginning of the fracturing job t 0  up to the end T. Every subsequent time step the fracture propagation problem is solved conventionally ( 201 ) such as there is no interaction with interfaces in the rock. Next, provided that the HF has contacted or crossed any rock interfaces the fracturing fluid leakoff module ILeak ( 202 ) is called to update the HF fluid volume, flowrate, and fluid pressure variations within the HF and infiltrated interfaces. Next, if the HF tip reaches any interface, the FracT module ( 203 ) is assessing the potential fracture tip arrest or crossing of the interface at the given time step. If the fracture tip is arrested, it remains non-propagating for the next time step. Otherwise, if the HF is crossing the interface or not contacted, it increments its length and goes to the next time step. 
         [0066]    The ILeak module ( 202 ) will be explained with more detail as follows. The input information includes the interface, contact pressure, fluid viscosity, and time step. The module operates at every change in time for all contacted or crossed interfaces. The module assumes no elastic interaction and that there is leakoff of fracturing fluid in the interfaces. The module computes the increment of fluid percolation with the given interface for a change in time and provides the fluid front, leaked off volume, and the flow rate into the interface. 
         [0067]    Consider an orthogonal junction of the vertical hydraulic fracture and a horizontal interface. The interface of finite thickness w int  is filled by a permeable material. The intrinsic permeability of the filling material in intact interface parts is κ i . Suppose that a certain segment of the interface, −b s &lt;x&lt;b s , nearby the junction is activated by shear displacement as a result of mechanical interaction with the hydraulic fracture. It results in the damage of the filling material within this segment and a change of its permeability to κ s  ( FIG. 11 ).  FIG. 11  shows the horizontal interface crossed by the vertical hydraulic fracture (top), and schematic distribution of the percolated fluid pressure along the interface (bottom). 
         [0068]    In tight formations κ i  can be negligibly small. This condition κ i =0) can be used later on to simplify the leak-off model. On the contrary, the activated part of the interface can be substantially more permeable than the intrinsic part due to the crushed grains of the filling material or shear dilation. Sliding activation of mineralized interfaces can be a dominant mechanism for the fracturing fluid leak-off in ultra-low permeability tight rocks. 
         [0069]    Let us assume that the fracturing fluid flow along a permeable interface is one-dimensional, steady and laminar. In these conditions it can be described by the following Darcy law 
         [0000]    
       
         
           
             
               
                 
                   
                     q 
                      
                     
                       ( 
                       x 
                       ) 
                     
                   
                   = 
                   
                     
                       - 
                       
                         w 
                         int 
                       
                     
                      
                     
                       κ 
                       μ 
                     
                      
                     
                       
                          
                         p 
                       
                       
                          
                         x 
                       
                     
                   
                 
               
               
                 
                   ( 
                   1 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where q(x) is the 2D rate of fluid percolation within the material of permeability κ, μ is the viscosity of the fluid, and p(x) is the fluid pressure distribution along the interface ( FIG. 11 , bottom). It is sometimes convenient to replace the product w int κ by the hydraulic conductivity of the interface c, typically measurable in laboratory (and use c s  and c i  notations hereafter, respectively). 
         [0070]    The total rate of the fracturing fluid leak-off from the hydraulic fracture into the particular interface at the junction point q L  is doubled due to symmetrical fluid diversion into both sides of the interface 
         [0000]        q   L =2(0)  (2)
 
         [0000]    Due to the symmetry of the fluid percolation in both sides of the interface, in what follows we obtain the solution only for the positive OX direction (x&gt;0). The Darcy law (1) establishes relationship between the local flow rate q and associated fluid pressure decay dp/dx at every point of a permeable material infiltrated by fluid. We write this law first for the flow rate q s  and pressure decay p s  within the activated (sheared) part as 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       
                         q 
                         s 
                       
                        
                       
                         ( 
                         x 
                         ) 
                       
                     
                     = 
                     
                       
                         - 
                         
                           
                             c 
                             s 
                           
                           μ 
                         
                       
                        
                       
                         
                            
                           
                             
                               p 
                               s 
                             
                              
                             
                               ( 
                               x 
                               ) 
                             
                           
                         
                         
                            
                           x 
                         
                       
                     
                   
                   , 
                   
                     x 
                     ≤ 
                     
                       min 
                        
                       
                         ( 
                         
                           
                             b 
                             f 
                           
                           , 
                           
                             b 
                             s 
                           
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
           
         
       
     
         [0000]    and for the fluid rate q i  and pressure p i  within the intact part of the interface 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       
                         q 
                         i 
                       
                        
                       
                         ( 
                         x 
                         ) 
                       
                     
                     = 
                     
                       
                         - 
                         
                           
                             c 
                             i 
                           
                           μ 
                         
                       
                        
                       
                         
                            
                           
                             
                               p 
                               i 
                             
                              
                             
                               ( 
                               x 
                               ) 
                             
                           
                         
                         
                            
                           x 
                         
                       
                     
                   
                   , 
                   
                     
                       b 
                       s 
                     
                     &lt; 
                     x 
                     &lt; 
                     
                       b 
                       f 
                     
                   
                 
               
               
                 
                   ( 
                   4 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where b f  is the front of percolated fluid. Outside of the zone of penetrated fluid we assume the in-situ pore pressure condition, i.e. 
         [0000]      ( x )=0, ( x )= p   p   , x≧b   f   (5)
 
         [0000]    The solution must include the position of the percolated fluid front b f  and the pressure profile (x) at every time of the leak-off process. 
         [0071]    From the fluid mass balance equation written for incompressible fluid within an interface with impermeable walls (except at the junction point) 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       
                         ∂ 
                         
                           ( 
                           
                             φ 
                              
                             
                                 
                             
                              
                             
                               w 
                               int 
                             
                           
                           ) 
                         
                       
                       
                         ∂ 
                         t 
                       
                     
                     + 
                     
                       
                         ∂ 
                         q 
                       
                       
                         ∂ 
                         x 
                       
                     
                   
                   = 
                   0 
                 
               
               
                 
                   ( 
                   6 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where φ is porosity of the filling material or natural interface asperities, q=q s (x) for x≦b s  and q=q i (x) for x&gt;b s , it follows that if the width w int  is constant (dw int /dt=0), the flow rate q has uniform value along the interface coordinate being only a function of time, i.e. 
         [0000]      ( x,t )=( x,t )= q ( x,t )=const( t )  (7)
 
         [0072]    Taking into account (7) and boundary condition (5) at x=b f , the solution of (3)-(4) for the distribution of the percolated fluid pressure (x) along the interface indicates a linear decay shown in  FIG. 12 .  FIG. 12  provides a profile of fluid pressure along the interface for the “in-slip” (top) and “out-of-slip” (bottom) regimes of percolation. 
         [0073]    The solution for the pressure profile is written separately for two regimes of fluid percolation into the interface: “in-slip” percolation, when the leaked fluid is totally contained within the slipped zone of the interface, i.e. b f ≦b s , and “out-of-slip” percolation into the intact interface zone also, i.e. b f &gt;b s . For the “in-slip” leak-off ( FIG. 12 , upper), we obtain the following linear pressure profile 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       p 
                        
                       
                         ( 
                         x 
                         ) 
                       
                     
                     = 
                     
                       
                         
                           p 
                           c 
                         
                         - 
                         
                           
                             ( 
                             
                               
                                 μ 
                                 
                                   κ 
                                   s 
                                 
                               
                                
                               
                                 
                                   b 
                                   . 
                                 
                                 f 
                               
                             
                             ) 
                           
                            
                           x 
                         
                       
                       = 
                       
                         
                           p 
                           c 
                         
                         - 
                         
                           
                             ( 
                             
                               
                                 
                                   p 
                                   c 
                                 
                                 - 
                                 
                                   p 
                                   p 
                                 
                               
                               
                                 b 
                                 f 
                               
                             
                             ) 
                           
                            
                           x 
                         
                       
                     
                   
                   , 
                   
                     x 
                     ≤ 
                     
                       b 
                       f 
                     
                     ≤ 
                     
                       b 
                       s 
                     
                   
                 
               
               
                 
                   ( 
                   8 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where p c =p(0) is the fluid pressure at “contact” with a hydraulic fracture, i.e. x=0. For the “out-of-slip” leak-off ( FIG. 12 , lower), we obtain the following broken line profile 
         [0000]    
       
         
           
             
               
                 
                   
                       
                   
                    
                   
                     
                       
                         p 
                          
                         
                           ( 
                           x 
                           ) 
                         
                       
                       = 
                       
                         
                           
                             p 
                             c 
                           
                           - 
                           
                             
                               ( 
                               
                                 
                                   μ 
                                   
                                     κ 
                                     s 
                                   
                                 
                                  
                                 
                                   
                                     b 
                                     . 
                                   
                                   f 
                                 
                               
                               ) 
                             
                              
                             x 
                           
                         
                         = 
                         
                           
                             p 
                             c 
                           
                           - 
                           
                             
                               ( 
                               
                                 
                                   
                                     p 
                                     c 
                                   
                                   - 
                                   
                                     p 
                                     1 
                                   
                                 
                                 
                                   b 
                                   s 
                                 
                               
                               ) 
                             
                              
                             x 
                           
                         
                       
                     
                     , 
                     
                       x 
                       ≤ 
                       
                         b 
                         s 
                       
                       ≤ 
                       
                         b 
                         f 
                       
                     
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       p 
                        
                       
                         ( 
                         x 
                         ) 
                       
                     
                     = 
                     
                       
                         
                           p 
                           1 
                         
                         - 
                         
                           
                             ( 
                             
                               
                                 μ 
                                 
                                   κ 
                                   i 
                                 
                               
                                
                               
                                 
                                   b 
                                   . 
                                 
                                 f 
                               
                             
                             ) 
                           
                            
                           
                             ( 
                             
                               x 
                               - 
                               
                                 b 
                                 s 
                               
                             
                             ) 
                           
                         
                       
                       = 
                       
                         
                           p 
                           1 
                         
                         - 
                         
                           
                             ( 
                             
                               
                                 
                                   p 
                                   1 
                                 
                                 - 
                                 
                                   p 
                                   p 
                                 
                               
                               
                                 
                                   b 
                                   f 
                                 
                                 - 
                                 
                                   b 
                                   s 
                                 
                               
                             
                             ) 
                           
                            
                           
                             ( 
                             
                               x 
                               - 
                               
                                 b 
                                 s 
                               
                             
                             ) 
                           
                         
                       
                     
                   
                   , 
                   
                     
                       b 
                       s 
                     
                     &lt; 
                     x 
                     &lt; 
                     
                       b 
                       f 
                     
                   
                 
               
               
                 
                   ( 
                   10 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where p 1 =p(b s ) is the fluid pressure at the slippage zone tip. In (8)-(10) we take into account that 
         [0000]        q=φw   int   {dot over (u)}=φw   int   {dot over (b)}   f   (11)
 
         [0000]    where {dot over (u)} is the lengthwise fluid velocity (upper dot stands for the differentiation with respect to time) equal to the velocity of the percolated fluid propagation {dot over (b)} f . Therefore, from (8)-(10) we obtain the following ordinary differential equations for the propagation of the fluid front (t) right after the contact (t&gt;t c ) for “in-slip” fluid penetration: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       
                         
                           b 
                           . 
                         
                         f 
                       
                        
                       
                         ( 
                         t 
                         ) 
                       
                     
                     = 
                     
                       
                         
                           κ 
                           s 
                         
                         μ 
                       
                        
                       
                         
                           
                             
                               p 
                               c 
                             
                              
                             
                               ( 
                               t 
                               ) 
                             
                           
                           - 
                           
                             p 
                             p 
                           
                         
                         
                           
                             b 
                             f 
                           
                            
                           
                             ( 
                             t 
                             ) 
                           
                         
                       
                     
                   
                   , 
                   
                     
                       b 
                       f 
                     
                     ≤ 
                     
                       b 
                       s 
                     
                   
                 
               
               
                 
                   ( 
                   12 
                   ) 
                 
               
             
           
         
       
     
         [0000]    and for “out-of-slip” penetration: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       
                         
                           b 
                           . 
                         
                         f 
                       
                        
                       
                         ( 
                         t 
                         ) 
                       
                     
                     = 
                     
                       
                         
                           
                             κ 
                             s 
                           
                           μ 
                         
                          
                         
                           
                             
                               
                                 p 
                                 c 
                               
                                
                               
                                 ( 
                                 t 
                                 ) 
                               
                             
                             - 
                             
                               
                                 p 
                                 1 
                               
                                
                               
                                 ( 
                                 t 
                                 ) 
                               
                             
                           
                           
                             b 
                             s 
                           
                         
                       
                       = 
                       
                         
                           
                             κ 
                             i 
                           
                           μ 
                         
                          
                         
                           
                             
                               
                                 p 
                                 1 
                               
                                
                               
                                 ( 
                                 t 
                                 ) 
                               
                             
                             - 
                             
                               p 
                               p 
                             
                           
                           
                             
                               
                                 b 
                                 f 
                               
                                
                               
                                 ( 
                                 t 
                                 ) 
                               
                             
                             - 
                             
                               b 
                               s 
                             
                           
                         
                       
                     
                   
                   , 
                   
                     
                       b 
                       f 
                     
                     &gt; 
                     
                       b 
                       s 
                     
                   
                 
               
               
                 
                   ( 
                   13 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where the fluid pressure at the slip zone tip p 1 =p(b s ) is found as 
         [0000]    
       
         
           
             
               
                 
                   
                     p 
                     1 
                   
                   = 
                   
                     
                       p 
                       p 
                     
                     + 
                     
                       
                         ( 
                         
                           
                             p 
                             c 
                           
                           - 
                           
                             p 
                             p 
                           
                         
                         ) 
                       
                        
                       
                         
                           
                             b 
                             f 
                           
                           - 
                           
                             b 
                             s 
                           
                         
                         
                           
                             b 
                             f 
                           
                           - 
                           
                             
                               ( 
                               
                                 1 
                                 - 
                                 
                                   κ 
                                   is 
                                 
                               
                               ) 
                             
                              
                             
                               b 
                               s 
                             
                           
                         
                       
                        
                       
                         h 
                          
                         
                           ( 
                           
                             
                               b 
                               f 
                             
                             - 
                             
                               b 
                               s 
                             
                           
                           ) 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   14 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where κ is =κ i /κ s , and H(x) is the Heaviside step function (zero for negative, and one for positive arguments respectively). 
         [0074]    The solution of (12)-(13) is found for both regimes of fluid penetration as follows 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       
                         b 
                         f 
                       
                        
                       
                         ( 
                         t 
                         ) 
                       
                     
                     = 
                     
                       
                         
                           b 
                           
                             f 
                              
                             
                                 
                             
                              
                             1 
                           
                         
                          
                         
                           ( 
                           t 
                           ) 
                         
                       
                       = 
                       
                         
                           
                             
                               2 
                                
                               
                                   
                               
                                
                               
                                 κ 
                                 s 
                               
                             
                             μ 
                           
                            
                           
                             
                               ∫ 
                               
                                 t 
                                 c 
                               
                               t 
                             
                              
                             
                               Δ 
                                
                               
                                   
                               
                                
                               
                                 
                                   p 
                                   c 
                                 
                                  
                                 
                                   ( 
                                   
                                     t 
                                     ′ 
                                   
                                   ) 
                                 
                               
                                
                               
                                   
                               
                                
                               
                                  
                                 
                                   t 
                                   ′ 
                                 
                               
                             
                           
                         
                       
                     
                   
                   , 
                   
                     
                       b 
                       f 
                     
                     ≤ 
                     
                       b 
                       s 
                     
                   
                 
               
               
                 
                   ( 
                   15 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       
                         b 
                         f 
                       
                        
                       
                         ( 
                         t 
                         ) 
                       
                     
                     = 
                     
                       
                         
                           
                             κ 
                             is 
                           
                            
                           
                             ( 
                             
                               
                                 
                                   b 
                                   
                                     f 
                                      
                                     
                                         
                                     
                                      
                                     1 
                                   
                                   2 
                                 
                                  
                                 
                                   ( 
                                   t 
                                   ) 
                                 
                               
                               - 
                               
                                 
                                   ( 
                                   
                                     1 
                                     - 
                                     
                                       κ 
                                       is 
                                     
                                   
                                   ) 
                                 
                                  
                                 
                                   b 
                                   s 
                                   2 
                                 
                               
                             
                             ) 
                           
                         
                       
                       + 
                       
                         
                           ( 
                           
                             1 
                             - 
                             
                               κ 
                               is 
                             
                           
                           ) 
                         
                          
                         
                           b 
                           s 
                         
                       
                     
                   
                   , 
                   
                     
                       b 
                       f 
                     
                     &gt; 
                     
                       b 
                       s 
                     
                   
                 
               
               
                 
                   ( 
                   16 
                   ) 
                 
               
             
           
         
       
     
         [0075]    where t c  is the time at the beginning of the fracture-interface contact, Δp c (t′)=p c (t′)−p p  is the differential fluid pressure at the interface. The evolution of the differential pressure with time therefore dictates the leak-off process in the given contacted interface. 
         [0076]    Consider a vertical plane-strain fracture pumped by a constant injection rate and growing symmetrically upward and downward in a homogeneous rock. Let a permeable interface be placed at some distance y=h c  from the injection point y=0. Once the height of the fracture reaches h=h c , the fluid begins to percolate into the interface. At time t=t c , the fracture may stop or continue growing with given leak-off as shown in  FIG. 13 .  FIG. 13  shows hydraulic fracture propagating upward and downward in plane-strain geometry (vertical cross-section). There are three distinct stages: (left) pre-contact with growing fracture without leak-off, (middle) early contact with non-growing fracture with leak-off, and (right) late contact with growing fracture with leak-off. 
         [0077]    We will suppose that prior to a direct contact with an interface at t=t c  the hydraulic fracture propagates without any elastic or hydraulic interaction. The remotely placed permeable interface is not mechanically activated due to the approaching fracture and thus it does not change the stress state around. Before the contact, the injected fluid is totally contained within the fracture, as the medium is supposed impermeable. Right after the contact with the interface (t=t c ), the fluid flows within the interface and causes a loss of fluid volume stored in the hydraulic fracture. The fracture continues grow once the fluid volume loss is compensated by the injected volume at a later time t=t r &gt;t c . We provide a detailed example of the mechanics of the fracture propagation affected by the presence of a hydraulically conductive interface on the path of its height growth on  FIG. 14 . 
         [0078]      FIG. 14 . The injected, fracture and leaked-off fluid volumes (top), net pressure (middle), and hydraulic fracture halfheight (bottom) during the whole cycle of fluid injection into the fracture. The left time region shaded in blue is the pre-contact stage. The middle time region shaded in orange is the early contact stage. The right time region shaded in green is the late contact stage. At the very beginning (in blue shaded time stage) the hydraulic fracture propagates without interaction and leak-off. The net pressure decline and fracture height growth follow the expected behavior. Right after the contact with a permeable plane (yellow shaded time stage), the leak-off starts following the known asymptotic behavior. Initially it dominates over the injection as predicted from leak-off equation above, and the fracture fluid volume ν partly drops. The rate of leak-off into the interface gradually reduces with time of percolation. During the early contact stage the leakoff rate becomes smaller than the injection rate into the fracture. This restores the fluid volume increase within the hydraulic fracture that it had lost at the moment of contact. When the fluid volume losses due to the leak-off are totally compensated by the postcontact injection into the fracture, the critical net pressure is achieved within the fracture again and it reinitiates its vertical growth (green shaded time region). At a late contact stage, the fracture growth takes place with continuing leak-off. The rate of fracture volume pumping is therefore less than it was before contact, so the decay of net pressure and velocity of fracture height growth are also smaller. If the leak-off takes place only into one interface, the rates of fracture growth will return back to initial values with time when the leak-off becomes negligibly small, and can be fully neglected in simulations. 
         [0079]    Next, we discuss the methods, inputs, and outputs of the FracT module ( 203 ). The inputs include the upper or lower tip coordinates, pressure profile, formation layers and interfaces, and the index of the interface at a T-shaped contact. The module provides a slip boundary, residual slip, and interface state of intact, T-shaped, or crossed. The FracT module is called for every interface at a T-shaped contact with the fracture tip and includes elastic interaction and crossing criterion and re-initiation past-interface. 
         [0080]    Consider the vertical cross-section of a hydraulic fracture growing in height ( FIG. 15 , left). Suppose that both fracture tips simultaneously reach two pre-existing horizontal interfaces above and below. After the contact, the interfaces slip and arrest further fracture tip propagation in the vertical direction ( FIG. 15 ).  FIG. 15  provides a two-sided contact of a vertically growing fracture and weak horizontal interfaces (left), interface activation, and fracture tip blunting as a result of the contact with the interfaces (right). 
         [0081]    At the point of contact, the problem becomes the one of the orthogonal contact between a pressurized fracture and two weak interfaces, shown in  FIG. 15  (right). To solve this problem, we first need to obtain the modified fracture characteristics such as the fracture volume, the opening (width), the blunting characteristics of the tip, the extent of the interfacial slip zone b s , and the associated drop of the net pressure within the fracture after the contact. Next, we need to evaluate the minimum buildup of net pressure necessary to cross the interfaces. This criterion of interface crossing can then be used, for example, in rigorous 3D fracture propagation models, where it will quantify the time delay of the fracture height growth due to interfacial contacts (i.e. from the moment of fracture contact with the interface and its arrest to the subsequent crossing of the interface to continue the propagation). 
         [0082]    The problem of an elasto-frictional fracture contact can be solved rigorously numerically. Here, we use an approximate analytical solution of this problem, described in more detail in SPE-173337, “Hydraulic Fracture Height Containment by Weak Horizontal Interfaces,” February 2015, by Dimitry Chuprakov and Romain Prioul, which is incorporated by reference herein. The analytical model facilitates parametrical insights into the fracture contact problem. We focus on the following characteristics of the fracture-interface contact: (i) the extent of the interface activation in shear b s , (ii) the associated hydraulic fracture opening w T  (width) at the junction with the interface, and (iii) the post-contact fracture volume V in the vertical cross section. These characteristics are found to be functions of the fracture net pressure p′, the critical shear stress at the slipping part of the horizontal interface {tilde over (τ)} m , the interface fracture toughness κ IIC   (INT) , and the half-height of the pressurized vertical fracture L. To facilitate the formulation of the problem in dimensionless form, we introduce the relative length of the interface activation β s =b s /L, the modified fracture opening at the contact Ω T =w T  E′/4, and the modified fracture volume ν=VE′/(2π), where E′=E/(1−ν 2 ) is the modified planestrain Young modulus, and they can be expressed as 
         [0000]      β s =(Π,κ IIC ), Ω T =Ω m   Ω (Π,κ IIC ), ν=ν 0   ν (Π,κ IIC )  (17)
 
         [0000]    where ν 0 =p′L 2  is the modified fracture volume, and Ω m =p′L is the maximum modified fracture opening at the middle of the fracture prior to contact. The two dimensionless parameters are the relative net pressure Π=p′/τ m  and the dimensionless interface toughness κ IIC =κ IIC   (INT) /(τ m √{square root over (πL)}), where τ m =λσ′ v , λ is the coefficient of friction, and σ′ v =σ v −p int  is the effective vertical stress at the interface with interstitial fluid pressure p int . Initially, p int  equals the pore pressure; after fracturing fluid penetration into the interface, it represents the pressure of the penetrated fluid. 
         [0083]    The magnitude of the relative net pressure Π defines the magnitude of these characteristics. The size of the interface activation monotonically increases with Π. It is small when the net pressure p′ is small or frictional stress τ m  is large. In most practical cases, when the net pressure is small relative to the frictional stress (Π=p′/τ m &lt;&lt;1), the activated zone obeys the following asymptote 
         [0000]    
       
         
           
             
               
                 
                   
                     β 
                     s 
                   
                   ≅ 
                   
                     
                       
                         π 
                         2 
                       
                       8 
                     
                      
                     
                       
                         ( 
                         
                           Π 
                           - 
                           
                             κ 
                             IIC 
                           
                         
                         ) 
                       
                       2 
                     
                   
                 
               
               
                 
                   ( 
                   18 
                   ) 
                 
               
             
           
         
       
     
         [0000]    In the opposite limit of relatively high net pressures (Π&gt;&gt;1), we arrive at the following linear asymptote 
         [0000]    
       
         
           
             
               
                 
                   
                     β 
                     s 
                   
                   ≅ 
                   
                     
                       2 
                       π 
                     
                      
                     Π 
                   
                 
               
               
                 
                   ( 
                   19 
                   ) 
                 
               
             
           
         
       
     
         [0000]    A similar trend is observed for the fracture opening (width)  Ω   T =Ω T /Ω m  at the junction. The fracture tends to close at the contact with the interface, if Π−κ IIC &lt;&lt;1, following the asymptote 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       Ω 
                       _ 
                     
                     T 
                   
                   ≅ 
                   
                     
                       π 
                       4 
                     
                      
                     
                       
                         
                           Π 
                           2 
                         
                         - 
                         
                           κ 
                           IIC 
                           2 
                         
                       
                       Π 
                     
                   
                 
               
               
                 
                   ( 
                   20 
                   ) 
                 
               
             
           
         
       
     
         [0000]    In the opposite limit (Π&gt;&gt;1), the opening at the junction is of the same order of magnitude as the maximum opening, Ω m . It changes logarithmically with Π, as follows 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       Ω 
                       _ 
                     
                     T 
                   
                   ≅ 
                   
                     
                       2 
                       π 
                     
                      
                     
                       ln 
                        
                       
                         ( 
                         
                           
                             2 
                             π 
                           
                            
                           Π 
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   21 
                   ) 
                 
               
             
           
         
       
     
         [0084]    In the case of simultaneous fracture contact with two weak interfaces, the profile of the fracture opening widens as a function of Π as shown in  FIG. 16  (left).  FIG. 16  provides profiles of the vertical fracture opening at the contact with two cohesionless interfaces (grey) for the relative net pressure Π equal to 0.1 (black), 1 (blue), and 10 (red) (left), and the relative net pressure Π in the fracture prior to (dashed line) and after (solid lines) the contact with interfaces versus the normalized fracture volume ν/(τ m L 2 ) for the case of two-sided fracture contact (right). Black lines represent the normalized fracture toughness along interfaces κ IIC =0, and red lines are for κ IIC =0.1. Blue arrows denote associated pressure drop within the fracture at the moment of contact with interfaces. 
         [0085]    The larger the relative net pressure Π, the wider fracture opens along the entire vertical cross section, as expected. The effect of interfaces on the elastic fracture opening resembles a sudden change of the elastic compliance of the rock. Indeed, the weakness planes represent two compliant planes in a stiff rock. When the fracture establishes contact with them, it is obvious that the elastic response of the fracture must become more compliant. This effect of abrupt fracture widening at the moment of contact with weak interfaces may result in an abrupt drop of the fracture pressure. Fast increase of the fracture volume must lead to an associated fast decrease in the fluid pressure. We performed additional investigations of the net pressure drop at the moment of the fracture contact with two weak interfaces.  FIG. 16  (right) shows the magnitude of the relative net pressure drop for the given volume of injected fluid within the fracture immediately prior to the contact with the interfaces. When the relative net pressure is small (Π&lt;1), the pressure drop is small and not detectable. For large relative net pressures (Π&gt;1), the pressure within the fracture drops noticeably. Herein, the fracture opening profile is found as a part of the problem solution. 
       Fracture Reinitiation Problem: Crossing of the Interfaces. 
       [0086]    The interface activation generates a localized tensile stress field on the opposite side of the interface ( FIG. 17 ). High tensile stresses are concentrated close to the junction point and can exceed tensile strength of the formation. In the most stress-perturbed region, the maximum principal tensile stress component is parallel to the interface. The contact-induced stresses favor the initiation of a new tensile crack in the intact rock in a direction normal to the interface (see arrows in  FIG. 17 ). A similar problem has been solved analytically assuming uniform opening of the fracture.  FIG. 17  includes the maximum tensile stress component generated on the opposite side of the cohesionless (left) and cohesional interface with κ IIC =1 (right). Vertical and horizontal white solid lines depict the fracture and interface, respectively. White arrows point out local directions of the maximum principal compressive stress (perpendicular to the maximum principal tensile stress). The coordinate scales are all normalized on the extent of sliding zone b s    
         [0087]    To initiate a new crack and cross the interface, sufficient elastic strain energy must also be accumulated in the rock. Critical stress and critical elastic energy release are both required for crack initiation in solids. To use this mixed stress-and-energy criterion for the fracture reinitiation, we derive and evaluate the initiation stress intensity factor K ini  within the critical stress zone as a function of the problem parameters. Then, we introduce the following crossing function, Cr, as the ratio of the initiation stress intensity factor K ini  and the fracture toughness of the rock behind the interface K IC   (2) , where the crack is to be initiated: 
         [0000]    
       
         
           
             
               
                 
                   Cr 
                   = 
                   
                     
                       
                         K 
                         ini 
                       
                       
                         K 
                         IC 
                         
                           ( 
                           2 
                           ) 
                         
                       
                     
                     = 
                     
                       
                         
                           K 
                           IC 
                           
                             ( 
                             1 
                             ) 
                           
                         
                         
                           K 
                           IC 
                           
                             ( 
                             2 
                             ) 
                           
                         
                       
                        
                       
                         
                           Cr 
                           ^ 
                         
                          
                         
                           ( 
                           
                             Π 
                             , 
                             
                               κ 
                               IIC 
                             
                             , 
                             α 
                           
                           ) 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   22 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where α=σ h /τ m  is the relative minimum horizontal stress σ h  in the layer behind the interface. The crossing function Cr is greater than 1 if the crossing criterion is satisfied, otherwise the fracture is arrested at the interface. The contrast of fracture toughnesses on both sides of the interface, K IC   (1) /K IC   (2)  plays an important role as expected. The fracture growth into a weaker formation is less resistant as opposed to the growth into a stronger rock. We further consider a particular case of identical rock toughnesses on both sides of the interface (K IC   (1) =K IC   (2) ). To understand the possible delay of the fracture tip growth at the interface, we investigate the dependence of the modified crossing function Cr=C{circumflex over (r)} on the dimensionless parameters of the problem: Π, κ IIC  and α. 
         [0088]    Consider the initial moment of the contact with the interface. It appears that for all values of the dimensionless parameters of the problem, the crossing function is initially less than 1. This means that interface can never be crossed straight away as a continuous fracture propagation process. The fracture tip is arrested by the interface until the net pressure builds up sufficiently to raise the value of crossing function to 1. One can understand this from a mechanical fracture energy perspective. The noninteractive fracture tip requires additional injected fluid energy to grow. Once the contact with the interface is established, part of the fracture energy is consumed into the energy required for the interface slippage. Therefore, the crossing of the interface requires more energy than is required in the noninteraction case. This explains the abrupt stop of the fracture tip at a weak interface. 
         [0089]    The above results on the interface crossing pertain to the two-sided hydraulic fracture contact problem. In the considered examples, the fracture half height L is therefore assumed fixed after the contact. In the general case, the fracture can interact with only one interface while the other vertical fracture tip continues to grow. This general case has been solved using a similar technique and shows that containment at the interfaces will follow the same trends in net pressure behavior. 
         [0090]    Intermittent Fracture Propagation Through Interfaces (LamiFrac Model) 
         [0091]    Next, we explore the impact of the previous mechanism on the 3D planar hydraulic fracture propagation from a horizontal well in a multilayer formation with horizontal weak interfaces on both sides of the well (we consider a symmetric case for simplicity, although the methodology is general). Within each layer, the stresses, rock elastic and strength properties do not change but they are allowed to vary between layers. The fracture propagation starts from a small circular fracture. Please refer back to  FIG. 1  which illustrates the geometry of the layers and interfaces and the hydraulic fracture. 
         [0092]    Initially the hydraulic fracture propagates equally in the upper vertical, lower vertical, and horizontal directions (i.e., as a radial fracture at the start). Then, following the contact with the interfaces, the propagation in the horizontal and vertical directions becomes different. For the sake of the demonstration, here we use an approximate solution of the 3D fracture problem based on the solution for an elliptical crack. The fracture geometry keeps an elliptical shape given the unequal growth in the three directions (two vertical and horizontal one). The modeling algorithm consists of three computational components. The first one computes the elastic fracture response to the injected fluid pressure and in-situ stress. It accounts for the fracture interaction with the interfaces as presented above. The second component solves for the simultaneous fracture tip growth in all three directions. The third component finds the fluid pressure within the fracture and all contacted interfaces, given the conditions for the fluid injection rate, the leakoff along the conductive interfaces, and the viscous fluid friction within the fracture. The latter obeys the known lubrication law for Newtonian fluids. 
         [0093]    In the simulations, we first prescribe the parameters of the rock and fluid injection in the borehole. Then, we compute the evolution of the fracture propagation geometry for the prescribed conditions, which enables us to investigate the impact of the pre-existing horizontal interfaces on the fracture containment. 
         [0094]    The qualitative picture of the fracture propagation looks similar in all simulations and can be described as follows. Once the vertical tips reach the upper and lower interfaces, their propagation stops for some time. The fracture still continues to propagate in the horizontal direction. At this stage, the net pressure in the fracture builds up (in similar fashion as one would observe in a PKN-type fracture). Once the net pressure has increased up to a critical value, the fracture has enough energy to break the interfaces. After the crossing of the interfaces, the fracture immediately contacts the next interfaces. As the fracture jumps vertically from one interface to another, the net pressure drops. As a result, the fracture growth temporarily ceases in all directions. Under further pressure increase, the fracture continues to grow in the horizontal direction again while it is still arrested in the vertical direction, and this growth leads to additional pressure buildup. The crossing of the interfaces and next cycle of pressure drop repeats itself. Such intermittent fracture propagation continues as long as the fracture interacts with horizontal interfaces. 
         [0095]      FIG. 18  illustrates the described mechanics of the fracture tip propagation and the pressure oscillations. It shows the results of two simulations with small and large injected fluid viscosity (1 cP and 10000 cP, respectively). The spacing between the interfaces is 0.1 m. For simplicity, the rock and interface properties within each layer are identical in these runs. These simulations show ( FIG. 18 , top) that the hydraulic fracture&#39;s vertical growth is inhibited due to the presence of weak interfaces. 
         [0096]    As a result, the fracture grows preferentially in the horizontal direction. The increased viscosity in the fluid injected into the fracture favors the interface crossing, which is well known. This explains why the containment effect is less prominent with larger fluid viscosity ( FIG. 18 , top right).  FIG. 18  shows fracture tip propagation (top) and inlet pressure decline (bottom) in the case of an elliptical fracture with Newtonian fluid with viscosity of 1 cP (left) and 10000 cP (right), respectively. Constant rate of the fluid injection into the fracture is 0.001 m 2 /s. The radius of initial fracture is 1 cm. Spatial spacing between horizontal interfaces is 0.1 m. The interfaces are cohesionless with the coefficient of friction 0.6 and pore pressure 12 MPa. Vertical in-situ stress is 20 MPa, minimum horizontal in-situ stress is 15 MPa. The fracture toughness of the rock is K IC =1 MPa*m 1/2 , tensile strength 5 MPa, E′=10 GPa. 
         [0097]    In the limiting case of a finely laminated structure, the pressure oscillations and tip jumps become vanishingly small. The fracture growth then represents a continuous process. The description of the fracture propagation in these rocks can be similar to that in a homogeneous rock, with the only difference being that the fracture toughness in the vertical direction across the interfaces has an increased “effective” value. The envelopes of the pressure curves for an “effective” finely laminated structure with weak interfaces and a continuous homogeneous rock without interfaces are plotted in  FIG. 18  (red and green curves, respectively). These pressure curves clarify the difference between the effect of the fracture toughness across a multilaminated/multilayered formation and the one without interfaces. 
         [0098]    Using the model above, we obtain the “effective” fracture toughness for laminated formations. The steady fracture propagation criterion requires that the stress intensity factor K I  at the tip equals the fracture toughness of the rock K IC : 
         [0000]        K   I   =K   IC   (23)
 
         [0000]    In a laminated formation, the steady growth in height means that the vertical tip constantly crosses the infinitesimally close interfaces, so that Cr=1 (Eq. 22). Rewriting this equation in terms of the stress intensity factor at the vertical tip, we have 
         [0000]        K   I   =K   IC   (eff)   (24)
 
         [0000]    where K IC   (eff) =K IC K I /K ini  is the “effective” fracture toughness. It is always larger than K IC  and depends on the mechanical properties of the interfaces, such as cohesion, friction coefficient, and hydraulic conductivity. This result is in agreement with the laboratory measurements of in- and cross-layer toughness used in the previous models. 
         [0099]      FIG. 19  builds a workflow for the conventional HF propagation solver ( 201 ) such as if there is no interaction with rock interfaces (but it includes the stress and strength constrast mechanism 1). The coupled solid-fluid HF solver ( 211 ) is called for every guessed increment of the fracture tip to output the solution for the stress intensity factor (SIF) K I  at the HF tip. The SIF is then compared with the fracture toughness of the present rock layer K IC  to find out if the fracture tip is stable or not. The loop is reinitiated every time unless the current increment of the HF tip is stable, and outputs the found solution. 
         [0100]      FIG. 20  builds a workflow for the sub-component ( 211 ) of the HF propagation solver ( 201 ) above. It represents a coupled solid-fluid HF solver for the given placement of the HF tips. It takes the solution for the HF at the previous time step ( 2111 ), finds out the coupled solution of elasticity ( 2112 ) and fluid flow ( 2113 ) at the next new time step and new fracture tips, and outputs it ( 2114 ). The coupled solution for the elasticity ( 2112 ) and fluid flow ( 2113 ) requires additional iterations (horizontal arrow between  2112  and  2113 ). 
         [0101]      FIG. 21  shows the output sub-modules of the main workflow ( 300  at  FIG. 9 ). They are geometrical ( 301 ), e.g. HF height and length, informational about the affected rock interfaces ( 302 ), e.g. coordinates of the crossed interfaces and generated slips at each of them, and mechanical ( 303 ), e.g. fluid pressure and fracture aperture.