DETERMINING HYDRAULIC FRACTURE PROFILES IN MULTI-LAYERED SUBSURFACE FORMATIONS

Description

TECHNICAL FIELD

The present disclosure describes systems and methods for hydrocarbon extraction. The systems and processes are configured for determining hydraulic fracture profiles in multi-layered formations accounting for non-constant horizontal in-situ stresses and hydrostatic fluid pressures in each layer to accurately represent subsurface conditions in multi-layered subsurface formations during fracturing operations.

BACKGROUND

Tight and shale reservoirs have become important because of their vast distribution and hydrocarbon accumulation. Hydraulic fracturing treatment can be used to stimulate production from tight reservoirs with low permeability, like shale, tight sandstone, etc. An issue related to hydraulic fracturing is fracture height and fracture width. Fracture height can be used to optimize the landing depth of horizontal wells in multi-layered formations. Accurate control of the hydraulic fracture height growth is important not only for maximizing the production performance but also for preventing the fractures from growing into non-production zones. If uncontrolled, hydraulic fractures can propagate into overlying and underlying layers and lead to negative consequences on the effectiveness of hydraulic fracturing operations and an uneconomical well.

SUMMARY

The systems and methods disclosed herein are configured to control systems for hydrocarbon extraction based on predictions of hydraulic fracture profiles (for example, fracture height and fracture width) in multi-layered formations. The predictions account for depth-dependent horizontal in-situ stresses exerted on each layer of the formation and depth-dependent hydrostatic fluid pressures exerted on the surfaces of the fracture profile. A data processing system (for example, a computer) achieves the accuracy of the predictions by discretizing the multi-layer formation in the vertical direction using one-dimensional high-order elements (for example, having at least second order interpolation). The data processing system uses the one-dimensional elements to calculate stress intensity factors at two fracture tips (for example, lower and upper tips or portions).

Based on the predicted fracture profiles, the data processing system can optimize a hydraulic fracturing pump schedule, control fracturing height within a production zone, predict fracture width, and select an optimal proppant size and type for admittance into the hydraulic fracture surface. For example, based on the predicted fracture profiles, the data processing system can determine optimal locations for perforations and where to extract hydrocarbons from the formation. Improperly selected perforation locations coupled with improper fluid pressures can lead to fractures reaching overburden or leaving the intended production zone which can have negative consequences.

The systems and methods described herein can accommodate subsurface scenarios including non-constant horizontal in-situ stress distributions and non-constant hydrostatic pressure distributions to determine fracture profile. In some examples, predictions can have limited value when constant horizontal in-situ stress and a linear distribution hydrostatic pressure within the fracture profile is assumed, which may not and usually cannot reflect the real subsurface conditions. The systems and methods disclosed herein increase the chances of success of hydraulic fracturing design for well stimulation by accounting for depth-dependent horizontal in-situ stresses and depth-dependent hydrostatic fluid pressures.

Hydraulic fracturing treatment stimulates shale and tight reservoirs when these reservoirs cannot produce hydrocarbons naturally after drilling. The data processing system described herein is configured to simulate and cause hydraulic fractures to form within a pay zone such that the hydraulic fractures propagate within the pay zone as long as possible. Accurate control of the hydraulic fracture height growth is necessary not only for maximizing the production performance but also for preventing the hydraulic fractures from growing into non-production zones. The data processing system is configured to predict hydraulic fracture height growth to facilitate hydraulic fracturing design. Many factors have been identified to impact the hydraulic fracture height growth, such as the minimum horizontal stress contrast between different formation layers, fracture toughness, bedding planes, mechanical properties, proppant distribution. Among these factors, the minimum horizontal stress contrast between the layers is often considered the dominant factor in controlling the hydraulic fracture height.

Numerical modeling of the hydraulic fracture height and width growing in multi-layered formations has been a challenging task in the oil and gas industry due to the complex subsurface conditions. The hydraulic fracture height growth and lateral propagation in a laminated reservoir is generally a fully three-dimensional (3D) problem that involves rock deformation of the fracture coupled with fluid flow inside the fracture. Fully coupled rock deformation and fluid flow for fracture propagation requires very fine grids (for example, finite element meshes), which can be computationally very expensive to solve. Therefore, the problem is often simplified and reduced to the so-called pseudo three-dimensional (P3D) problem, which includes two parts: hydraulic fracture height growth and the lateral propagation. The main assumptions of the models are that each cross-section perpendicular to the lateral propagation of the hydraulic fracture is in a plane strain condition and that the clastic properties throughout the layers are uniform. To estimate the fracture height in multi-layered formations, the solutions to the equilibrium height models have been investigated. Typically, these solutions start from the simplest symmetric three-layer problem and evolve to a complex asymmetric multiple-layer problem considering in-situ stress and other parameters such as fracture toughness and fluid density. The systems and methods disclosed herein do not need to assume linear hydrostatic pressure distribution and constant minimum horizontal stress within each layer.

The systems and methods disclosed herein can be used to predict the fracture height growth and fracture width versus increasing fracture mouth pressure for hydraulic fracturing treatment in multi-layered formations. The hydrostatic pressure and minimum horizontal stress within each layer need not be limited to be constant or linear distribution. The systems and methods dynamically discretize the multilayer formations along the vertical direction using onc-dimensional high-order elements, which can result in highly accurate stress intensity factors at the two fracture tips (for example, lower and upper tips). The number of one-dimensional elements is dependent on the minimum horizontal stress distribution along the vertical depth and the type of clements being used (for example, linear, quadratic, etc.). The fracture height model can be used to determine the landing depths of horizontal wells in tight and shale gas developments.

The systems and methods calculate hydraulic fracture height and fracture width for fracturing a multi-layer formation. The methods are applicable to nonlinear distributions of hydrostatic pressure and minimum horizontal stress along a vertical direction. The accuracy can be achieved using high order one-dimensional clements with smaller clement size to discretize the multi-layered formations along the vertical direction. A fracture height map of a given hydraulic fracture growing with increasing perforation pressure in the multi-layered formation can be displayed, which offers important information for hydraulic fracturing design. The systems and methods can be used to optimize the hydraulic fracturing pump schedule design to control the downhole pressure so that hydraulic fracture height contained within the production zone and the fracture width. The systems and methods can be used to predict the fracture width for selecting the right proppant size and type for admittance into the fracture surface, which determines how much of the hydraulic fracture width created by the fracturing treatment will be retained as propped width once the hydraulic fractures closes. The systems and methods can be used to determine where the horizontal well should be landed in the reservoirs for good production.

The details of one or more implementations of the subject matter described in this disclosure are set forth in the accompanying drawings and the description below. Other features, aspects, and advantages of the subject matter will become apparent from the description, the drawings, and the claims.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is diagram of an example wellbore system.

FIG. 2 is a plot of a hydraulic fracture propagating in a multi-layered formation.

FIG. 3 is a plot of perforation pressure and hydrostatic pressure distribution along the fracture height.

FIG. 4 is a plot of a discretized hydraulic fracture from the lower tip to the higher tip.

FIG. 5 is a plot of a one dimensional second order element.

FIG. 6 is a flowchart of a method to predict a hydraulic fracture profile for a perforation location.

FIGS. 7A and 7B are plots of hydraulic fracture height map versus perforation pressure for constant in-situ stresses within each layer and non-constant in-situ stresses within each layer, respectively.

FIGS. 8A-8C are plots of hydraulic fracture widths for perforation pressures of 41.4 MPa (6000 psi), 44.1 MPa (6400 psi), and 46.5 MPa (6750 psi), respectively.

FIG. 9A is an illustration of a hydraulic fracture without proppant.

FIG. 9B is an illustration of a hydraulic fracture with proppant.

FIG. 10 is a flowchart of an example method for fracturing a formation based on a predicted fracture profile.

FIG. 11 is a schematic of an example controller for determining perforation locations and a suggested fracturing method for each perforation location.

DETAILED DESCRIPTION

The systems and methods disclosed herein are configured to predict hydraulic fracture profiles in multi-layered formations. The hydraulic fracture profiles account for depth-dependent horizontal in-situ stresses exerted on each layer of the formation and depth-dependent hydrostatic fluid pressures exerted on the surfaces of the hydraulic fracture profile. In some examples, the data processing system predicts fracture profiles by discretizing the multi-layer formation in the vertical direction using one dimensional high-order elements to calculate stress intensity factors at two fracture tips.

Wellbore System

FIG. 1 is a diagram of a wellbore system 100. The wellbore system 100 includes a well 101 with a wellbore 102 that has been drilled into one or more layers (for example, layers 104A-104E) of a multi-layered formation 105 (e.g., by a drilling rig). The wellbore 102 includes a substantially vertical section 102A and a substantially horizontal section 102B. In some examples, the wellbore 102 is drilled such that the horizontal section 102B lies in, or near, a reservoir. Hydrocarbons (e.g., oil and gas) are extracted from the reservoir to a ground surface of the well 101.

Wellbore system 100 includes a perforation location 106 at a depth D1 within the wellbore 102. D1 can represent the true vertical depth (TVD) in the wellbore 102 or it can represent the measured depth in the wellbore 102 (e.g., as measured by a downhole logging device 103). In the example shown, depth represents true vertical depth in the wellbore 102. While wellbore 102 is illustrated to include one perforation location 106, some wells have more than one (e.g., 2-10) perforation locations. In some examples, one or more perforations in the wellbore wall are formed by charged explosives at the perforation locations.

Determining the perforation location in the wellbore 102 is an important aspect of well placement and hydraulic fracturing design. In some examples, it is preferable to include as many perforation locations as possible to increase the channel size and number of flow paths for hydrocarbons to flow from the formation 105, to the wellbore 102, and to the ground surface of the well 101. However, determining where to locate the perforation location and fracturing stages and tunnels is nontrivial. For example, formation lithology, porosity, permeability, breakdown pressure, and fracture toughness, can vary significantly with depth.

In the example of FIG. 1, the formation 105 has been fractured in the form of a fracture profile 110. In some examples, a pump 107 pumps hydraulic fluid into the wellbore and through perforations to propagate fractures in the reservoir. The hydraulic fracture profile extends vertically towards a lower portion or tip 112 and an upper portion or tip 114 and has a depth-dependent width W extending between the lower tip 112 and the upper tip 114. The fracture profile 110 has a height H defined by the vertical distance between the lower tip 112 and the upper tip 114. For example, when the pump 107 pumps hydraulic fluid into the wellbore and through perforations, a fracture will propagate from the wellbore perforation location towards a lower tip location and/or upper tip location to form the fracture in the formation. The propagation speed of the fraction in the upper and downward directions can be different depending on the stresses and fracture toughness of the formation layers.

The hydraulic fracture profile extends upwards from the perforation location 106 to the upper tip 114 by a second distance D2. The fracture profile extends downwards from the perforation location 106 to the lower tip 112 by a third distance D3. In this example, the third distance is less than the second distance. In general, the shape of the fracture above the perforation location 106 does not need to be the same as the fracture below the perforation location 106. This allows the systems and methods described herein to be able to predict fracture profiles that are asymmetric in vertical depth relative to the perforation location 106.

In some examples, the difference in upper versus lower shape of the hydraulic fracture profile is due to differences in in-situ stresses and fracture toughness of each layer of the formation. For example, if layer 104H has a higher in-situ stress and fracture toughness than layer 104C, then it will be more difficult to fracture layer 104H than it will be to fracture layer 104C, resulting in asymmetries in the fracture profile. In general, all layers (for example, layers 104A-104I) can have different in-situ stresses and fracture toughness material properties.

Wellbore system 100 includes a data processing system 120 configured predict the hydraulic fracture profile 110 in the multi-layered formation 105 by determining locations or depths of the lower tip 112 and the upper tip 114 (for example, by determining the distances D1-D3), the fracture height H, and the depth-dependent width W. In some implementations, the data processing system 120 includes a computer or computer system 500 as described with reference to FIG. 11. The data processing system 120 is configured to receive one or more inputs (for example, data of the formation 105), predict the hydraulic fracture profile 110, and, in some implementations, fracture the formation 105 using hydraulic fracturing equipment so that the formation 105 is fractured substantially similarly as the predicted hydraulic fracture profile 110. In some examples, “substantially similarly” means within 10% of the predicted values (for example, D1-D3, W, and H). In some examples, “substantially similarly” means within 5% of the predicted values (for example, D1-D3, W, and H). In some examples, “substantially similarly” means within 2% of the predicted values (for example, D1-D3, W, and H).

In some implementations, the data processing system 120, the logging device 103, and the pump 107 are configured to communicate with each other (for example, by transmitting and receiving data using one or more transceivers). In some implementations, the pump 107 is located proximate the wellbore for being in fluid communication with the wellbore 102A for injecting hydraulic fluid into the formation 105. In some implementations, the data processing system 120 is at a location different than the wellbore (for example, at a data analysis facility, or at a research facility). In some examples, the data processing system 120 is configured to be in communication with the logging device 103 and the pump 107 using a telecommunications network.

Method to Predict Hydraulic Fracture Height in Multi-Layered Formations

FIG. 2 is a plot of a hydraulic fracture propagating in a multi-layered formation 105. FIG. 2 represents a schematic of FIG. 1. Generally, reservoirs are characterized by multiple layers, which have different mechanical properties due to the diagenetic process and in-situ stress profile.

Many factors affect hydraulic fracture height evolution in multi-layered formations 105. Important parameters are minimum horizontal in-situ stress σ_i, fracture toughness K_ICi, layer thickness h_i, weak layer interfaces, and fracturing fluid leak-off. The index i denotes the i^th formation layer, which is penetrated by the hydraulic fracture. Fluid injection through the perforation needs to break down the rock first and then propagates the hydraulic fracture within the multi-layered formation 104.

The systems and methods disclosed herein obtain the full height map versus hydrostatic pressure for a given hydraulic fracture growing in a multi-layered formation 105, which utilizes the linear elastic fracture mechanics classical solution for non-uniformly loaded fracture surfaces in an opening mode. For example, the data processing system 120 is configured to execute software to obtain the full height map versus hydrostatic pressure for a given hydraulic fracture growing in a multi-layered formation 105 using the methods described herein. In subsurface, hydraulic fractures propagate laterally along the preferred far-field maximum horizontal stress direction in an infinite multi-layered formation 105. Therefore, a pressurized vertical fracture in plane-strain condition is reasonably assumed.

As shown in FIG. 2, the fracture extends from −c to +c along the vertical z-axis in an infinite domain. h_perf denotes the distance from lower fracture tip to the perforation location or point 106, P_perf denotes the fluid pressure at the perforation location 106, d_fc denotes the distance from lower fracture tip to the fracture center, and P_fc represents the fluid pressure at the fracture center.

The z-axis is assumed in the opposite direction of depth (for example, towards the ground surface of the formation). Assume the fluid pressure P_perf at the perforation location is known, which can be calculated using Bernoulli's equation. The fluid pressure within the open hydraulic fracture is assumed to be a linear distribution, which is a function of fluid density.

FIG. 3 is a plot of perforation pressure and hydrostatic pressure distribution along the fracture height and shows the fluid pressure distribution within the open fracture, which is assumed to be linear in this case. FIG. 3 represents a schematic of FIG. 1. With this assumption, the pressure at the fracture center can be calculated by:

P fc = P perf + ρ ⁢ g ⁡ ( h perf - d fc ) ( 1 )

where ρ is the fracturing fluid density, and g is gravitational acceleration (for example, 9.81 m/s²). The hydrostatic pressure within the open hydraulic fracture can be estimated by:

P ⁡ ( z ) = P fc - pgz ( 2 )

where z is the coordinate of a point along the vertical fracture height direction with respect to the z-axis origin at the fracture center.

Hydraulic fracture growth criteria is fracture toughness-based. The hydraulic fracture propagates when the stress intensity factor K_I around the fracture tip is larger than the fracture toughness K_ICi of the hosting layer. The consensus of the so-called “equilibrium height” is the hydraulic fracture height when the stress intensity factors at the upper and lower fracture tips 114, 112 are equal to the fracture toughness of the hosting layers at the same time. The hydraulic fracture is tensile fracture, which is classified as mode I type. The stress intensity factor for the lower and upper tensile fracture tips 112, 114 is defined by:

K I UFT = 1 π ⁢ c ⁢ ∫ - c c P net ( z ) ⁢ c + z c - z ⁢ dz ( 3 ) K I LFT = 1 π ⁢ c ⁢ ∫ - c c P net ( z ) ⁢ c - z c + z ⁢ dz ( 4 )

where P_net(z) is the net-pressure profile acting on the open hydraulic fracture surfaces, and c is half of the fracture height. The net pressure is the driving pressure for fracture growth and can be expressed as:

P net ( z ) = P perf + ρ ⁢ g ⁡ ( h perf - d fc ) - ρ ⁢ gz - σ i ( z ) ( 5 )

The minimum horizontal stress represents a resistance for further fracture openings. Unlike the assumptions used in analytical models, the minimum horizontal stress σ_i(z) does not need to be piece-wise constant in the numerical approaches described herein. The minimum horizontal stress σ_i(z) can be non-constant or a non-linear distribution. The systems and methods disclosed herein can account for both in-situ stress and hydrostatic pressure being non-constant.

The fracture width for any point along the fracture height can be expressed using different coordinate system as follows:

W ⁡ ( h , x i ) = 4 E ′ ⁢ π ⁢ ∫ 0 h [ P perf + ρ ⁢ g ⁡ ( h perf - x ) - σ i ( x ) ] × cosh - 1 [ y i ( h - 2 ⁢ x h ) + x ❘ "\[LeftBracketingBar]" y i - x ❘ "\[RightBracketingBar]" ] ⁢ dx ( 6 )

Making the following substitutions:

h = 2 ⁢ c , x = z + c ( 7 )

results in integral intervals of x∈[0, 2c] and z∈[−c, c], respectively. Using this integral transformation, then Eq. (6) is changed to the following form:

W ⁡ ( h , z i ) = 4 E ′ ⁢ π ⁢ ∫ - c c P net ( z ) × cosh - 1 [ ( z i + c ) ⁢ ( h - 2 ⁢ z - 2 ⁢ c 2 ⁢ c ) + z + c ❘ "\[LeftBracketingBar]" z i - z ❘ "\[RightBracketingBar]" ] ⁢ dz ( 8 ) W ⁡ ( h , z i ) = 4 E ′ ⁢ π ⁢ ∫ - c c P net ( z ) × cosh - 1 [ - z i ⁢ z c + c ❘ "\[LeftBracketingBar]" z i - z ❘ "\[RightBracketingBar]" ] ⁢ dz E ′ = E 1 - v 2 ( 9 )

where E′ is the plane strain Young's modulus for a plane strain problem and is one of the most important rock parameters for predicting hydraulic fracture geometry. E′ controls the fracture width and the value of the net pressure. E is the Young's modulus and v is the Poisson's ratio, which needs to be identified in each layer of the multi-layered formation 105. In order to calculate the fracture width profile, the systems and methods disclosed herein first loop over different values of z_i, which ranges from [0, h] and x=z+c at a coordinate point.

To find analytical solutions to Eqs. (3), (4), and (8), some simplifications are taken by assuming each formation layer has constant in-situ stress and fluid pressure. However, in the systems and methods of this disclosure, we do not limit the in-situ stress and fluid pressure within each layer to be piecewise constants. Instead, we take a numerical approach utilizing high order one-dimensional elements to discretize the potential vertical fracture height from the lower fracture tip to the upper fracture tip.

FIG. 4 is a plot of a discretized fracture from the lower tip 112 to the higher tip 114. FIG. 4 represents a schematic of FIG. 1. FIG. 4 shows one-dimensional elements 130 that discretize the potential vertical fracture height from the lower fracture tip 112 to the upper fracture tip 114. The elements used in this example are the same ones used in conventional finite element methods. While one element 130 is shown per layer, in practice multiple elements 130 are used per layer depending on the layer depth. For example, more elements can be used in layers with a larger depth and less elements can be used in layers with a smaller depth. In some examples, one element spans between 0.3 m (1 ft) and 3 m (10 ft) of the layer depth.

The accuracy of the results can be improved by increasing the order of the interpolation functions, which is also called shape functions and used to determine the values of the variables within an element by interpolating the nodal values. Polynomial type functions have been most widely used as they can be integrated, and differentiated, and the accuracy of the results can be improved by increasing the order of the polynomial. In practice, the systems and methods disclosed herein use polynomials of finite order as an approximation. Quadratic or higher order clements (for example, cubic (third order) elements, quartic (fourth order) elements, quintic (fifth order) elements etc.) can give a sufficiently high accuracy that is generally of higher accuracy than linear (first order) elements. Importantly, the in-situ stress and fluid pressure can be non-constant when calculating the equilibrium fracture height such that real world loading conditions acting on the fracture surfaces can be accounted for.

FIG. 5 is a plot of the example one dimensional second order element or quadratic element 130 shown in FIG. 4, which is merely an example of possible elements that can be used with the systems and methods disclosed herein. The element 130 has three nodes and the interpolation function is:

{ N 1 = 1 2 ⁢ ξ ⁡ ( ξ - 1 ) N 2 = 1 - ξ 2 N 3 = 1 2 ⁢ ξ ⁡ ( ξ + 1 ) ( 10 )

The vertical coordinate, net pressure, and in-situ stress, respectively, for each one-dimensional element 130 can be expresses as:

z ⁡ ( ξ ) = ∑ k = 1 N ⁢ o ⁢ d ⁢ e z k ⁢ N k ( ξ ) = z 1 ⁢ N 1 ( ξ ) + z 2 ⁢ N 2 ( ξ ) + z 3 ⁢ N 3 ( ξ ) ( 11 ) P net ( z ) = ∑ k = 1 N ⁢ o ⁢ d ⁢ e P net k ⁢ N k ( ξ ) = P net 1 ⁢ N 1 ( ξ ) + P net 2 ⁢ N 2 ( ξ ) + P net 3 ⁢ N 3 ( ξ ) ( 12 ) σ i ( z ) = ∑ k = 1 N ⁢ o ⁢ d ⁢ e σ i k ⁢ N k ( ξ ) = σ net 1 ⁢ N 1 ( ξ ) + σ net 2 ⁢ N 2 ( ξ ) + σ net 3 ⁢ N 3 ( ξ ) ( 13 )

A Gauss integration scheme is an efficient method to perform numerical integration over intervals and produces accurate results when the function to be integrated is polynomial of appropriate degree. In this disclosure, we use Gauss integration to calculate the stress intensity factors for both the lower and upper fracture tips 112, 114 as follows:

K I UFT = 1 π ⁢ c ⁢ ∫ - c c P net ( z ) ⁢ c + z c - z ⁢ dz = ∑ IEL = 1 NEL ( 1 π ⁢ c ⁢ ∑ j = 1 NG ❘ "\[LeftBracketingBar]" J j ❘ "\[RightBracketingBar]" ⁢ w j ⁢ P net ( z j ) ⁢   c + z ⁡ ( ξ j ) c - z ⁡ ( ξ j ) ) ( 14 ) K I LFT = 1 π ⁢ c ⁢ ∫ - c c P net ( z ) ⁢ c - z c + z ⁢ dz = ∑ IEL = 1 NEL ( 1 π ⁢ c ⁢ ∑ j = 1 NG ❘ "\[LeftBracketingBar]" J j ❘ "\[RightBracketingBar]" ⁢ w j ⁢ P net ( z j ) ⁢   c - z ⁡ ( ξ j ) c + z ⁡ ( ξ j ) ) ( 15 )

where NEL is the total number of one-dimensional elements used along the vertical fracture height, NG represents the number of Gauss integration points for each one-dimensional element, w_j is the weight at the Gauss integration point, and J_j is the Jacobian matrix which represents a scaling factor between the derivatives of the natural and the local coordinate systems. To improve the accuracy, we can use three or four integration points. For a one-dimensional element, the matrix J_j is given by:

J j = [ ∑ k = 1 Node ∂ N k ∂ ξ ⁢ z k ] = [ ∂ N 1 ∂ ξ ⁢ z 1 + ∂ N 2 ∂ ξ ⁢ z 2 + ∂ N 3 ∂ ξ ⁢ z 3 ] ( 16 )

The determinant of the Jacobian matrix J is given by:

❘ "\[LeftBracketingBar]" J j ❘ "\[RightBracketingBar]" = ∂ N 1 ∂ ξ ⁢ z 1 + ∂ N 2 ∂ ξ ⁢ z 2 + ∂ N 3 ∂ ξ ⁢ z 3 ( 17 )

The differences between the fracture tip stress intensity factors and fracture toughness for the upper and lower fracture tips can be calculated using Eq. 18:

{ Δ ⁢ K U = K Ii UFT - K IC = 0 Δ ⁢ K L = K Ii LFT - K IC = 0 ( 18 )

Solving the two nonlinear equations of Eq. (18) can give admissible equilibrium fracture tip locations for a corresponding perforation pressure. The top equation ΔK_U is for the upper fracture tip 114 propagation and the bottom equation ΔK_L is for the lower fracture tip 112 propagation. One or both fracture tips will continue to grow if the associated differences of ΔK_U and ΔK_L are positive. Eq. (18) can calculate the scenario that one fracture tip is stable while the other fracture tip grows and vice versa. The equilibrium fracture height is dependent on the comprehensive effect of perforation pressure, fluid density, the in-situ stresses, fracture toughness barriers, and landing depth.

Once a potential fracture height and location are determined, the fracture width profile along the vertical height can be calculated similarly using Gauss quadrature. It can be given by:

W ⁡ ( h , z i ) = ∑ IEL = 1 NEL ∑ j = 1 NG 4 E ′ ⁢ π ⁢ ❘ "\[LeftBracketingBar]" J j ❘ "\[RightBracketingBar]" ⁢ w j ⁢ P net ( z ⁡ ( ξ j ) ) × cosh - 1 ( - z i ⁢ z ⁡ ( ξ j ) c + c ❘ "\[LeftBracketingBar]" z i - z ⁡ ( ξ j ) ❘ "\[RightBracketingBar]" ) ( 19 ) cosh - 1 ⁢ x = In ⁡ ( x + x 2 - 1 ) ,   x ≥ 1 ( 20 )

where z_i can be any point elevation measured from the fracture lower tip 112 along the vertical fracture height direction where the fracture width requires calculating. And it should be fixed when using Eq. 18 to calculate the fracture width at a specified point y_i. z_i can be defined by:

z i = z ⁡ ( ξ ) + c = ∑ k = 1 Node z k ⁢ N k ( ξ ) + c = z 1 ⁢ N 1 ( ξ ) + z 2 ⁢ N 2 ( ξ ) + z 3 ⁢ N 3 ( ξ ) + c ( 21 )

Computational Workflow to Calculate Hydraulic Fracture Profile

The fracture height should be controlled within the desired pay zones as much as possible. The fluid pressure at the perforation should be well controlled, otherwise hydraulic fracture can penetrate into the overburden and underlying zones. To cover a large fracture height range, we choose to start from a small fracture with the lower and the upper tips within the same hosting layer, where perforations are created. Within this workflow, we ramp up the fluid pressure at the perforation location to predict the formation layers that the fracture upper and lower tips reach. After running a series of different fluid pressure for the perforation location, the systems and methods visualize the different required fluid pressures to break through the targeted formation layers, which might represent the pay zone, overburden and underlying zone. The results can indicate the maximum bottom hole pressure at the perforation to be controlled for hydraulic fracturing treatment.

The systems and methods predict hydraulic fracture height growth within the multi-layered formation 105 versus the perforation pressure. In this section, a strategy to obtain the fracture upper and lower tip locations 114, 112 and the perforation pressure at the fracture mouth (for example at the perforation location 106) are disclosed. The fracture height is determined by the fracture upper and lower tip locations 114, 112. These three parameters specify the predicted equilibrium state and height of a hydraulic fracture propagating in a multi-layered formation 105.

FIG. 6 is a flowchart 150 of a method to predict a hydraulic fracture profile for a perforation location. In some examples, the method predicts fracture upper and lower tip locations 114, 112 versus perforation pressure. For example, the data processing system 120 is configured to execute software to predict the fracture upper and lower tip locations 114, 112 versus perforation pressure using the method according to flowchart 150.

At block 152, the data processing system 120 collects data. In some examples, the data includes a wellbore trajectory of a wellbore, well log data of the wellbore, drilling data of the wellbore, formation tops, and height of each formation layer. In some examples, this data is used to prepare a one-dimensional mechanical earth model of the formation. In some examples, the one-dimensional mechanical earth model contains the mechanical properties and in-situ stresses along the well trajectory. For example, the logging device 103 is lowered into the wellbore 102 and is used to acquire a Young's modulus, a Poisson's ratio, and a layer thickness for each layer of the formation 105. In some examples, this data is used to determine a fracture toughness of each layer.

At block 154, the data processing system 120 divides the multi-layered formation 105 into one or more layers. In some examples, the formation layers labeled from the bottom of the formation to the top of the formation based on minimum horizontal stresses and fracture toughness for each layer.

At block 156, the data processing system 120 iterates (or loops) over a range of perforation pressures from a minimum perforation pressure to a maximum perforation pressure. For example, referring to FIG. 7A, the minimum pressure can be 39.3 MPa (5700 psi) and the maximum pressure can be 48.3 MPa (7000 psi). Generally, the fluid pressure should start from a smaller number at the beginning, otherwise it can lead to no stable fracture tip locations, which means one or both tips propagating uncontrolled within the multilayer formation.

At block 158, the data processing system 120 initializes variables i and j to correspond to the respective formation layer of interest (for example, i=4, j=4 where layer 4 is the layer containing the perforation location 106.

At block 160, the data processing system 120 generates the numerical model using high order one dimensional elements (for example, second order elements) to discretize the trial interval from the bottom of the TVD of the j^th trial stress layer to the top TVD of the i^th trial stress layer. For example, the data processing system 120 dynamically discretizes each layer along the vertical direction using one-dimensional high-order elements upon each iteration (or call) of block 160. In some examples, however, the data processing system 120 discretizes all layers initially before block 156.

In oil and gas geoscience, the formation can be divided into different formation layers first based on sedimentation time and rock facies. Each formation layer can range from several feet to several hundred feet. To get reliable solutions, the mesh throughout the layers should be fine enough. In some examples, “fine enough” means that one element spans between 0.3 m (1 ft) and 3 m (10 ft) of the layer depth. In some examples, each layer of the formation 105 is referred to as a stress layer.

At block 162, the data processing system 120 determines the potential fracture center and fracture height, and the pressure equation along the potential fracture height. For example, the data processing system 120 evaluates Eq. (1) to determine the pressure equation along the potential fracture height. For example, the data processing system 120 determines the potential fracture center and fracture height based on the trial stress layer j^th to the trial stress layer i^th and based on the perforation TVD. In some examples, the data processing system 120 determines the fracture height based on the upper i^th layer and lower j^th layer.

At block 164, the data processing system 120 determines the stress intensity factors for the fracture upper tip

K I UFT

and lower tip

K I UFT

using Gauss integration. For example, the stress intensity factor for the fracture upper tip K_I^UFT is determined using Eq. (14) and the stress intensity factor for the fracture lower tip

K I UFT

is determined using Eq. (15).

At block 166, the data processing system 120 checks (for example, determines) whether the fracture upper tip will be stable or not in the trial formation layer based on a failure criterion. For example, the failure criterion (for example, the upper fracture tip propagation condition) can be when the fracture upper tip

K I UFT

exceeds the fracture toughness K_ICi in the trial formation layer. If this criterion is met, then the layer number is incremented (for example, i=i+1) assuming that the max number of layers N has not been reached and the method goes back to block 160 to evaluate a different formation layer. In some implementations, the data processing system 120 determines whether the fracture upper tip will be stable or not in the trial formation layer is in accordance with the criterion of Eq. (18). For example, in accordance with Eq. (18), if the upper tip stress intensity factor is larger than the fracture toughness, the upper tip fracture will continue to propagate upward.

At block 168, the data processing system 120 checks (for example, determines) whether the fracture lower tip will be stable or not in the trial formation layer based on a failure criterion. In some examples, the failure criterion (for example, the lower fracture tip propagation condition) is when the fracture lower tip

K I UFT

exceeds the fracture toughness K_ICi in the trial formation layer. If this criterion is met, then the data processing system 120 decrements the layer number (for example, j=j−1) assuming that layer 1 has not been reached and the data processing system 120 goes back to (e.g., or reverts to) block 160 to evaluate a different formation layer. In some implementations, the data processing system 120 determines whether the fracture lower tip will be stable or not in the trial formation layer is in accordance with the criterion of Eq. (18). For example, in accordance with Eq. (18), if the lower tip stress intensity factor is larger than the fracture toughness, the lower tip fracture will continue to propagate downward.

When the data processing system 120 determines that the failure criteria of blocks 166 and 168 are not met (for example, denoting that the fracture does not propagate completely through the trial formation layers), the data processing system 120 proceeds to block 170.

At block 170, the data processing system 120 calculates (for example, determines) the fracture height and fracture center for the fracture profile based on the final upper stress layer's top TVD and the bottom stress layer's bottom TVD. For example, the data processing system 120 determines an admissible fracture upper tip location (for example, the top TVD at which the fracture upper tip currently reaches) and the lower tip location (for example, the bottom TVD at which the fracture lower tip currently reaches), when both ΔK_U≤0 and ΔK_L≤0, for the upper and lower fractures respectively, are met.

At block 172, the data processing system 120 calculates (for example, determines) the fracture width profile along the stable fracture height. In some implementations, the data processing system 120 evaluates Eqs. (19), (20), and (21) to determine the fracture width profile.

In some examples, the data processing system 120 optimizes a hydraulic fracturing pump schedule design to control the downhole pressure based on the targeted fracture height so that hydraulic fracture height is contained within the production zone. For example, the data processing system 120 determines a hydraulic fracturing pump schedule based on the fracture height based on the lower fracture tip location and the upper fracture tip location. In some examples, the data processing system 120 selects (for example, determines) a proppant size and type for admittance into fracture surface to guide how much of the hydraulic fracture width created by the fracturing treatment will be retained as a propped width once the hydraulic fractures closes. In some examples, the data processing system 120 determines a location of the well based on the fracture height. For example, the data processing system 120 can determine where the horizontal well should be landed in the reservoirs for good production based on the targeted fracture height.

Example #1

This section includes an example for simulating fracture height growth which initiates in a low stress layer. Table 1 provides the data for the example study which includes the layer thicknesses of the respective layer of the formation, minimum horizontal in-situ stresses exerting on the respective layer as a function of depth within the respective layer (assumed constant for this example), and fracture toughness for each layer. The layers are numbered subsequently from the bottom #1 to the top #7 and can correspond to layers 104I-104A of FIG. 1. The perforation location (which can correspond to perforation location 106) is given at TVD of 2821 m (9255 ft), in which the layer has the lowest minimum horizontal stress 39.3 MPa (5700 psi). In some implementations, the data processing system 120 acquires the data represented in Table 1 from wellbore logging device 103 and stores the data in memory for retrieval when the data processing system 120 executes software to predict fracture profiles.

TABLE 1
Formation layers with piecewise constant in-situ stresses.

Layer	Top TVD	Bottom TVD	Top Layer	Bottom Layer	Fracture toughness
No.	(ft)	(ft)	Shmin (psi)	Shmin (psi)	(psi√{square root over (in)})	Perforation

7	8700	8990	6500	6500	2000	False
6	8990	9080	7150	7150	2000	False
5	9080	9170	7150	7150	2000	False
4	9170	9340	5700	5700	2000	True
3	9340	9380	7350	7350	2000	False
2	9380	9455	5800	5800	2000	False
1	9455	9650	8200	8200	2000	False

FIG. 7A includes a plot 200A of the horizontal in-situ stress distribution across the formation. As shown, the horizontal in-situ stresses are constant within each layer in this example. Plot 200A also indicates the perforation location 106 of the wellbore. FIG. 7A also includes a plot 200B of a fracture height map versus perforation pressure according to workflow 150. The equilibrium height map includes the fracture upper-and lower-tip locations 114, 112 vs. perforation pressure. In some implementations, the data processing system 120 generates plots 200A and 200B based on the data from Table 1 after predicting one or more fracture profiles.

A first curve 206 represents a dependence of a first vertical depth (for example, the depth of the upper tip 114) on perforation pressure and is generally different than a second curve 208 representing a dependence of the second vertical depth (for example, the depth of the lower tip 112) on perforation pressure. As shown in box 202, the progression of the lower tip 112 drops abruptly from layer #4 through layer #3 at approximately 45.12 MPa (6,550 psi). In this example, layer #2 has much lower in-situ stress of 40.0 MPa (5,800 psi) compared with 50.7 MPa (7,350 psi) and 56.5 MPa (8,200 psi) in the neighboring layers.

A jump in fracture height occurs when the lower tip 112 completely crosses layer #3 and is about to enter layer #2 that is a low-stress layer with the same toughness as other layers. This reflects the fact that fracture propagation much easy in a lower layer of lower propagation resistance after the lower tip 112 reaches the layer #2. For example, the fracture can quickly propagate to the bottom of layer #2. Vertical lines typically mean low in-situ stress and fracture toughness of formation such that when the specific pressure is applied, fracture propagates very easily. Hence, the more vertical a line appears, the easier the fracture is to propagate through the respective formation layer. For example, a pressure of 45.2 MPa (6550 psi) is a critical pressure which enables the fracture to penetrate to the middle of in-situ stress layer 2 (not half of layer 2, full layer 2). The data processing system 120 determines whether equilibrium has been achieved for each admissible fracture tip locations based on the criteria.

The data processing system 120 uses the profile of the fracture tip locations versus fracture height to generate plot 200B to visualize the contribution of each fracture tip 112, 114 to the fracture height growth (for example, the growth of a vertical distance between the lower fracture tip 112 and the upper fracture tip 114). FIG. 7A indicates that the slope of fracture height growth over bottom hole pressure is highly dependent on each layer minimum horizontal stress. As evident in FIG. 7A, higher minimum horizontal stress leads to a slower fracture height growth as pressure is increased. In some implementations, the data processing system 120 uses the height maps to define boundaries on the vertical depths of a modeled height of the wellbore in complex three-dimensional models based on a specific treating pressure. For example, the data processing system 120 can use the height map as a maximum value (or limit) for the hydraulic fracture height to propagate from a wellbore for a specific treating pressure. In some implementations, the data processing system 120 plots horizontal lines 204 representing the top and bottom of each layer of the formation 105 to aid in visualizing the locations of the upper and lower tips 114, 112. In some implementations, the data processing system 120 plots horizontal lines 204 representing the formation layers in both plots 200A and 200B.

Example #2

The systems and methods disclosed herein can also account for a complex distribution of in-situ stress in each formation layer and fluid pressure profiles in hydraulic fractures. Example #2 assumes non-constant in-situ stress for each formation layer according to Table 2. In some implementations, the data processing system 120 acquires the data represented in Table 1 from wellbore logging device 103 and stores the data in memory for retrieval when the data processing system 120 executes software to predict fracture profiles.

TABLE 2
Formation layers with non-constant minimum in-situ stresses.

No. of	Top TVD	Bottom	Top Layer Shmin	Bottom Layer Shmin	Fracture toughness
Layer	(ft)	TVD (ft)	(psi)	(psi)	(psi√{square root over (in)})	Perforation

6	8990	9080	7050	7100	2000	False
5	9080	9170	7100	7150	2000	False
4	9170	9340	5600	5700	2000	True
3	9340	9380	7300	7350	2000	False
2	9380	9455	5800	5900	2000	False
1	9455	9650	8100	8200	2000	False

FIG. 7B includes a plot 220A of the horizontal in-situ stress distribution across the formation. As shown, the horizontal in-situ stresses are not constant within each layer in this example. Plot 220A also indicates the perforation location 106 of the wellbore. FIG. 7B also includes a plot 220B of a fracture height map versus perforation pressure according to workflow 150. The equilibrium height map includes the fracture upper-and lower-tip locations 114, 112 vs. perforation treating pressure. In some implementations, the data processing system 120 generates plots 220A and 220B based on the data from Table 2 after predicting one or more fracture profiles.

A first curve 222 represents a dependence of a first vertical depth (for example, the depth of the upper tip 114) on perforation pressure and is generally different than a second curve 224 representing a dependence of the second vertical depth (for example, the depth of the lower tip 112) on perforation pressure. Further, there is a difference between curves 222, 224 and curves 206, 208 due to the data processing system 120 accounting for the depth-dependent horizontal in-situ stresses in the fracture profile prediction.

In plot 220A, the non-constant in-situ stresses are visualized as a function of depth. This is in contrast to plot 200A, which shows straight vertical lines representing constant in-situ stresses as a function of depth. Non-constant hydrostatic pressure is also accounted for in the examples of FIG. 7A and 7B using interpolation and gauss quadrature.

By including a plurality of layers in the fracture profile prediction, the data processing system 120 predicts a full-height map as shown in plots 200B, 220B of FIGS. 7A and 7B. In turn, the data processing system 120 uses the full-height map versus perforation pressure to determine fracturing locations and pressure ranges used during hydraulic fracturing. For example, the data processing system 120 uses the full-height map to determine locations where the hydraulic fractures penetrate within the greatest number of layers of a target or pay zone (for example, shale and sandstone) in the multi-layered formation while avoiding layers outside the target or pay zone (for example, overburden or underlying layers above and below the target or pay zone, respectively) within practical pressure limits achievable by the pump 107 (for example, up to 48.3 MPa (7,000 psi) or 68.9 MPa (10,000 psi)).

In some implementations, the data processing system 120 uses the full-height map to optimize or determine preferable pump schedules for hydraulic fracturing design. Once the data processing system 120 determines the fracture height based on the perforation pressure, the data processing system 120 calculates (for example, determines) the fracture aperture profile (for example, according to block 172).

FIGS. 8A-8C are plots of fracture widths for perforation pressures of 41.4 MPa (6000 psi), 44.1 MPa (6400 psi), and 46.5 MPa (6750 psi), respectively. The results shown in FIGS. 8A-8C reflect the results of FIG. 7B at these specific perforation pressures. As with the results of FIGS. 7A and 7B, the data processing system 120 plots indications of the upper and lower boundaries of the layers to assist the user in identifying the location of the fracture.

FIG. 8A shows the fracture aperture at the perforation pressure of around 41.4 MPa (6000 psi), which indicates the fracture height is contained within the targeted zone. The targeted zone is defined as the reservoirs that produce oil or gas. FIG. 8B shows the fracture will breakthrough the layer #4 when the fluid pressure at the perforation point reaches 44.1 MPa (6400 psi). FIG. 8C displays the fracture aperture profile when the fluid pressure at the perforation reaches 46.5 MPa (6750 psi). As shown in FIG. 8C, the fracture aperture at the high stress layer exhibits a pinching behavior compared to the lower stress layer. Plot 200B includes horizontal lines 204 representing the top and bottom of each layer of the formation to aid in visualizing the location of the upper and lower tips 114, 112.

The data processing system 120 uses the fracture width to determine a preferable hydraulic fracturing design. When the pump 107 pumps hydraulic fluid into the formation, the fracture width will generally be dependent on the in-situ stresses, formation properties, injection rate, and fluid type, which can decrease substantially after flowback (for example, by 50% or more). In some examples, the fracture width during fracturing treatment indicates the potential fracture surface areas where the selected proppants can reach. The fracture width can also dictate the final propped fracture width after accounting for fracturing flowback. The fracture width after fracturing treatment is usually substantially smaller than the fracture width during hydraulic fracturing fluid injection. The retained fracture width can be highly dependent on the proppant coverage, proppant size, proppant type, and formation embedment strength. The post fracturing production is limited by the fracture width.

In some implementations, the data processing system 120 determines proppant coverage, proppant size, and proppant type based on the fracture width by accounting for formation embedment strength and expected fracturing flowback. In some implementations, the data processing system 120 automatically selects (for example, determines) the appropriate proppant size based on the predicted fracture width profile along the fracture height. In some implementations, the data processing system 120 determines which proppant to use based on conductivity of the predicted fracture width at in-situ stress conditions. In general, larger proppant size yields better fracture conductivity after fracturing treatment, but size should be checked against proppant admittance criteria, both through the perforations and inside the hydraulic fracture. Too large proppant size can cause early screen out, which can lead to an unsuccessful hydraulic fracturing.

FIG. 9A is an illustration of a fracture 350 without proppant and FIG. 9B is an illustration of a fracture 360 with proppant 362. In some examples, the data processing system 120 determines the diameter of the proppant 362 based on the predicted fracture profile. For example, if the fracture width is bigger based on calculation for an in-situ stress and treating pressure, the data processing system 120 determines that the proppant size should be increased relative to a baseline size, otherwise the data processing system 120 determines that the proppant size should be decreased relative to the baseline size.

Proppant placement and coverage within hydraulic fracture remains a challenging question given that proppant transport is a complicated phenomenon coupling several physics such as fracture growth, fluid flow inside fracture, fluid leak-off, and proppant settling and jamming. Proppant density should be as small as possible, which facilitates proppant transport relatively casy within fractures. Also the proppant size can be selected based on the predicted fracture width profile during hydraulic fracturing treatment. In some implementations, the data processing system 120 determines that a larger proppant size should be used when the predicted fracture width is large (for example, above a threshold). In some implementations, the data processing system 120 determines that a smaller proppant size should be used when the predicted fracture width is small (for example, below the threshold).

As an effective computer simulation tool, the systems and methods described herein can be used to optimize the pump schedule before hydraulic fracturing treatment. Pump schedules typically define the injection rate, perforation fluid pressure, proppant selection, proppant size, and proppant concentration for fracturing a wellbore. In some implementations, the data processing system 120 predicts the fracture height for various downhole pressures at the perforation based on surface treating pressures, injection rates, fluid pressure friction losses from the wellhead to the downhole perforation locations, fluid types, proppant sizes, and concentration, for example.

Examples of Fracturing a Wellbore based on a Predicted Fracture Profile

FIG. 10 is a flowchart of an example method 400 for fracturing the wellbore 102. In some implementations, the data processing system 120 includes one or more processors configured to execute software to perform the steps of the method 400.

At block 402, the data processing system 120 receives data representing a depth-dependent horizontal in-situ stress for at least two layers of the formation and a fracture toughness for the at least two layers of the formation. For example, in-situ stresses and fracture toughness can be measured or determined using one or more logs or acquired from a database for a formation of interest. In some examples, logging device 103 is lowered into the wellbore 102 to acquire data determining the in-situ stresses and fracture toughness of the layers of the formation 105. In some examples, logging device 103 acquires data described with reference to block 152 of FIG. 6.

At block 404, the data processing system 120 predicts a hydraulic fracture profile of the formation for at least one perforation pressure. In some implementations, the hydraulic fracture profile of the formation is predicted by the data processing system 120 performing blocks 404A-404F.

At block 404A, the data processing system 120 determines an upper fracture tip stress intensity factor within the formation for the at least one perforation pressure, the determining being based on the depth-dependent horizontal in-situ stress of one or more layers of the at least two layers of the formation. For example, the data processing system 120 performs the steps of blocks 154-164 described with reference to FIG. 6 to determine an upper fracture tip stress intensity factor K_I^UFT for an upper portion (for example, tip) of the predicted hydraulic fracture profile. In some examples, the upper portion is above a perforation location of a wellbore.

In some implementations, the data processing system 120 evaluates Eq. (14) to determine the upper fracture tip stress intensity factor of the upper portion. For example, in some implementations, the data processing system 120 determines the upper fracture tip stress intensity factor by evaluating:

K I UFT = ∑ IEL = 1 NEL ( 1 π ⁢ c ⁢ ∑ j = 1 NG ❘ "\[LeftBracketingBar]" J j ❘ "\[RightBracketingBar]" ⁢ w j ⁢ P net ( z j ) ⁢   c + z ⁡ ( ξ j ) c - z ⁡ ( ξ j ) )

where NEL is the total number of one-dimensional elements used along the height of the fracture,

c = h 2

where h is the height of the fracture, NG is the number of Gauss integration points for each one-dimensional element, J_j is the Jacobian matrix, w_j is the weight at the jth Gauss integration point, P_net(z_j) is the net pressure acting on open hydraulic fracture surfaces at depth z of the jth Gauss integration point, and z(ξ_j) is the depth z evaluated at ξ_j.

At block 404B, the data processing system 120 determines a first vertical depth within the wellbore at which the upper fracture tip stress intensity factor satisfies an upper fracture tip propagation condition based on the fracture toughness of one or more layers of the at least two layers of the formation. In some implementations, the data processing system 120 performs the steps of blocks 166-170 described with reference to FIG. 6 to determine a first vertical depth as the depth of the final upper stress layer's top TVD. In some examples, the data processing system 120 evaluates Eq. (18) to determine the first vertical depth.

At block 404C, the data processing system 120 determines a lower fracture tip stress intensity factor within the formation for the at least one perforation pressure, the determining being based on the depth-dependent horizontal in-situ stress of one or more layers of the at least two layers of the formation. For example, the data processing system 120 performs the steps of blocks 154-164 described with reference to FIG. 6 to determine a lower stress intensity factor KLFT for a lower portion (for example, tip) of the fracture profile. In some examples, the lower portion is below the perforation location of the wellbore.

In some implementations, the data processing system 120 evaluates Eq. (15) to determine the lower fracture tip stress intensity factor of the lower portion. For example, in some implementations, the data processing system 120 determines the lower fracture tip stress intensity factor for the lower tip of the predicted fracture profile by evaluating:

K I LFT = ∑ IEL = 1 NEL ( 1 π ⁢ c ⁢ ∑ j = 1 NG ❘ "\[LeftBracketingBar]" J j ❘ "\[RightBracketingBar]" ⁢ w j ⁢ P net ( z j ) ⁢   c - z ⁡ ( ξ j ) c + z ⁡ ( ξ j ) )

where NEL is the total number of one-dimensional elements used along the height of the fracture,

c = h 2

where h is the height of the fracture, NG is the number of Gauss integration points for each one-dimensional element, J_j is the Jacobian matrix, w_j is the weight at the jth Gauss integration point, P_net(z_j) is the net pressure acting on open hydraulic fracture surfaces at depth z of the jth Gauss integration point, and z(ξ_j) is the depth z evaluated at ξ_j.

At block 404D, the data processing system 120 determines a second vertical depth within the wellbore at which the lower fracture tip stress intensity factor satisfies a lower fracture tip propagation condition based on the fracture toughness of one or more layers of the at least two layers of the formation. In some implementations, the data processing system 120 performs the steps of blocks 166-170 described with reference to FIG. 6 to determine a second vertical depth as the depth of the final lower stress layer's bottom TVD. In some examples, the data processing system 120 evaluates Eq. (18) to determine the second vertical depth.

At block 404E, the data processing system 120 determines a height of the hydraulic fracture profile for the at least one perforation pressure, the height extending from the first vertical depth to the second vertical depth. In some implementations, the data processing system 120 performs the steps of block 170 described with reference to FIG. 6 to determine the height of the hydraulic fracture profile. In some examples, the data processing system 120 solves Eq. (18) to determine the height of the hydraulic fracture profile.

At block 404F, the data processing system 120 determines a depth-dependent width of the hydraulic fracture profile at the at least one perforation pressure, the depth-dependent width extending along the height of the hydraulic fracture profile. In some implementations, the data processing system 120 performs the steps of block 172 described with reference to FIG. 6 to determine the depth-dependent width of the hydraulic fracture profile. The predicted hydraulic fracture profile for the at least one perforation pressure is at least partially defined by the determined height and the determined depth-dependent width for the at least one perforation pressure.

In some implementations, the data processing system 120 evaluates Eqs. (19)-(21) to determine a depth-dependent width of the hydraulic fracture profile. For example, in some implementations, the data processing system 120 determines the depth-dependent width of the hydraulic fracture profile along the height of the hydraulic fracture profile by evaluating:

W ⁡ ( h , z i ) = ∑ IEL = 1 NEL ∑ j = 1 NG 4 E ′ ⁢ π ⁢ ❘ "\[LeftBracketingBar]" J j ❘ "\[RightBracketingBar]" ⁢ w j ⁢ P net ( z ⁡ ( ξ j ) ) × cosh - 1 ( - z i ⁢ z ⁡ ( ξ j ) c + c ❘ "\[LeftBracketingBar]" z i - z ⁡ ( ξ j ) ❘ "\[RightBracketingBar]" )

where h is the height of the fracture, z_i is the depth of the fracture where the depth-dependent width is to be determined, NEL is the total number of one-dimensional elements used along the height of the fracture, NG is the number of Gauss integration points for each one-dimensional clement, E′ is the plane strain Young's modulus, J_j is the Jacobian matrix, w_j is the weight at the jth Gauss integration point, P_net(z(ξ_j)) is the net pressure acting on open hydraulic fracture surfaces at depth z evaluated at ξ_j, and

c = h 2

where h is the height of the fracture.

At block 406, the data processing system 120 determines a perforation pressure value based on at least one of the predicted hydraulic fracture profiles of the formation. For example, the data processing system 120 identifies the top surface of the target zone and selects the perforation pressure in which the upper tip of the hydraulic fracture profile does not go above the top surface of the target zone. In some examples, the data processing system 120 identifies the bottom surface of the target zone and selects the perforation pressure in which the lower tip of the fracture profile does not go below the bottom surface of the target zone. In some examples, a petrophysical engineer manually selects the perforation pressure based on a desired hydraulic fracture profile. In some implementations, the data processing system 120 determines the perforation pressure value automatically without user assistance.

At block 408, a pump injects a fluid into the wellbore at the determined perforation pressure to fracture the formation. For example, the pump 107 injects a hydraulic fluid into the wellbore 102 such that the hydraulic fluid flows through the perforation locations 106 and enters the formation 105. Due to the accuracy of the systems and method described herein, it is expected that the formation 105 will fracture in a similar way as the predicted fracture profile when subjected to the perforation pressure.

In some implementations, the depth-dependent width of the hydraulic fracture profile is determined using Gauss quadrature. For example, the processor uses Eq. (19) to determine the depth-dependent width of the fracture profile.

In some implementations, the upper and lower fracture tips' stress intensity factors are determined using Gauss quadrature. For example, the processor uses Eqs. (14) and (15) to determine the upper and lower fracture tip stress intensity factors.

In some implementations, a logging device measures a Young's modulus, a Poisson's ratio, and a layer thickness for each layer of at least two layers of the formation. In such cases, the data processing system 120 predicts the hydraulic fracture profile based on the Young's modulus, the Poisson's ratio, and the layer thickness for each layer of the at least two layers of the formation. In some examples, the logging device 103 is lowered into the wellbore to measure the Young's modulus, the Poisson's ratio, and the layer thickness for each layer of at least two layers of the formation 105.

In some implementations, predicting the hydraulic fracture profile of the formation includes determining, by the processor on a layer-by-layer basis, whether the upper fracture tip stress intensity factor satisfies the upper fracture propagation condition and whether the lower fracture tip stress intensity factor satisfies the lower fracture tip propagation condition. For example, as described with reference to FIG. 6, the data processing system 120 can test each formation layer on a layer-by-layer basis (for example, layer i=1, followed by layer i=2, etc.).

In some implementations, predicting the hydraulic fracture profile of the formation includes dynamically discretizing, by the processor on the layer-by-layer basis, each layer of the at least two layers of the formation into a plurality of one dimensional second order (or higher order) finite elements. For example, as described with reference to FIG. 6, the data processing system 120 can generate a one-dimensional mesh on each layer on a layer-by-layer basis (for example, layer i=1, followed by layer i=2, etc.). In some examples, the data processing system 120 discretizes each layer using one dimensional second order (or higher order) finite elements 130 described with reference to FIGS. 4 and 5.

In some implementations, the depth-dependent width of the hydraulic fracture profile along the height of the fracture profile is determined using a plane strain Young's modulus. For example, in some implementations, the data processing system 120 uses Eq. (9) to account for a plane strain Young's modulus.

In some implementations, the at least two layers of the formation includes at least five layers of the formation. For example, the formation can include eight layers as shown in FIG. 1 or seven layers as shown in FIGS. 7A-7B and 8A-8C. In general, there is no upper limit on the number of layers that can be considered by the systems and methods disclosed herein.

In some implementations, the upper and lower fracture tip stress intensity factors are determined based on a depth-dependent hydrostatic pressure within the formation. For example, the data processing system 120 determines the upper and lower fracture tip stress intensity factors by evaluating Eqs. (14) and (15) which account for depth-dependent hydrostatic pressure within the formation.

In some implementations, the data processing system 120 plots a fracture height map on a first plot representing a progression of the height of the predicted hydraulic fracture profile as a function of the at least one perforation pressure. For example, in some implementations, the data processing system 120 plots the fracture height maps shown in FIGS. 7A and 7B. In some examples, the data processing system 120 plots the fracture height map with a range of perforation pressures on a first axis and locations of the upper and lower tips a second axis. For example, as shown in FIGS. 7A and 7B, the range of perforation pressures are on a first axis (for example, the horizontal axis), and the location within or depth of the formation is on a second axis (for example, the vertical axis). The first vertical depth (for example, corresponding to the upper tip 114) and the second vertical depth (for example, corresponding to the lower tip 112) are shown as a function of perforation pressure.

In some implementations, the data processing system 120 plots depth indications on the first plot. The depth indications represent the top and bottom boundaries for each layer of the at least two layers of the formation. For example, the data processing system 120 generates a plot with horizontal lines 204 representing the top and bottom of each layer of the formation.

In some implementations, the data processing system 120 plots a first curve representing a dependence of the first vertical depth on perforation pressure, and a second curve representing a dependence of the second vertical depth on perforation pressure. In some examples, the first and second curves are different from one another due to differences in (i) the fracture toughness of each layer of the at least two layers of the formation, (ii) the depth-dependent horizontal in-situ stresses of each layer of the at least two layers of the formation, and (iii) the depth-dependent hydrostatic pressure acting on the surfaces of the hydraulic fracture profile. For example, curve 222 is different than curve 224 of FIG. 7B and these differences are due to differences in fracture toughness, depth-dependent horizontal in-situ stress, and depth-dependent hydrostatic pressure.

In some implementations, the data processing system 120 plots the predicted hydraulic fracture profile at the perforation pressure on a second plot with the depth-dependent width of the predicted hydraulic fracture profile on a first axis and the height of the predicted fracture profile on a second axis such that a two-dimensional shape of the predicted fracture profile at the perforation pressure is visible. For example, the data processing system 120 generates plots 300B, 320B, and/or 340B shown in FIGS. 8A-8C, respectively. In such examples, the depth-dependent fracture width of the predicted fracture profile is on a first axis (for example, the horizontal axis), and the height of the predicted profile is on a second axis (for example, the vertical axis).

In some implementations, the data processing system 120 plots depth indications on the second plot. For example, the data processing system 120 generates a plot with horizontal lines 204 representing the top and bottom of each layer of the formation such that intersections between the predicted hydraulic fracture profile and the at least two layers are visible.

In another example implementation for fracturing a formation, which can be combined with the implementation described above and/or used standalone, the data processing system 120 (i) receives data representing fracture toughness properties for a plurality of layers of the formation, and depth-dependent horizontal in-situ stresses exerted on the plurality of layers of the formation, and (ii) automatically generates a plurality of one-dimensional elements along a vertical direction of a first layer of the formation. For example, the data processing system 120 automatically generates one-dimensional elements 130 spanning the entire vertical dimension of the first layer. In some examples, the first layer is layer 4 shown in FIG. 7A and represents the layer in which the perforation locations are located. In some examples, the first layer is layer 104E described with reference to FIG. 1. In either example, the data processing system 120 receives an indication that the perforation location is in the first layer and automatically generates the one-dimensional elements across the vertical dimension of the first layer. In some implementations, the data processing system 120 generates the one-dimensional elements such that a vertical length of each element is between 0.3 m (1 ft) and 3 m (10 ft) of the vertical dimension of the layer.

The data processing system 120 numerically predicts whether a fracture is expected to occur within the first layer based on (i) a perforation pressure, (ii) a fracture toughness of the layer based on the received fracture toughness properties, and (iii) a depth-dependent horizontal in-situ stress of the layer based on the received depth-dependent horizontal in-situ stresses. For example, in some implementations, the data processing system 120 evaluates one or more equations of Eqs. (1)-(18) to predicts whether a hydraulic fracture is expected to occur within the first layer. In this implementation, the data processing system 120 may ignore the geometry of the surrounding layers since those layers do not necessarily need to be meshed.

In some implementations, when the data processing system 120 determines that the hydraulic fracture is expected to occur in the first layer, the data processing system 120 controls a pump to inject a fluid into a wellbore at the perforation pressure to fracture the first layer of the formation or defines a pump schedule that includes injecting the fluid into the wellbore at the perforation pressure to fracture the first layer of the formation. For example, in this example implementation, the data processing system 120 determined the hydraulic fracture profile and perforation pressure without needing to model or mesh the surrounding layers. Thus, the data processing system 120 predicts the hydraulic fracture profile in a computationally efficient manner. In some cases, the injected fluid causes a hydraulic fracture to propagate through at least a portion of the first layer of the formation and stop within the first layer of the formation such that the hydraulic fracture does not extend into other layers of the formation beyond the first layer.

In some implementations, when the data processing system 120 determines that the hydraulic fracture is not expected to occur in the first layer, the data processing system 120 automatically generates a plurality of one-dimensional elements along the vertical direction of a second layer of the formation and numerically predicts whether a hydraulic fracture is expected to occur within the second layer. For example, the data processing system 120 generates the one-dimensional elements on the second layer after the data processing system 120 predicts that the hydraulic fracture profile is likely to extend beyond (for example, either above or below) the first layer. Thus, the data processing system 120 automatically updates the meshed domain for the model when the data processing system 120 determines that a larger domain is required to predict the fracture profile.

In some implementations, when the data processing system 120 determines that the hydraulic fracture is expected to occur in the second layer, the data processing system 120 controls the pump to inject the fluid into the wellbore at a perforation pressure at which the hydraulic fracture is expected to occur in the second layer to fracture the second layer of the formation or defines the pump schedule such that the pump schedule includes pumping the fluid into the wellbore at the perforation pressure to fracture the second layer of the formation. For example, in this example implementation, the data processing system 120 determined the hydraulic fracture profile and perforation pressure without needing to model or mesh the layers other than the first and second layers. In some cases, the injected fluid causes a hydraulic fracture to propagate through the first layer of the formation and at least a portion of the second layer and stop within the second layer of the formation such that the hydraulic fracture does not extend into other layers of the formation beyond the first and second layers.

In some implementations, when the data processing system 120 determines that the hydraulic fracture is not expected to occur in the second layer, the data processing system 120 automatically generates a plurality of one-dimensional elements along the vertical direction of a third layer of the formation and numerically predicts whether a hydraulic fracture is expected to occur within the third layer. For example, the data processing system 120 generates one-dimensional elements on the third layer after the data processing system 120 predicts that the hydraulic fracture profile is likely to extend beyond the first and second layers.

In some implementations, when the data processing system 120 determines that the hydraulic fracture is expected to occur in the third layer, the data processing system 120 controls the pump to inject the fluid into the wellbore at a perforation pressure at which the hydraulic fracture is expected to occur in the third layer to fracture the third layer of the formation or defines the pump schedule such that the pump schedule includes injection the fluid into the wellbore at the perforation pressure to fracture the third layer of the formation. For example, in this example implementation, the data processing system 120 determined the hydraulic fracture profile and perforation pressure without needing to model or mesh the layers other than the first, second, and third layers. In some implementations, this process continues (for example, by considering the fourth layer, fifth layer, etc.) until the fracture profile is predicted. In some cases, the injected fluid causes a hydraulic fracture to propagate through the first and second layers of the formation and at least a portion of the third layer and stop within the third layer of the formation such that the hydraulic fracture does not extend into other layers of the formation beyond the first, second, and third layers.

Systems and methods configured to predict hydraulic fracture height and fracture width profile for fracturing wells in multi-layered formations are disclosed. These systems and methods are applicable to real in-situ stress distributions and fluid pressure distributions in hydraulic fractures. The systems and methods do not require simplification of in-situ stress distributions and/or fluid pressure distributions for the formation layers.

FIG. 11 is a schematic illustration of an example controller 500 (or control system) for predicting hydraulic fracture profiles according to the present disclosure. For example, the controller 500 can include or be part of the data processing system 120. The controller 500 is intended to include various forms of digital computers, such as printed circuit boards (PCB), processors, digital circuitry, or otherwise parts of a system for determining a subsurface formation breakdown pressure. Additionally the system can include portable storage media, such as, Universal Serial Bus (USB) flash drives. For example, the USB flash drives may store operating systems and other applications. The USB flash drives can include input/output components, such as a wireless transmitter or USB connector that may be inserted into a USB port of another computing device.

The controller 500 includes a processor 510, a memory 520, a storage device 530, and an input/output device 540 (for displays, input devices, example, sensors, valves, pumps). Each of the components 510, 520, 530, and 540 are interconnected using a system bus 550. The processor 510 is capable of processing instructions for execution within the controller 500. The processor may be designed using any of a number of architectures. For example, the processor 510 may be a CISC (Complex Instruction Set Computers) processor, a RISC (Reduced Instruction Set Computer) processor, or a MISC (Minimal Instruction Set Computer) processor.

In one implementation, the processor 510 is a single-threaded processor. In another implementation, the processor 510 is a multi-threaded processor. The processor 510 is capable of processing instructions stored in the memory 520 or on the storage device 530 to display graphical information for a user interface on the input/output device 540.

The memory 520 stores information within the controller 500. In one implementation, the memory 520 is a computer-readable medium. In one implementation, the memory 520 is a volatile memory unit. In another implementation, the memory 520 is a non-volatile memory unit.

The storage device 530 is capable of providing mass storage for the controller 500. In one implementation, the storage device 530 is a computer-readable medium. In various different implementations, the storage device 530 may be a floppy disk device, a hard disk device, an optical disk device, or a tape device.

The input/output device 540 provides input/output operations for the controller 500. In one implementation, the input/output device 540 includes a keyboard and/or pointing device. In another implementation, the input/output device 540 includes a display unit for displaying graphical user interfaces.

The features described can be implemented in digital electronic circuitry, or in computer hardware, firmware, software, or in combinations of them. The apparatus can be implemented in a computer program product tangibly embodied in an information carrier, for example, in a machine-readable storage device for execution by a programmable processor; and method steps can be performed by a programmable processor executing a program of instructions to perform functions of the described implementations by operating on input data and generating output. The described features can be implemented advantageously in one or more computer programs that are executable on a programmable system including at least one programmable processor coupled to receive data and instructions from, and to transmit data and instructions to, a data storage system, at least one input device, and at least one output device. A computer program is a set of instructions that can be used, directly or indirectly, in a computer to perform a certain activity or bring about a certain result. A computer program can be written in any form of programming language, including compiled or interpreted languages, and it can be deployed in any form, including as a stand-alone program or as a module, component, subroutine, or other unit suitable for use in a computing environment.

Suitable processors for the execution of a program of instructions include, by way of example, both general and special purpose microprocessors, and the sole processor or one of multiple processors of any kind of computer. Generally, a processor will receive instructions and data from a read-only memory or a random access memory or both. The essential elements of a computer are a processor for executing instructions and one or more memories for storing instructions and data. Generally, a computer will also include, or be operatively coupled to communicate with, one or more mass storage devices for storing data files; such devices include magnetic disks, such as internal hard disks and removable disks; magneto-optical disks; and optical disks. Storage devices suitable for tangibly embodying computer program instructions and data include all forms of non-volatile memory, including by way of example semiconductor memory devices, such as EPROM, EEPROM, and flash memory devices; magnetic disks such as internal hard disks and removable disks; magneto-optical disks; and CD-ROM and DVD-ROM disks. The processor and the memory can be supplemented by, or incorporated in, ASICs (application-specific integrated circuits).

To provide for interaction with a user, the features can be implemented on a computer having a display device such as a CRT (cathode ray tube) or LCD (liquid crystal display) monitor for displaying information to the user and a keyboard and a pointing device such as a mouse or a trackball by which the user can provide input to the computer. Additionally, such activities can be implemented via touchscreen flat-panel displays and other appropriate mechanisms.

The features can be implemented in a control system that includes a back-end component, such as a data server, or that includes a middleware component, such as an application server or an Internet server, or that includes a front-end component, such as a client computer having a graphical user interface or an Internet browser, or any combination of them. The components of the system can be connected by any form or medium of digital data communication such as a communication network. Examples of communication networks include a local area network (“LAN”), a wide area network (“WAN”), peer-to-peer networks (having ad-hoc or static members), grid computing infrastructures, and the Internet.

While this specification contains many specific implementation details, these should not be construed as limitations on the scope of any inventions or of what may be claimed, but rather as descriptions of features specific to particular implementations of particular inventions. Certain features that are described in this specification in the context of separate implementations can also be implemented in combination in a single implementation. Conversely, various features that are described in the context of a single implementation can also be implemented in multiple implementations separately or in any suitable subcombination. Moreover, although features may be described above as acting in certain combinations and even initially claimed as such, one or more features from a claimed combination can in some cases be excised from the combination, and the claimed combination may be directed to a subcombination or variation of a subcombination.

Similarly, while operations are depicted in the drawings in a particular order, this should not be understood as requiring that such operations be performed in the particular order shown or in sequential order, or that all illustrated operations be performed, to achieve desirable results. In certain circumstances, multitasking and parallel processing may be advantageous. Moreover, the separation of various system components in the implementations described above should not be understood as requiring such separation in all implementations, and it should be understood that the described program components and systems can generally be integrated together in a single software product or packaged into multiple software products.

A number of implementations have been described. Nevertheless, it will be understood that various modifications may be made without departing from the spirit and scope of the disclosure. For example, example operations, methods, or processes described herein may include more steps or fewer steps than those described. Further, the steps in such example operations, methods, or processes may be performed in different successions than that described or illustrated in the figures. Accordingly, other implementations are within the scope of the following claims.

EXAMPLES

In some implementations, methods for fracturing a formation include: receiving, at a processor, data representing a depth-dependent horizontal in-situ stress for at least two layers of the formation and a fracture toughness for the at least two layers of the formation; predicting, by the processor, a hydraulic fracture profile within the formation for at least one perforation pressure by: determining an upper fracture tip stress intensity factor within the formation for the at least one perforation pressure, the determining being based on the depth-dependent horizontal in-situ stress of one or more layers of the at least two layers of the formation; determining a first vertical depth within the wellbore at which the upper fracture tip stress intensity factor satisfies an upper fracture tip propagation condition based on the fracture toughness of one or more layers of the at least two layers of the formation; determining a lower fracture tip stress intensity factor within the formation for the at least one perforation pressure, the determining being based on the depth-dependent horizontal in-situ stress of one or more layers of the at least two layers of the formation; determining a second vertical depth within the wellbore at which the lower fracture tip stress intensity factor satisfies a lower fracture tip propagation condition based on the fracture toughness of one or more layers of the at least two layers of the formation; determining a height of the hydraulic fracture profile for the at least one perforation pressure, the height extending from the first vertical depth to the second vertical depth; and determining a depth-dependent width of the hydraulic fracture profile at the at least one perforation pressure, the depth-dependent width extending along the height of the hydraulic fracture profile, the predicted hydraulic fracture profile for the at least one perforation pressure being at least partially defined by the determined height and the determined depth-dependent width for the at least one perforation pressure; determining, by the processor, a perforation pressure value based on at least one of the predicted hydraulic fracture profiles of the formation; and injecting, by a pump, a fluid into the wellbore at the determined perforation pressure to hydraulically fracture the formation.

In an example implementation combinable with any other example implementation, the upper fracture tip stress intensity factor, the lower fracture tip stress intensity factor, and the depth-dependent width of the hydraulic fracture profile are determined using Gauss quadrature.

In an example implementation combinable with any other example implementation, the method includes measuring, by a logging device, data representing a Young's modulus, a Poisson's ratio, and a layer thickness for each layer of at least two layers of the formation, wherein predicting the hydraulic fracture profile is based on the Young's modulus, the Poisson's ratio, and the layer thickness for each layer of the at least two layers of the formation.

In an example implementation combinable with any other example implementation, predicting the hydraulic fracture profile of the formation comprises determining, by the processor on a layer-by-layer basis, whether the upper fracture tip stress intensity factor satisfies the upper fracture tip propagation condition and whether the lower fracture tip stress intensity factor satisfies the lower fracture tip propagation condition.

In an example implementation combinable with any other example implementation, predicting the hydraulic fracture profile of the formation comprises dynamically discretizing, by the processor on the layer-by-layer basis, each layer of the at least two layers of the formation into a plurality of one dimensional finite elements having an order greater than one.

In an example implementation combinable with any other example implementation, the depth-dependent width of the hydraulic fracture profile along the height of the hydraulic fracture profile is determined using a plane strain Young's modulus.

In an example implementation combinable with any other example implementation, injecting the fluid into the wellbore causes a hydraulic fracture to propagate until a height of the hydraulic fracture reaches the determined height of the hydraulic fracture profile.

In an example implementation combinable with any other example implementation, the at least two layers of the formation comprises at least five layers of the formation.

In an example implementation combinable with any other example implementation, the upper and lower fracture tip stress intensity factors are determined based on a depth-dependent hydrostatic pressure within the formation.

In an example implementation combinable with any other example implementation, the method includes plotting, by the processor, a fracture height map on a first plot representing a progression of the height of the predicted hydraulic fracture profile as a function of the at least one perforation pressure.

In some implementations, methods for fracturing a formation include: receiving, at a processor, data representing a depth-dependent horizontal in-situ stress and a fracture toughness for one or more layers of the formation; predicting, by the processor, a hydraulic fracture profile of the formation by: determining a location of an upper portion of the hydraulic fracture profile when an upper fracture tip stress intensity factor satisfies an upper fracture tip propagation condition based on the depth-dependent horizontal in-situ stress and the fracture toughness for the one or more layers of the formation, the upper portion being above a perforation location of a wellbore; determining a location of a lower portion of the hydraulic fracture profile when a lower fracture tip stress intensity factor satisfies a lower fracture tip propagation condition based on the depth-dependent horizontal in-situ stress and the fracture toughness for the one or more layers of the formation, the lower portion being below the perforation location of the wellbore; and determining a depth-dependent width of the hydraulic fracture profile based on the determined location of the upper portion and the determined location of the lower portion; and injecting, by a pump, a fluid into the wellbore to hydraulically fracture the formation at the perforation location based on the predicted hydraulic fracture profile.

In an example implementation combinable with any other example implementation, the depth-dependent width of the hydraulic fracture profile is determined using a plane strain Young's modulus.

In an example implementation combinable with any other example implementation, predicting the hydraulic fracture profile of the formation comprises using one-dimensional finite clements to discretize each layer of the one or more layers of the formation.

In an example implementation combinable with any other example implementation, the one-dimensional finite clements are at least second order elements.

In an example implementation combinable with any other example implementation, predicting the hydraulic fracture profile within the formation comprises accounting for a depth-dependent hydrostatic pressure acting on one or more surfaces of the hydraulic fracture profile.

In an example implementation combinable with any other example implementation, the location of the lower portion of the predicted hydraulic fracture profile is in a first layer of the one or more layers, and the location of the upper portion of the predicted hydraulic fracture profile is in a second layer of the one or more layers, the first layer being below the second layer.

In an example implementation combinable with any other example implementation, the upper portion of the predicted hydraulic fracture profile is a upper-most tip of the hydraulic fracture profile, and the lower portion of the predicted hydraulic fracture profile is a lower most-tip of the hydraulic fracture profile.

In an example implementation combinable with any other example implementation, the method includes determining, by the processor, a pump schedule for the pump based on the predicted hydraulic fracture profile, wherein the fluid is injected into the wellbore according to the pump schedule.

In some implementations, methods for fracturing a formation include: receiving, by a processor, data representing (i) fracture toughness properties for a plurality of layers of the formation, and (ii) depth-dependent horizontal in-situ stresses exerted on the plurality of layers of the formation; automatically generating, by the processor, a plurality of one-dimensional elements along a vertical direction of a first layer of the formation; numerically predicting, by the processor, whether a hydraulic fracture is expected to occur within the first layer based on (i) a perforation pressure, (ii) a fracture toughness of the layer based on the received fracture toughness properties, and (iii) a depth-dependent horizontal in-situ stress of the layer based on the received depth-dependent horizontal in-situ stresses; and when the hydraulic fracture is expected to occur in the first layer, injecting, by a pump, a fluid into a wellbore at the perforation pressure to hydraulically fracture the first layer of the formation such that a hydraulic fracture propagates through at least a portion of the first layer of the formation and stops within the first layer of the formation.

In an example implementation combinable with any other example implementation, when the hydraulic fracture is not expected to occur in the first layer, automatically generating, by the processor, a plurality of one-dimensional elements along the vertical direction of a second layer of the formation and numerically predicting whether a hydraulic fracture is expected to occur within the second layer; and when the hydraulic fracture is expected to occur in the second layer, injecting, by the pump, the fluid into the wellbore at a perforation pressure at which the hydraulic fracture is expected to occur in the second layer to hydraulically fracture the second layer of the formation such that a hydraulic fracture propagates through the first layer and stops within the second layer of the formation.

In an example implementation combinable with any other example implementation, when the hydraulic fracture is not expected to occur in the second layer, automatically generating, by the processor, a plurality of one-dimensional elements along the vertical direction of a third layer of the formation and numerically predicting whether a hydraulic fracture is expected to occur within the third layer; and when the hydraulic fracture is expected to occur in the third layer, injecting, by the pump, the fluid into the wellbore at a perforation pressure at which the hydraulic fracture is expected to occur in the third layer to hydraulically fracture the third layer of the formation such that a hydraulic fracture propagates through the first and second layers and stops within the third layer of the formation.