The Study of Options for Identification Stress Contrasts via Pumping History

A recorded pressure history is commonly the only quantitative information available as the result of a hydraulic treatment. It would be of value to use the history for decreasing uncertainty of in situ conditions, first of all, of in situ stresses. Yet the inversion of a pressure history in terms of in situ stresses is an ill-posed problem. This makes unclear the options for identification of the stresses. Of essence is to study when, to which extent and how the identification may be performed. The paper aims to meet this request. We start from formulation the inverse problem with employing truly and pseudo-three-dimensional models for generation experimental and theoretical histories. To simplify identification, we (1) reduce the number of input parameters present in these models by using specially normalized variables, and (2) distinguish three dimensionless key parameters R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R$$\end{document}, kasym\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{\mathrm{asym}}$$\end{document} and bP\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${b}_{\mathrm{P}}$$\end{document}. They define, respectively, the strength of a stress contrast, asymmetry of stress contrasts, and the prevailing type of fracture growth. The analysis, based on these parameters, gives ranges within which the identification is possible in principle. It also yields a sound initial guess to facilitate rapid identification of stress contrasts within the ranges distinguished. Numerical examples illustrate the options for identification of stress contrasts. They demonstrate agreement with the theoretical conclusions. Specifically, on the one hand, they confirm impossibility of identification beyond the ranges distinguished, and, on the other hand, within these ranges, they demonstrate rapid identification of stress contrasts by the method suggested. Study is simplified by using special normalizing and three dimensionless parameters. Ranges of stress contrast identification via pumping history are specified. Sound initial guess is given to facilitate rapid identification of stress contrasts. Numerical examples for pumping histories confirm theoretical conclusions. Method may be used in real time of hydraulic fracturing treatment. Study is simplified by using special normalizing and three dimensionless parameters. Ranges of stress contrast identification via pumping history are specified. Sound initial guess is given to facilitate rapid identification of stress contrasts. Numerical examples for pumping histories confirm theoretical conclusions. Method may be used in real time of hydraulic fracturing treatment.


Introduction
Widely used operation of hydraulic fracturing (HF) is performed in field under quite uncertain data on the properties and state of rock stratum at the depth of a pay-layer. To make a HF treatment a success, it is important to improve the reliability of parameters which influence initiation, propagation and performance of a fracture. Such are in situ stresses in a pay-layer and in situ stresses in its neighbor layers. Thus, for years, great attention has been paid to finding the tensor of in situ stresses or, at least, its minimal principal value [see, e.g. the classical (Haimson and Fairhurst 1969) and recent (Zhang and Yin 2014) reviews]. As regards HF, the most popular and advanced technique used in practice is minifrack treatments to estimate the minimal in situ stress in a pay-layer (e.g., Warpinski and Smith 1989). A mini-frack test refers to the initial stage of a HF propagation, and its results are interpreted in terms of the minimal in situ stress in the pay-layer.
The knowledge of this stress is vital for a proper choice of pumping pressure and pumping equipment. Yet, it is insufficient to conclude on a proper pumping regime, in particular, on a pumping rate. This impels complementing the in situ stress in the pay-layer with estimations of the in situ stresses in its two closest neighbors. The differences between these stresses and the stress in the pay-layer are stress contrasts. As follows from general considerations (e.g., Geertsma 1989;Mack and Warpinski 2000), they drastically influence HF propagation, its final sizes, opening and, consequently, the success of a treatment.
This stimulated developments and successive improvements of the pseudo-three-dimensional (P3D) model (e.g., Clifton and Abou-Sayed 1979;Clifton 1989;Morales 1989;Mack and Warpinski 2000;Settari and Cleary 1986;Rahman and Rahman 2010;Dontsov and Peirce 2015;Linkov and Markov 2020) of HF. As the output, the model provides the net-pressure and opening distributions, net-pressure history at the inlet, average openings and fracture footprints as functions of time. Thus, it connects a net-pressure history with stress contrasts. This suggests its employing for identification of stress contrasts by comparing an observed pressure history with histories calculated for arbitrary in situ stresses. However, as usual for inversion problems and illustrated in the pioneering paper by Mack and Warpinski (2000), an identification process may diverge or give quite different solutions for very small changes of input parameters. In particular, when using a P3D model, close net-pressure histories may correspond to quite different stress contrasts. Then the inversion would lead to wrong practical decisions.
Still, an observed pressure history is an important source of information on a HF propagation available in the field (e.g., Nolte and Smith 1981;Nolte 1989). Therefore, of essence is to study when, to which extent and how the identification of stress contrasts via a pressure history may be performed. The paper aims at this goal. For reaching the goal, we (i) reduce the number of input parameters by special normalizing variables; and (ii) distinguish three key parameters R , k asym and b P , which characterize, respectively, the strength of a stress contrast, asymmetry of stress contrasts, and the prevailing type of fracture growth (channelized or contrast-free). The analysis, employing their limiting values, gives the ranges of possible identification and a proper initial guess for rapid identification of stress contrasts within the ranges.

Inverse Problem of Identification Stress Contrasts
Let p ex (t) be a given experimental pressure history, observed in the field or laboratory for a pay-layer between neighbor layers (Fig. 1). Since the history was influenced by in situ stresses in these layers, it would be of value to extract these stresses from the curve p ex (t) . This is a classical inverse problem (e.g., Aster et al. 2018), solved, when possible, numerically by employing a finite number N of experimental values. When this number is greater than the number m of identified parameters ( N > m ), the inversion is performed by the standard least-square approach. In our case, N is the number of instances t i with sampled values P i ex = p ex (t i ) of a pressure history. They are assembled as a vector-column P ex with N components. The number N is taken large enough to reproduce a pressure history to an acceptable accuracy. Normally, it is sufficient to take N > 10 ; the number N = 100 is well above this value. The number m of identified in situ stresses is notably less; in the problem considered, it is two or three. Thus, we always have N > m , which suggests using the least-square method. By this method, we need to have a mathematical model, which generates theoretical pressure histories p th (t) for a wide range of stress contrasts. A theoretical history, when discretized at the same time instances t i ( i = 1, … , N ), yields a vector of theoretical pressure history P th with components P i th = p th (t i ) . The vector of residuals r = P ex − P th gives the misfit between the observed and theoretical histories. The identified parameters are found to minimize Fig. 1 Scheme of a planar fracture the Euclidean norm of the misfit r . In general, their best-fit values are searched by using a specially constructed procedure. If a success, the found best-fit parameters yield the best-fit theoretical pressure history P bf . The residual vector r bf = P ex − P bf , corresponding to P bf , reveals if the theoretical model can be considered as realistic or not. The search of the best-fit values is commonly performed by using the leastsquares method in Newton-Raphson iterations, employed to account for non-linearity. Its application to the inversion of a pressure history is outlined in Appendix 1. The appendix contains also comments on its numerical implementation for identification of stress contrasts.
True difficulties arise when starting practical calculations. They are caused by the basic mathematical fact: the inverse problem is ill-posed. For this reason, the iterations often diverge or give quite different solutions for very small changes of input parameters.
Thus there are two major problems to be solved. The first is to find the ranges, within which the inversion in terms of stress contrasts is possible in principle. The second is to suggest a sound initial guess, which would facilitate convergence of iterations within these ranges. The discussion in the following sections is focused on solving these two problems.
Clearly, before starting the general study, we need to specify the experimental p ex (t) and theoretical p th (t) pressure histories, which are basic for an identification of stress contrasts. The next two subsections contain these prerequisites.

Model for Generation Experimental Pressure History
Pressure histories, observed in field, are strongly influenced by a number of factors other than in situ stresses.
Hence, their numerical values should be present in a theoretical model, serving for generation the theoretical history p th (t) . Unfortunately, some of these factors, similar to in situ stresses, are also unavailable or quite uncertain in field conditions. Thus, we need either to exclude these factors from the model, or to add them to the identified parameters. The first option makes the model inadequate. The second makes the inversion problem practically unsolvable. As shown by the analysis presented in following sections, even when knowing exact values of the input parameters additional to in situ stresses, and even when employing the theory, facilitating identification, the inversion is impossible in some ranges of input parameters. It would become absolutely impossible if trying to identify all the parameters unknown in field conditions. This excludes using field pressure histories as the source of assigning experimental pressure histories p ex (t) for the general study aimed. The reasonable compromise for assigning experimental pressure histories p ex (t) consists of using laboratory experiments. For them, all the input data, needed to generate experimental histories, are known in advance. However, the number of available experiments is quite limited because of their complexity and cost, which drastically reduces options for the general analysis. The options may be extended by employing a truly 3D model, which accurately reproduces the experimental results available.
Consider a quite comprehensive scheme, used for modeling truly 3D planar hydraulic fractures. Its exposition includes brief comments on the physical essence of underlying equations and on input parameters significant for pressure histories. The equations, being well-known, are not reproduced explicitly; they may be found elsewhere (e.g., Adachi et al. 2007;Peirce and Detournay 2008;Peirce 2015;Linkov 2015;Chen et al. 2020).
Geometry: A planar fracture is driven by a fluid, pumped since the initial time t = 0 at the inlet in the middle plane of a horizontal pay-layer (Fig. 1). The inlet is associated with a source point. It is taken as the origin O of the Cartesian coordinates x, y, z . Thus the source is assumed pointed. The thickness of the pay-layer is H ; the neighbor layers are half-spaces.
In situ stresses, stress contrasts and net-pressure: The fracture propagates in the plane xz orthogonal to the direction y of the minimal in situ stress (z) = yy (z) . In the paylayer ( −H∕2 ≤ z ≤ H∕2 ), the in situ stress is P = yy (0) . It is found in special experiments (e.g., Haimson and Fairhurst 1969). The in situ stress l in the lower ( z < −H∕2 ) and u in the upper ( z > H∕2 ) layer may be arbitrary. For convenience, compressive stresses are assumed positive.
The stress contrasts Δ l and Δ u are defined as the differences between the in situ stresses in neighbor layers and that in the pay-layer. Similarly, the fluid net-pressure p net is defined as the difference between the pumping pressure p and the in situ stress P : When an in situ stress in a neighbor layer ( |z| > H∕2 ) is greater than that in the pay-layer, the stress contrast is positive ( Δ > 0 ), and the corresponding boundary is a barrier for fracture propagation. In the opposite case ( Δ < 0 ), the boundary is an accelerator of fracture propagation.
Fluid properties and equations for flow in narrow channel: A fracturing fluid is assumed to have the power-law rheology with the behavior index n and the consistency index M (see, e.g., Montgomery 2013). For a flow in a narrow channel with the width of the order of the fracture opening w , the fluid is assumed incompressible. Its flow in such a channel is described by quantities (opening w , net-pressure p net , particle velocity v, flux q = wv ), averaged over the width of the channel. The continuity and movement equations are formulated through these quantities as the well-known scalar equation of volume conservation and the 2D Poiseuille-type 1 3 movement equation (e.g., Adachi et al. 2007;Peirce and Detournay 2008;Peirce 2015). For the power-law fluid, the single dimensional characteristic of fluid rheology, entering fluid equations, is � = M , where = 2 2(2n + 1)∕n n is dimensionless.
Rock properties and dependence of net-pressure on opening: The rock mass is assumed isotropic linear-elastic medium with the Young's modulus E and the Poisson's ratio . In general, these constants may be different in pay-and neighbor layers. Then the connection between the net-pressure p net and the opening w involves three boundary integral equations (BIE) in 3D vectors of tractions and displacement discontinuities at the fracture surface. For a planar fracture, the connection between p net and w is reduced to a single elasticity BIE by suggestion that the difference between elastic properties of the layers may be neglected as compared with the influence of differences between in situ stresses in the layers (e.g., Peirce 2015; Zia and Lecampion 2019; Chen et al. 2020). The constants E and enter the BIE joined as the plane-strain elasticity modulus E � = E∕(1 − 2 ).
Leak-off losses: In general, rock mass is permeable, and a fracturing fluid leaks to pores and natural cracks. Commonly, the leak-off losses are accounted by means of the Carter's leak-off coefficient C leak (Carter 1957). They may strongly influence the fracture propagation (e.g., Howard and Fast 1970;Nordgren 1972;Lenoach 1995;Madiarova 2003;Mack and Warpinski 2000;Adachi et al. 2007). When accounting for leak-off losses in numerical modeling, the continuity (volume conservation) equation includes the term with the factor C leak .
Fracture (rupture) condition: A fracture is initiated and propagates because the strength limit of a real material is finite. Above the limit, a material is ruptured. In continuum description, this means arising of a surface element with displacement discontinuity, that is a fracture. Thus to initiate and maintain fracture propagation, a fracture condition has to be met. For hydraulic fractures, this condition presumes the tensile mode of rupture. It is taken in the simplest form of linear fracture mechanics (e.g., Rice 1968) as K I ≥ K IC . Herein, K I is the normal stress intensity factor, defined by external loading, K IC is its critical value (toughness). The critical SIF K IC is a physical property of a material, characterizing its strength limit. When modeling HF propagation, it is always assumed, while sometimes implicitly, that the strength condition is fulfilled. Consequently, the critical SIF K IC is always assigned as an input parameter, while for some propagation regimes, it is implicitly supposed that K IC = 0.
Speed equation: When a fracture condition is met, a fracture grows with a propagation speed v * of its contour (front). The growth may occur in various regimes, depending on the relative input of factors resisting or favorable to front advance. Such are limiting regimes of dominated toughness, viscosity, or leak-off, and regimes intermediate between these extremes (e.g., Rice 1968;Spence and Sharp 1985;Desroches et al. 1994;Garagash et al. 2011;Linkov 2015). In any case, tracing a surface propagation always involves a speed equation (e.g., Sethian 1999). For a hydraulic fracture, the speed v * equals to the normal component of the fluid particle velocity, when a particle tends to a front point (e.g., Kemp 1990;Linkov 2015). This implies that the propagation speed is entirely defined by the leading asymptotics of the opening at near-front zone. Consequently, a speed equation does not contain input parameters additional to those already listed.
Pumping rate: After fracture initiation, its surface grows to sizes notably greater than the borehole diameter. This justifies modeling the fluid influx from the borehole by a pointed source of a given intensity Q(t) . As mentioned, for certainty, it is located at the origin. A pumping regime is assigned as Q(t) = Q 0 f 0 (t) where Q 0 is a typical pumping rate Q 0 , and f 0 (t) is a dimensionless function of time. If the pumping rate is constant, f 0 (t) = 1 . The time t is counted from the pumping start. Finally, the fluid influx enters the continuity equation as the term Q 0 f 0 (t) (x) (z) , where (x) and (z) are Dirac's delta-functions. Emphasize that the only dimensional parameter of a pumping regime is Q 0 .
Boundary and initial conditions: Since the source is pointed and pumping starts at t = 0 , the initial condition is w = 0 in the entire plane of perspective propagation. For any time instant t > 0 , the only boundary condition is: w = 0 at points of a current fracture contour (the conditions at the interfaces z = ±H∕2 are accounted via stress contrasts).
Boundary value problem and its solution: The equations and conditions listed above define a truly 3D comprehensive model of hydraulic fracturing. Its five equations (volume conservation, two components of the 2D movement equation, elasticity, and speed equation) comprise the system of five equations in five unknown quantities: the opening w ; two components of the particle velocity v (or the flux q = wv ); the fluid pressure p (or the corresponding net-pressure p net = p − P ); and the propagation speed v * . The equations are complemented with the rupture, initial and boundary conditions. The system includes the ten input parameters defined above. These are: the pay-layer thickness H ; in situ stresses P , l , u (and respective stress contrasts Δ l , Δ u ); fluid behavior index n and generalized consistency index ′ ; plane-strain modulus E ′ ; Carter's leak-off coefficient C leak ; fracture toughness K IC ; and typical pumping rate Q 0 .
This comprises the boundary value problem (BVP) for the truly 3D model of hydraulic fracturing. Its solution includes dependences of the pressure (accordingly, netpressure) on time at points of the fracture surface. Presently, the BVP may be solved by the contemporary methods of implicit or explicit temporal integration of an algebraic system, obtained after spatial discretization of its equations (e.g., Peirce and Detournay 2008;Peirce 2015;Zia and Lecampion 2019;Linkov 2019;Chen et al. 2020). However, a numerical solution obtained by these methods is quite time-costly. Therefore, while the truly 3D solution is appropriate for generation individual pressure histories, needed for the general study, it is inappropriate for multiple runs, performed in iterations with changed input data. Adequacy of the model; generation experimental histories: Available comparisons of solutions, obtained by using the truly 3D model, with results of laboratory physical experiments show good agreement (Dontsov and Peirce 2015;Zia and Lecampion 2019). This suggests using the truly 3D model for generation pressure histories, observed in laboratory experiments. For these histories, all the input parameters are known in advance. This makes them applicable when studying options for inversion of a pressure history in terms of in situ stresses. Thus in the following study, an experimental pressure history p ex (t) is associated with that obtained for the truly 3D model with fixed values of the input parameters listed. The input parameters include stress contrasts Δ l , Δ u to be identified via p ex (t).
Numerical implementation: We have developed a computer code for the truly 3D model. Specifically, following Chen et al. (2020), the solutions to the truly 3D problem are found by using the explicit stabilized Runge-Kutta-Legendre (RKL) method of second order (RKL2) (Meyer et al. 2014) for integration of ordinary differential equations (ODE), resulting spatial discretization. In difference with Chen et al. (2020), the system of ODE includes speed equations in accordance with Linkov (2019). The examples, studied in Chen et al. (2020), served us to check the accuracy of our results; they coincide to the accuracy of graphical representations used in the cited paper.

Model for Generation Theoretical Pressure Histories
As mentioned, calculations by using the truly 3D model are time-costly. Thus, for studies, involving multiple runs with changing input parameters, it is reasonable to follow the line by Mack and Warpinski (2000) and to employ a P3D model. This model facilitates drastic reduction of time expenses without significant losses in the accuracy of calculated quantities for fractures lengthened along the x-axis.
Below it is used in the form of the improved (IP3D) model (Linkov and Markov 2020), which includes and carefully takes into account all the input parameters listed for the truly 3D model. Numerical implementation: The special code, developed for this study to solve equations of the IP3D model, employs explicit integration by the simplest forward Euler method. There is no need to employ the stabilized RKL2 method, because for the IP3D model, the CFL stability condition is not too restrictive.
Ranges of applicability IP3D model: The theoretical model has to be realistic. Thus, when using the IP3D model, it is necessary to distinctly define the ranges of its applicability to the identification problem considered.
This model, as any P3D model, is based on the main assumption that the fracture size along the x-axis exceeds that along the z-axis, so that cross sections are in the plane-strain state at some distance for the front. Recall that the fluid source is located at the center of a pay-layer with constant in situ stress. Hence, till the time t H , at which a fracture reaches the boundaries of a pay-layer, the propagation is axisymmetric; the fracture is circular. Clearly, the main assumption is not met for a circular footprint. Hence, a P3D model becomes applicable only for time exceeding the time t H . This implies that an inversion by means of the IP3D model should involve the part of a pressure history well-beyond the initial interval [ 0, t H ]. Normally, t H is much less than the final time t F of a hydraulic treatment: t H is of order of minutes, while t F has the order of hundreds minutes. Thus the major part of a pressure history may be available for using the IP3D model. The question is: for which stress contrasts, the model adequately describes the pressure history on this part? Some preliminary conclusions on the applicability of the IP3D model follow from the publications Dontsov and Peirce (2015) and Linkov and Markov (2020), containing the comparison of fracture footprints for the truly 3D and IP3D models. It is established that for symmetric stress barriers, the solutions become in good agreement when the fracture half-length exceeds the pay-layer thickness. This implies, that the IP3D model is applicable for solving the inverse problem, when (i) the stress contrasts are positive and symmetric, and (ii) the time exceeds first ten minutes. In the limiting case of practically impenetrable contrasts, corresponding to the Nordgren (1972) problem, numerical results, summarized in Appendix 2, show very high accuracy provided by the code for IP3D model.
Clearly, the main assumption is violated in cases when stress contrasts are zero, negative or strongly non-symmetric. Then any P3D model becomes inapplicable for accurate quantitative studies. Despite of the obvious deficiency, the IP3D model still may provide illuminating qualitative conclusions and rough estimations regarding to footprints and stress contrasts. This is due to strong influence of the mass (volume) conservation law on the solution. Numerical realizations of both truly 3D and IP3D models meet this law identically at each time instant and any time interval. This important feature suffices qualitative and rough quantitative similarity of the results. A detailed study, also summarized in Appendix 2, shows that even in unfavorable limiting cases, when one or both stress contrasts are zero, the IP3D model, besides providing illuminating qualitative conclusions, may serve for rough quantitative estimations regarding footprints and stress contrasts.
Pseudo-experimental histories: Calculations by employing the codes for truly and IP3D models provide a priory quantitative information on the accuracy to which the theoretical pressure history, obtained by means of the IP3D model, agrees with the experimental history, obtained by means of the 3D model, under the same input data. Indeed, for a set of used time instances t i ( i = 1, … , N ), and fixed input parameters (including those identified), we have graphical representations of both the experimental P ex and pseudoexperimental P pex histories. Thus, the corresponding misfit r = P ex − P pex is known in advance. In numerical examples of Sect. 7, this beneficial option is offered by plotting the curves of experimental and pseudo-experimental pressures. Their comparison gives one-glance knowledge on a time interval (if any), in which the IP3D model is applicable for identification of stress contrasts for a particular set of input parameters. The Euclidean norm of r gives a quantitative measure of the misfit.
Comment 1: As agreed, an experimental net-pressure history p ex (t) is found for the truly 3D model, in which the fluid supply is modeled by the Dirac's delta-function. Using the latter simplifies the study of fluid injection into a borehole of a small radius r bh ( r bh ≪ H) . However, in cases when rocks, embedding a pay-layer of height H , are practically impenetrable for fracture propagation, a pointed source generates a boundary layer effect. A perturbation occurs in an area of the characteristic size H near the source. An asymptotic analysis of the truly 3D solution for this extremal case has shown that the net-pressure, calculated at the cross section of the source point, is unavoidably less than that for the corresponding P3D (Nordgren's) case. Thus, it may be expected that, in cases of strong barriers, the pseudo-experimental net-pressure histories p pex (t) will systematically give greater values of the net-pressure than corresponding experimental histories p ex (t) . The calculations, performed by the authors and discussed below, clearly reveal this tendency in cases of strong contrasts, when the fracture propagation is notably channelized.

Extension to Multilayered Formation
The problem formulation above has been presented for the classical three-layer scheme with the pay-layer between half-spaces. The formulation may be promptly extended to a multilayer formation with arbitrary stress contrasts. Indeed, the formulation employs the assumption that the difference between elastic properties of the layers may be neglected as compared with the influence of differences between in situ stresses. Thus the medium itself is considered homogeneous. This assumption notably simplifies computations by avoiding the need in building, storing and repeatedly using a special Green function for layered medium. Then implicit, semi-explicit and explicit algorithms, mentioned in Sect. 2.2 and serving for generation of an experimental pressure history, stay applicable for a package of an arbitrary number of layers between embedding half-spaces. The layers in the package include the pay-layer and may have arbitrary stress contrasts. Examples given in the papers Peirce (2015Peirce ( , 2016, Zia and Lecampion (2019) and Chen et al. (2020), illustrate this option.
The IP3D model, discussed in Sect. 2.3 and serving for rapid generation of theoretical net-pressure histories, is also of immediate use for any package of layers with various stress contrasts. An example, given in Linkov and Markov (2020), illustrates this option.
Specifically, the both codes (for truly 3D and IP3D models), developed by the authors for generation of experimental and theoretical pressure histories, are applicable for any number of layers having the same elastic modules and arbitrary stress contrasts. Thus, formally any number of stress contrasts may be included into the set of identified parameters. The identification, for instance, may be performed by the standard least-square method, used in this paper. However, as mentioned, the true problems arise when starting practical calculations.
Inclusion of additional stress contrasts into the set of identified parameters, as any extension of the set, aggravates the consequences of the problem ill-posedness. As clear from the further discussion and examples, this commonly leads to divergence of iterations if not having a sound initial guess. Thus it is reasonable to start from the simplest threelayer scheme, for which the ranges of convergence and the choice of a proper initial guess are studied in the following sections. Then we at least obtain either a clear indication of impossibility of stress contrast identification, or the best-fit estimation of average contrasts in the closest layers. In the first case, extending the set of identified parameters looks senseless. In the second (favorable) case, the found best-fit values may serve as an appropriate initial guess for solving the problem with extended number of identified parameters, in particular, with additional stress contrasts in a package of, say, five layers.

Input Parameters. Reducing Their Number
The input parameters, listed above, are the same for the both models. In total, their number is ten. The number of their various numerical combinations, which are to be studied, remains very large even when some parameters are fixed by their typical values. This notably complicates the study of the inversion problem. To simplify the study and to obtain results in a general form, of value is to reduce the number of input parameters to minimum. The reduction is achieved by a proper normalizing physical quantities. The analysis below employs the normalized variables, using of which excludes three of the ten listed parameters from formulations of both the truly 3D and IP3D models. The appropriate variables are (Linkov 2016b(Linkov , 2019: with the normalizing coefficients This normalizing excludes the parameters ′ , E ′ , and, what is of essence, Q 0 . Of importance is also that, although the normalized variables have unusual dimensions, the normalizing does not change the time units. This is especially convenient when considering functions of time, in particular, pressure histories and footprints. In the normalized variables, the equations of the truly 3D and IP3D models become those, considered above in physical variables, when setting � = 1, E � = 1 , Q 0 = 1 , and replacing the physical quantities with their normalized counterparts. Further reduction of input parameters is achieved by neglecting the contribution of the toughness K IC into the resistance to fracture propagation, as compared with the resistance caused by the viscosity ′ . Indeed, as shown in Savitski and Detournay (2002), the viscous resistance dominates over toughness during the major time of a fracture propagation. Thus normally it is assumed (sometimes implicitly) that K IC = 0 . Consequently, K ′′ IC = 0 , as well. This assumption excludes the need to perform calculations for various values of the toughness K IC .
The most of the following discussion refers equally to Newtonian and non-Newtonian (commonly thinning) fluids. For certainty and having in mind comparisons with published results for planar fractures, we shall focus on the case of a Newtonian fluid ( n = 1 ). Still, the cases of perfectly plastic ( n = 0 ) and thinning ( 0 < n < 1 ) fluids are sometimes included into discussion, as well. Further simplification follows from the fact that both the truly 3D and IP3D models do not contain the in situ stress P in the pay-layer. Merely the stress contrasts Δ l , Δ u and the net-pressure p net enter the planar models. Thus their solutions do not depend on the reference value P of the in situ stress; the latter may be arbitrary. Specifically, it can be found in field measurements by methods reviewed in Haimson and Fairhurst (1969). Its value becomes of essence on the next stage of an analysis, when the found net-pressure is summed with P to find the actual pressure of a pumped fluid. Thus, when studying the inversion problem in terms of stress contrasts, the in situ stress P drops out from consideration. However, since normally P ≫ p net , even small errors in an estimation of the in situ stress P may strongly influence the applicability of found stress contrasts to field conditions. In Sect. 6.2, we shall specially discuss this issue.
Finally, for a constant pumping rate, the number N in of input parameters, which influence a solution, reduces to four only ( N in = 4 ). These are: the normalized height H ′′ , the normalized stress contrasts Δ l�� , Δ u�� , and the normalized Carter's leak-off coefficient C �� leak . Some of them, say Δ l�� , Δ u�� , are to be identified. If wanted, the normalized in situ stress �� P may be added to known or identified parameters; then N in = 5.
Comment 2: Obviously, the two constants ′ , E ′ may be promptly excluded by introducing the normalized time t �� = t∕t n , where t n = � ∕E � 1∕n is the intrinsic time of the system. However, this normalizing does not exclude the pumping rate Q 0 . Besides, for usual values of ′ and E ′ , the intrinsic time t n is less than 10 −7 s, which causes extremely inconvenient scale of the normalized time t ′′ .

Propagation Through Boundary of Pay-Layer: Parameter of Contrast Intensity
Recall that a hydraulic fracture propagates only when a rupture condition is met at points of its front. At the interface of two layers, the properties, asymptotics of fields and propagation conditions abruptly change. This strongly influences further propagation of a crack. In frames of linear fracture mechanics, when a cohesion zone is small as compared with linear sizes of a crack, the studies are always performed asymptotically (e.g., Cook and Erdogan 1972;He and Hutchinson 1989). It is assumed that a crack is semiinfinite, and a vicinity of its tip is in the plane-strain state.
In cases of dissimilar elastic materials, the penetration or deflection conditions are essentially defined by the dimensionless parameter, which is the ratio of plane-strain modules of the two materials in contact. In our problem, when neglecting the difference in the modules, the only quantities, which determine the asymptotic fields around the tip, are the average net-pressure p av in the crack and the stress contrast Δ in the half-plane ahead of its tip (Fig. 2). Then the only non-dimensional parameter, characterizing hydrofracture penetration, is the ratio This parameter defines the actual intensity of a stress contrast in reference to the net-pressure activating the penetration. Clearly, when a stress contrast is zero ( Δ = Δ �� = 0 ) or infinite ( Δ = Δ �� = ±∞ ), the corresponding intensity is zero ( R = 0 ) or infinite ( R = ±∞ ), as well.
For a fracture, initiated and driven by a pointed source, located at the center of a pay-layer, the propagation is axisymmetric till the entire crack is within the pay-layer. Thus to apply the parameter R to estimation of contrast intensity in the problem considered, we need to estimate the average net-pressure p av at the instant t H , when the fracture front reaches the boundaries z = ±H∕2 . Both t H and p av (t H ) may be found from the self-similar solution to the axisymmetric problem. The technical details of the derivation are given in Appendix 3.
For the time of reaching the boundaries, the solution gives where * n has the meaning of the self-similar fracture radius; its values for thinning fluids are tabulated in Linkov (2016a). As mentioned above, for typical ranges of the input parameters, the time t H does not exceed minutes. Normally leakoff losses may be neglected for an initial interval of such duration (see, e.g., Zia and Lecampion 2019). This justifies using Eq. (5), derived for impermeable rock, in cases of permeable rock mass. For the parameter R , the derivation outlined in Appendix 3, gives where k R = (1∕P av ) 0.5∕ * n 1.5n∕(n+1) . The constant P av = 2∫ 1 0 P( ) d has the meaning of the average selfsimilar net-pressure. The distributions P( ) over the fracture radius for thinning fluids are given in Linkov (2016a). For a Newtonian fluid, P av = 0.479 , * n = 0.698 , and k R = 1.63 . As shown in Appendix 3, for power-law fluids, k R ≈ (2∕3)(1.43) 1.5(n+2∕(n+1) . For a Newtonian fluid ( n = 1 ), the approximate value is k R ≈ 1.49 ; it differs merely 8.6% from the given above exact value k R = 1.63.
With factor k R known, Eq. (6) entirely defines the intensities of the stress contrasts at the upper ( Δ �� = Δ u�� ) and lower ( Δ �� = Δ l�� ) boundaries of the pay-layer in terms of the normalized input parameters H ′′ and the fluid behavior index n . From the definitions (2), (3), it follows that in terms of physical quantities, Eq. (6) reads The right hand side includes two dimensionless ratios: Δ ∕E � and H 3 ∕(Q 0 t n ) . The first of them presents a rough estimation of the strain, caused by a stress contrast Δ in an elastic medium with the Young's modulus E . The second is t Q ∕t n , where t Q = H 3 ∕Q 0 is the time needed to fill a cube of the volume H 3 by a fluid pumped with the pumping rate Q 0 . This implies that in HF problems, the intensity R of a stress contrast is proportional to the n∕(2n + 2)-th degree of the ratio H 3 ∕Q 0 . For a Newtonian fluid, the degree is 1/4; then R is proportional to H 3∕4 . From (7), it is obvious that a positive (negative) value Δ of the stress contrast does not define the intensity of a barrier (accelerator) alone. Much depends also on the pumping rate Q 0 , the pay-layer height H and the intrinsic time t n . A barrier (accelerator) may be strong, weak or intermediate, depending on these parameters. In particular, with decreasing pumping rate, the intensity R goes to infinity, and whatever small barrier becomes impenetrable. In the opposite . Fig. 2 Asymptotic scheme explaining the stress contrast intensity case, with increasing pumping rate, whatever large barrier becomes "invisible" for fracture propagation. Quantitative estimations of strong (impenetrable) and weak ("invisible") contrasts are given in Sect. 6.1. Clearly, when having the intensity R defined, the corresponding stress contrast Δ becomes known, as well, Thus, identification of stress contrasts Δ u , Δ l (or their normalized values Δ u�� , Δ l�� ) is actually equivalent to identification of the intensities R l and R u at the lower and upper boundary, respectively.

Parameter of Contrasts Asymmetry. Reducing Analysis to Cases of Symmetric and Utmost Asymmetric Contrasts
The identification by means of the IP3D model is possible only, when both stress contrasts are non-negative. Then, for the scheme of a pay-layer between half-spaces ( Fig. 1), without loss of generality, we may assume that the non-negative stress contrast Δ l at the lower half-space is greater or equal to the non-negative stress contrast Δ u at the upper halfspace: Δ l ≥ Δ u ≥ 0 . In terms of the intensities, this reads Introduce the contrasts asymmetry parameter as the ratio (Linkov and Markov 2020): From now on, we are focusing on the practically significant cases of non-negative (non-accelerating) contrasts ( R l , R u ≥ 0 ). Then with changing the upper stress contrast Δ u (and consequently R u ) from its maximal value Δ u = Δ l (respectively, R u = R l ) to its minimal value Δ u = 0 (respectively, R u = 0 ), the asymmetry parameter (9) continuously changes from 0 to 1: From (10), by continuity, it appears that we need to study in detail merely two limiting cases of asymmetry, which correspond to the bounds of the interval. These are (i) symmetric contrasts: R u = R l ≥ 0 ; k asym = 0 , and (ii) utmost asymmetric contrasts: R l = ∞ , R u ≥ 0 ; k asym = 1 By continuity, for intermediate cases, the results are easily obtained by interpolations between these two. Thus considering the two limiting cases significantly simplifies theoretical and numerical studies. For each of the limiting cases, it is sufficient to change merely one value ( R u ) of the intensities from zero to infinity. In the first case, we obtain results for equal stress contrasts R l = R u ( k asym = 0) . In the second, the results correspond to very large values of R l ( R l ≫ R u , k asym ≈ 1 ); the exact meaning of "very large" values is established numerically as those, for which the numerical solution practically does not change with further growth of R l .

Parameter of Pressure History for Channelized and Contrast-Free Propagation
Consider three limiting cases, which do not include finite nonzero values of stress contrasts. The first two of them refer to symmetric, while the third refers to utmost asymmetric contrasts. The cases are: (CH) channelized propagation, when the both stress contrasts are infinite ( Δ l = Δ u = +∞), (CF) contrast-free propagation, when the both stress contrasts are zero ( Δ l = Δ u = 0 ), and (SCF) semi-contrast-free propagation, when one of the stress contrasts, for certainty lower, is infinite ( Δ l = +∞) , while the other (upper) is zero ( Δ u = 0).
In the case (CH), the fracture, after an initial interval of reaching the boundaries, propagates entirely along the pay-layer in accordance with the classical PKN (Perkins and Kern 1961;Nordgren 1972) model. The propagation is channelized in the pay-layer. Underline that this case includes both symmetric and utmost asymmetric propagation. Indeed, while both Δ l and Δ u go to plus infinity, their ratio Δ u ∕Δ l may be arbitrary in the interval (0, 1] . Consequently, by the definition (9), the asymmetry parameter k asym is zero, when Δ u ∕Δ l = 1 , and it may be whatever close to 1, when Δ u ∕Δ l ≪ 1 . Thus the study of this case provides the same asymptotics of the solution for both symmetric and utmost asymmetric contrasts.
In the case (CF), the problem is axisymmetric. It corresponds to symmetric contrasts ( k asym = 0).
In the case (SCF), the problem corresponds to the utmost asymmetric propagation ( k asym = 1 ), when one of the contacts is contrast-free. Remarkably, in this case, a solution to the hydrofracture problem is reduced to the corresponding solution for the case (CF) with appropriate interpretation of input and output values.
Indeed, the solution to the axisymmetric problem with the pointed source of intensity Q at the fracture center may be considered as the limit of the solution for the case when two equal pointed sources (of intensities Q∕2 ) with a small distance d between them are located symmetrically about a plane orthogonal to the fracture and passing through its center. At distances r from the center, the solutions will practically coincide when r ≫ d . The plane orthogonal to the fracture is the plane of symmetry for the problem with two sources. Hence the normal component 1 3 of the flux on it is zero. Then the solution corresponds also to a half-space with the source of intensity Q∕2 and with zero influx through its plane boundary. The latter is actually impenetrable for fluid. In terms of stress barriers, the impenetrable boundary may be attributed to the stress contrast of infinite intensity R.
This implies that the solution to the HF problem for a half-space with a pointed source of a given intensity Q 0 located at the distance d∕2 = H∕2 from the impenetrable ( R = ∞ ) plane boundary is approximately the same as the solution to the benchmark axisymmetric problem with the pointed source of doubled intensity 2Q 0 . With growing distance r from the source in the fracture plane, which is the case when employing the IP3D model, the difference between the solutions goes to zero. Thus with growing time, the problem for the semi-contrast-free (SCF) case becomes equivalent to the axisymmetric problem with doubled pumping rate. When using this implication in the normalized variables (2), (3), a "symmetric" solution for the fracture radius remains applicable with the additional factor 3 √ 2 in front of radius. The dependence of the net-pressure on time does not change at all. The conclusions on the net-pressure histories, obtained for the case (CF), equally refer to the case (SCF).
Thus it is sufficient to focus on the cases (CH) and (CF). For obtaining general conclusions, it is also reasonable to study these cases for the extreme values of the leak-off parameter C leak : (A) impermeable rock, when fluid losses are negligible ( C leak = 0), (B) highly permeable rock, when fluid losses are large ( C leak → ∞).
A solution for impermeable rock is obtained by neglecting the leak-off term in the mass conservation (continuity) equation. In its turn, a solution for highly permeable rock is obtained by neglecting in this equation the storage term, which actually accounts for the change of the fracture opening in time. Non-trivial special analysis is required to conclude if it is possible to neglect the influence of the loss/storage term as compared with the influence of the storage/loss term. In many practically significant cases, the both terms are to be kept. Much depends on particular input parameters and on the time of fracture propagation, which is of main interest for a problem considered. The discussion of the non-simple question when the fracture propagation will occur in the storage dominated regime, or in the leak-off dominated regime, or in an intermediate regime, is well apart from our theme. This specific issue may be clarified by employing the variables and analysis of the paper by Nordgren (1972). In numerical calculations, discussed below, the case of large fluid losses was characterized by those values of the normalized coefficient C �� leak , for which the numerical solution practically coincided with the asymptotic analytical solution, given in Appendix 4.
Thus we assume that the values of C leak , corresponding to negligible and large losses, are specified. Then known analytical solutions are available for the cases of small ((A-CH), (A-CF)) and large ((B-CH), (B-CF)) fluid losses. They are promptly expressed in terms of the normalized variables (2). Their general feature is that the solutions are given by functions monomial in time. Consider these cases.
(A-CH) Impermeable rock: channelized propagation: For impermeable rock, the extreme case of channelized propagation ( Δ l = Δ u = +∞ ) corresponds to the classical PKN model, whose exact mathematical formulation for a Newtonian fluid is found by Nordgren (1972). For an arbitrary power-law fluid, its self-similar analytical solution (Linkov 2013) gives for the normalized net-pressure history: where b p and W 0n depend merely on the behavior index n . Specifically, b p = 1∕(2n + 3) ; for a Newtonian fluid ( n = 1 ), W 0n = 1.326 ; for a perfectly plastic fluid ( n = 0 ), W 0n = 1.442 ; in intermediate cases, W 0n = 1.38 to the accuracy of 4.5%.
As expected, the dependence (11) of the net-pressure is monomial in the time. Of essence is that the dimensionless exponent b p does not depend on the normalized thickness H ′′ . This remarkable property of the normalized solution is confirmed by detailed calculations, discussed below. It drastically simplifies the general analysis of the inversion problem. In particular, it serves to obtain an appropriate initial guess to guarantee convergence of the identification process and to make it rapid.
Pay attention, that for the considered case of channelized propagation, the exponent b p of the net-pressure history is positive, being equal +1∕5 for a Newtonian fluid and +1∕3 for an ideally plastic fluid. It is a dimensionless parameter available from the net-pressure history. Its positive values show prevailing influence of channelized propagation.
(A-CF) Impermeable rock: contrast-free propagation: For impermeable rock, the extreme case of contrast-free propagation ( Δ l = Δ u = 0 ) corresponds to the axisymmetric HF problem. The fracture grows as a circle of radius r �� * (t) = x �� * (t) = z �� * (t) . Again the self-similar solutions for a Newtonian (Savitski and Detournay 2002) and power-law (Linkov 2016a) fluids yield dependencies monomial in the time. As shown in Appendix 3, for the normalized average net-pressure the solution yields where b p = −n∕(n + 2) ; P av is that defined in Appendix 3; thus approximately P av ≈ 3∕(16 * n 3 ) , and for the middle value * n = 0.716 , this gives P av ≈ 0.511 . For a Newtonian fluid ( n = 1 ), b p = −1∕3 , P av = 0.479.
p av (t) = P av t b p , Pay attention that for the considered case of contrastfree propagation, the exponent b p of the net-pressure history is non-positive, being equal to −1∕3 for a Newtonian fluid and 0 for an ideally plastic fluid. Again this parameter is independent of the normalized thickness H ′′ , and it is a dimensionless characteristic of the net-pressure history. This again indicates that the exponent b p may serve as a key parameter of a net-pressure history. Its negative values show prevailing influence of contrast-free propagation.
(B-CH) Highly permeable rock: channelized propagation: For highly permeable rock, the extreme case of channelized propagation ( Δ l = Δ u = +∞ ) is studied for a Newtonian fluid by Nordgren (1972, Appendix 2 of the cited paper). The extension to arbitrary power-law fluid is given in Appendix 4. In particular, in terms of the normalized variables (2), the net-pressure history is where b p = 1∕(4n + 4) ; similar to W 0n in (11), P 0n depends m e r e l y o n t h e f l u i d b e h a v i o r i n d e x : Newtonian fluid, b p = +1∕8 ; P 0n = 0.478 ; for a perfectly plastic fluid, b p = +1∕4 , P 0n = 0.450. Note that, similar to small losses, the exponent b p is positive for the channelized propagation. Again, it is independent of the pay-layer normalized thickness H ′′ .
(B-CF) Highly permeable rock: contrast-free propagation: For highly permeable rock, the extreme case of contrast-free propagation ( Δ l = Δ u = 0 ) is promptly obtained following the same path as that used in Appendix 4. Omitting technical details, for arbitrary power-law fluid, the normalized average net-pressure is where as in (12), the factor P av depends merely on the behavior index; b p = −(3∕8)n∕(n + 1) . Note that, similar to the case of small fluid losses, the exponent b p is non-positive for the contrast-free propagation. As above, it is independent of the normalized thickness H ′′ .
We see that in all the limiting cases considered, the exponent b p in a monomial representation of a net-pressure history, is independent of the normalized thickness H ′′ . For a given normalized leak-off coefficient C �� leak , the exponent depends merely on the intensities of stress contrasts, being positive for channelized and non-positive for contrast-free propagation. This suggests employing the exponent b p in a monomial approximation of a net-pressure history as its key parameter in general case of arbitrary stress contrasts.
p �� av (t) = P av t b p ,

Key Parameter of Pressure History in General Case
The analysis of the limiting cases implies that the sign of the net-pressure exponent b p is indicative as regards the important feature of HF propagation. Remarkably, the exponent is positive for channelized propagation, and it is non-positive for contrast-free propagation. For our theme, of prime significance is that the exponent b p is available from an observed (experimental) pressure history. If the pressure grows in time on the major part of the pressure history, the sign of b p is positive. This indicates that the propagation is of channelized type. In the opposite case, when the pressure decreases in time, the exponent is negative what corresponds to the contrast-free propagation. Thus, we obtain one-glance knowledge on the important feature of fracture propagation: channelized or contrast-free.
In the first case, the stress contrasts are large enough to attempt their identification by solving the inverse problem. If a success, the identified stress contrasts allow us to obtain valuable information on the footprint of the fracture and its opening.
In the second case, at least one of the stress contrasts is zero or negative. The identification of such contrasts is impractical, because of a priori knowledge that the major part of the fracture surface is beyond the boundaries of the pay-layer. The identification, if attempted, becomes complicated, because a P3D model, serving to speed up iterations, is inapplicable. This justifies focusing on positive stress contrasts, for which the exponent b p is non-negative.

Connection with Intensities of Stress Contrasts
The analysis above has also revealed that the exponent b p depends merely on the normalized stress contrasts Δ l�� , Δ u�� and the normalized leak-off coefficient C �� leak . It is independent of the normalized thickness H ′′ of the pay-layer. Then Eq. (6) implies that, for a given leak-off parameter C �� leak , the exponent b p is entirely defined by the intensities of the lower R l and upper R u stress contrasts. Hence, plotting the curves R u (b p ) for the limiting cases of symmetric ( k asym = 0 ) and utmost asymmetric ( k asym = 1 ) contrasts provides the range, in which the actual stress contrasts are to be searched. Thus, we obtain a means to solve the two major problems: (i) to define the ranges within which stress contrasts may be identified; and (ii) to choose a proper initial guess for starting iterations of an identification process. This makes the exponent b p the key parameter of a netpressure history. Clearly, this parameter, being dimensionless, does not depend on units, in which the pressure is measured. In particular, it is the same in normalized and non-normalized variables.

Numerical Calculation of Key Parameter b p
We could see that in all the limiting cases considered, the analytical solutions are monomial in time, Thus it may be expected that monomial representation of a net-pressure history may serve as its simple and accurate approximation.
The key parameter b p of the approximation is available by presenting a net-pressure history p net (t) in log-log scale. It is the slope of the plot in the log-log representation. In view of its significance, we sketch its evaluation.
Let p net (t) be a net-pressure, experimental or theoretical. We consider the time interval well-beyond the time t H , at which the fracture reaches the pay-layer boundary (this time is defined by Eq. (5)

). A monomial representation of p net (t) on this interval is
The best-fit values of b p and A p are easily found by the least-squares method after representing the net-pressure p net (t) by its N values p j = p net (t j ) at N time instances t j > t H . Explicit equations for b p , a p = logA p and A p are given in Appendix 5, which also illustrates the accuracy of approximation (15) in cases, intermediate between the limiting.
Comment 3: In fact, the key significance of the slope parameter b p has been revealed yet in early 80-ies by Nolte and Smith (1981). Since then, its using has gained recognition as the Nolte-Smith slope analysis. It looks appropriate to cite the prophetic assertion by Nolte (1989, p. 304): "…the log-log plot of net-pressure vs time is a basis for interpreting pressures during fracturing." Our analysis above gives an interpretation in terms of stress contrasts.

Limiting Curves
Suppose we have an experimental net-pressure history. Its key parameter b p is found as explained in Appendix 5. The parameter depends merely on the intensities of the lower R l and upper R u stress contrasts and on the normalized leak-off coefficient C �� leak . This implies that the value b p may serve to make a reasonable initial guess on the intensities R l and R u of the contrasts. By (6) and (2), assigning R l , R u is equivalent to assigning the normalized Δ l�� , Δ u�� and physical Δ l , Δ u (15) p net (t) ≈ A p t b p stress contrasts. Thus it is sufficient to have an initial guess in terms of the intensities R l and R u . The discussion above suggests the following approach.
For a given leak-off coefficient C �� leak , by using the IP3D model, we calculate the dependences R u (b p ) for each of the two limiting cases, distinguished in Sect. 3.2.2: (I) symmetric contrasts: R u = R l ; k asym = 0 , and (II) utmost asymmetric contrasts: R l = ∞ , R u ≥ 0 ; k asym = 1.
In graphical representation [see below pairs (AI, AII) and (BI, BII) in Figs. 3, 4], we obtain two limiting curves. The actual values of R u are in the area between the curves. For the found b p , the true value of R u is within the interval between the curves for this b p . Thus a proper initial guess to start identification may use R u and R l from this interval.

Numerical Calculation of Limiting Curves
Calculations of limiting curves are performed by employing the IP3D model for the symmetric and utmost asymmetric cases. They involve merely two numerical values: the intensity parameter R u and the normalized thickness H ′′ . Recall that the theoretical analysis of Sect. 3.2.3 has revealed that in the limiting cases of channelized ( R u = ∞) and contrast-free ( R u = 0) propagation, the parameter b p is independent of H ′′ . Thus, by continuity, it may be expected that the influence of the normalized thickness H ′′ on b p is negligible for intermediate values of R u , as well. Calculations, performed for various values of H ′′ , specified below, confirm this suggestion. For certainty, we present results for the case of a Newtonian fluid; the results, obtained for thinning fluids ( 0 < n < 1 ), Fig. 3 Limiting curves R u (b p ) for symmetric (AI, BI) and utmost asymmetric (AII, BII) stress contrasts are quite similar, and lead to similar conclusions concerning with identification. With a fixed value of the normalized leak-off parameter C �� leak , the calculations were performed as follows.
The dimensionless intensity parameter R u was changed from 0.01 to 100. (Actually, as follows from the numerical results obtained, it is sufficient to consider a narrower interval 0 < R u < 10 .) In the logarithmic scale, logR u was changed from − 2 to 2. The cycle in values of R u employed 21 values of logR u : (logR u ) j = −2 + 0.2j ( j = 0, … , 20) . The normalized thickness H ′′ was assigned to include the practically significant range 1.5 < H ′′ < 15 . We used three values of H ′′ : 1.572; 4.715; and 14.14. The near middle value ( H �� = 4.715 ) corresponds to the input data of the example given in Fig. 9 of the paper by Dontsov and Peirce (2015). The values threefold less ( H �� = 1.572 ) and threefold grater ( H �� = 14.14 ) refer to the minimal and maximal bounds of the interval. Finally the numerical results were compiled in 6 tables (for the three values of H ′′ in each of the two cases considered). They contained the values b p , A p , a p = logA p , found in accordance with Appendix 5 as functions of the discrete values of R u .

Limiting Curves for Impermeable and Highly Permeable Rock
For illustration, we present results of numerical calculations for the limiting curves in the extreme cases, distinguished in Sect. 3.2.3, when the leak-off parameter C leak is either very small or very large: (A) impermeable rock, when fluid losses are negligible ( C leak = 0), (B) highly permeable rock, when fluid losses are large ( C leak → ∞).
Such calculations were performed for each of the cases (A) and (B), corresponding, respectively, to impermeable (AI), (AII) and highly permeable (BI), (BII) rock. In the second case, the value C �� leak , for which the propagation is entirely dominated by fluid losses, was C �� leak = 0.045 . For this and greater values of C �� leak , the calculated parameter b p coincides with its value in the analytical solution (13). As mentioned, the results were compiled in 6 tables. For brevity, we present merely their summary, rather than the tables themselves.
Numerical results for the exponent b p : Noteworthy, in agreement with the expectation, the results, obtained in the calculations, show that the exponent b p does not depend on the normalized thickness H ′′ . Specifically, for a given R u , the values of b p , calculated for various H ′′ , stay the same to the fourth decimal digit. (Naturally, the values b p are different in the cases of symmetric and utmost asymmetric contrasts). As mentioned, this property, established theoretically and confirmed numerically, simplifies using a net-pressure history for identification of stress contrasts. It also suggests an appropriate choice of the initial guess to start iterations of an identification process. Note that in difference with the exponent b p , the factor A p in the monomial representations depends on the normalized thickness H ′′ of a pay-layer.
Numerical results for the limiting curves: The calculated dependences between values of R u and b p are plotted as the limiting curves R u (b p ) . Figure 3 presents the pairs of such curves for impermeable (A) and highly permeable (B) rock. In accordance with the theory, each of the pairs has vertical asymptotics, corresponding to entirely channelized propagation. For the latter, in the case of a Newtonian fluid, the exact theoretical values of b p are, respectively, 0.2 and 0.125.
From the asymptotic behaviors, demonstrated by Fig. 3, it is clear that when the intensity R u exceeds four, the propagation is, in fact, completely channelized. Thus, it is surely channelized, when b p ≥ 0.18. Fig. 4 The dependence between the stress intensity and net-pressure history parameters for the cases of symmetric (AI) and utmost asymmetric (AII) stress contrasts. The solid line presents average values of stress contrast 1 3

Sound Initial Guess to Start Identification of Stress Contrasts
For certainty, illustrate a proper choice of a sound initial guess by considering the case of impermeable rock. Its limiting curves, shown by the dashed lines in Fig. 3, are reproduced in Fig. 4. They are complemented with the solid line R mid (b p ) of middle values. With a given b p , defined by (15) for a given net-pressure history, it is reasonable to start identification of R l and R u from the initial pair ( R u 0 , R l 0 ), where R u 0 = R mid , and R l 0 is taken from the upper curve. The corresponding starting values of physical stress contrasts are defined by (8). The normalized stress contrasts, defined by (6),

Ranges of Strong, Weak and Intermediate Stress-Barriers
From the asymptotics of limiting curves, it could be seen that when the minimal of contrast intensities R u exceeds four, the propagation is practically channelized. Thus, for a given normalized leak-off coefficient C �� leak , the corresponding value of b p defines the border, starting from which the propagation is entirely channelized. Then there is no sense in solving the inversion problem to find exact values of the stress contrasts. The inversion is actually impossible: iterations, if a success, may converge to any pair of R l and R u , exceeding four. Below we shall illustrate this effect numerically.
On the other hand, small values of b p correspond to practically contrast-free propagation, for which finding exact values of stress contrasts is also of low significance. Besides, in this case it is impossible to have accurate results by employing a P3D model, because the assumption of plane-strain state is violated.
Thus, solving the inversion problem for identification of stress contrasts appears reasonable merely in the intermediate region between the channelized and contrast-free propagation. Including into an analysis the factor A p in the monomial representation (15) serves to deepen understanding typical features of the characteristic intervals.
Numerical results for the factor A p : In difference with the exponent b p , the factor A p depends on the normalized height H ′′ of a pay-layer. Since the dependence on H ′′ is quite smooth, some general conclusions may be obtained for the fixed normalized height H �� = 6.926 . As mentioned, this value corresponds to the middle of an interval of practical significance.
The monomial dependence (15) is linear in log-log coordinates. Thus it is convenient to use a p = logA p , rather than A p itself. The dependences a p on b p for the symmetric and utmost asymmetric contrasts in the case of impermeable rock are given in Fig. 5 by the dotted and dashed curve, respectively. The areas of strong and weak contrasts, are shown by lower and upper ellipses. These areas have specific features.
The area of strong contrasts ( b p > 0.18 ) is characterized by coalescence of the curves. The coalescence conforms with the expected entirely channelized propagation along the pay-layer. Figure 5 gives quantitative estimation to this tendency. Actually, the cases of symmetric and utmost asymmetric contrasts become indistinguishable. Accordingly, as established, the inversion becomes impossible.
In the area of weak contrasts ( b p < −0.05 ), the curves are nearly parallel. Even for the same limiting value b p = −1∕3 , corresponding to both contrast-free and semicontrast-free propagations, the curves are far from coalescence. Clearly, the disagreement with the theory is due to inapplicability of a P3D model to these cases. Thus, as expected, the inversion by means of the IP3D model is impossible in the area of weak contrasts.
The central area, shown by the middle ellipse, defines the range of intermediate contrasts. It defines the range of the net-pressure parameter b p , in which identification by means of the IP3D model is sensible. Comment 4: From Fig. 4, it appears that the distinguished intervals agree with the intervals of strong, weak and intermediate intensity R , found in the papers Dontsov and Peirce (2015), Gladkov and Linkov (2018). Thus, while a bit conditionally, the earlier concepts of strong, weak and intermediate contrasts are now quantified via the key parameter of a net-pressure history b p .

Influence of Errors in Evaluation In Situ Stresses on Identification Stress Contrasts
In field, merely total pressure history p(t) is recorded. The net-pressure p net (t) is a theoretical quantity, defined as the difference between the total pressure and the in situ stress P in a pay-layer. The latter is estimated by specially designed experimental methods (e.g., Haimson and Fairhurst 1969), and an estimation has an error P . This error may notably influence the accuracy of identification stress contrasts via a net-pressure history. Its influence grows with the depth of oil/gas recovery. Denote 0 = P ∕ P the relative error of the in situ stress P found experimentally, = p net ∕p net the error of a net-pressure p net defined by (1), = p net ∕ P is the ratio of the net-pressure to the in situ stress. The definition (1) yields = − 0 1 . This implies that for large P , when 0 > , the relative error of the net-pressure exceeds 100%. At usual depths, the in situ stress P is of orders 10-100 MPa. In difference, the net-pressure p net , typical for a HF treatment, is notably less; it is of order 1 MPa. Thus the ratio is of orders 0.1-0.01. Then, if the relative error 0 of finding the in situ stress exceeds 10%, the relative error of the net-pressure becomes greater 100%, when the ratio is less than 0.1. Then the identification of the both stress contrasts via the net-pressure history becomes impossible for whatever perfect theoretical model and whatever perfect method of identification.
Illustrate this conclusion for the case when the recorded pressure (after an initial period of first minutes) is nearly constant. In this case, total and net-pressure histories are nearly monomial with the exponent b p close to zero ( b p ≈ 0 ). Then (15) yields p net (t) ≈ A p . The relative error A p ∕A p of the identified factor A p is A p ∕A p ≈ p net ∕p net = . Thus, if the error exceeds 100%, the error of A p is the same: this factor is actually unknown.
With uncertain factor A p , identification of the both stress contrasts becomes impossible. Merely the known in advance exponent b p may be identified correctly as b p ≈ 0.
Still, the minimal of the contrasts, denoted Δ u , may be estimated avoiding iterations. It is obtained by direct employing the limiting curves of Fig. 3. Specifically, when b p ≈ 0 , the lowest of the limiting curves provides the low bound of the minimal stress contrast parameter: R u ≈ 0.9 for both impermeable and highly permeable rock. By continuity, this estimation refers to an arbitrary leak-off. Then for a Newtonian fluid ( n = 1 , k R = 1.63 ), from (8) it follows that the minimal contrast Δ u is estimated as The maximal stress contrast Δ l stays uncertain. Its value may be arbitrary in the interval Δ u ≤ Δ l < ∞ . This uncertainty is caused by the mentioned uncertainty of the netpressure, arising from the error 0 in the measured in situ stress P .
Comment 5: In the case of nearly constant total pressure, an error in evaluation of the in situ stress P may be decreased by including it as an additional parameter to be identified and by employing the total pressure history p(t) . Actually, this improves the accuracy of finding the in situ stress. The initial (experimental) value of P may serve as a sound initial guess to start iterations. When employing this option, the arising uncertainties in the leak-off coefficient and in the stress contrast Δ l may be damped by assuming the rock impermeable ( C leak = 0 ) and the contrasts symmetric ( Δ l = Δ u ).

Influence of Uncertain Interpretation of Recorded Pressure
In practice, the pressure recorded at wellhead or at the bottom of a wellbore notably differs from the pressure in a fracture, which is present in equations describing the fracture propagation. Furthermore, during pumping there occur changes of the recorded pressure caused by unpredictable oscillations of frictional losses. Thus interpretation of a recorded history in terms of the pressure in the fracture is quite challenging. The influence of arising uncertainty in the identification of stress contrasts is similar to that of errors in assigning the in situ confining pressure. Recall that merely the net-pressure p net and the stress contrasts Δ , defined by (1), rather than the p and P themselves, are present in the HF equations. Since p net = p − P , the error in identification of stress contrasts Δ is completely defined by the error in the difference p − P . Therefore, the estimations and conclusions of Sect. 6.2 for the influence of errors in the in situ stresses refer, as well, to errors caused by the uncertain interpretation of recorded pressure. The only change concerns with the definition of the relative error 0 . The reference value P and its error P are replaced with a reference pressure p ref and its error p ref . Actually, since the difference p ref − P is small as compared with the values of p ref and P themselves, the relative error may be defined as 0 = p ref ∕ P .
The emphasized similar influence of the in situ stress P and the reference pressure p ref suggests also an extended interpretation of the results, when the in situ stress P is included as an additional identified parameter. Then its bestfit value Pbf provides the correction to the whole difference p ref − P , rather than to the value P used as an initial guess. Consequently, if an experimental value of P is found accurately, then its difference with Pbf gives a correction to p ref .
Vice versa, for a reliably estimated reference pressure, the difference gives a correction to the in situ stress P .
The influence of white noise, caused by unpredictable changes of frictional pressure, may be dampened, at least partly, by employing the standard identification procedure of the Kalman's filter (e.g., Catlin 1989; Bolzon et al. 2002). It may replace or complement the least-square method employed in numerical calculations of the next section. Surely, its use requires preliminary analysis of statistical characteristics of random fluctuations observed.
Fortunately, the slope parameter b P is quite insensitive to errors in the reference pressure p ref and in situ stress P . Thus, with dampened influence of the fluctuations, this parameter becomes known and may serve for rough estimations of stress contrasts.

Outline of Calculations
In this section, we illustrate options for identification of stress contrasts by numerical examples. The calculations were performed following the path described in Sect. 2.
Initial guess: Normally, the initial guess is taken in accordance with Sect. 5.4 as a sound choice suggested by the theory. For comparison, other choices, made by rule of thumb, are also employed for the inversions.
Examples used: Most of the examples below employ the input data of papers with numerical results for truly 3D modeling.

Examples of Notably Channelized Fracture Propagation
Consider the example, given in the paper Dontsov and Peirce (2015) and corresponding to the laboratory experiment (Jeffrey and Bunger 2009). The physical input parameters are: n = 1 , = 30.2 Pa s, E = 3.3 GPa, = 0.4 , H = 0.05 m, Q 0 = 1.7⋅ 10 −8 m 3 /s, Δσ u ex = Δσ l ex = 4.3 MPa. Normalizing (2) yields the three input parameters H �� = 6.926 , Δσ u�� ex = Δσ l�� ex = 0.242 . The corresponding dimensionless intensities are R u ex = R l ex = 1.686 . From Fig. 4, it appears that the corresponding parameter of the net-pressure history is b P = 0.106 . Thus, the example refers to an intermediate case, which is closer to the channelized propagation than to contrast-free. The footprints, given in Dontsov and Peirce (2015) for truly 3D and P3D models, are quite close, and their form, stretched along the x-axis, confirms that the propagation is indeed notably channelized. Our calculations of footprints for the both models, performed by explicit integration, are indistinguishable from those in Fig. 4 of Dontsov and Peirce (2015).
However, the net-pressure histories at the source for these two models do not agree as good as the footprints. They are shown in Fig. 6 by solid and dotted lines, respectively for truly 3D and IP3D models. Evidently, at the interval of the main interest (after first 5 min), the slopes of the curves differ. With developing channelized propagation, the pseudoexperimental net-pressure p pex (t) = p th t, Δ u ex , Δ l ex becomes systematically greater than the experimental p ex (t) . As a result, the net-pressure history parameter for the experimental history is merely b P = 0.050 against the mentioned value b P = 0.106 for the pseudo-experimental history. This suggests that the accuracy of identification cannot be high. Comment 1 explains this systematical disagreement by the boundary layer effect, caused by the different simulation of the injection sources in the two models. The difference in the parameters b P implies that merely rough estimation of the minimal stress contrast Δσ u may be obtained. The maximal of the contrasts Δσ l is uncertain; it may be arbitrary in the range from Δσ u to infinity.
The identification process does not presume that the stress contrasts are symmetric. Thus the best-fit contrasts, obtained in the iterations for the experimental history p ex (t) , may be asymmetric. In the example considered, the contrasts, obtained in iterations, are close to the utmost asymmetric case: the calculated best-fit values are Δ u bf = 3.27 MPa and Δ l bf = 122470 MPa. The best-fit history parameter is b P = 0.07.
Notably, the found minimal best-fit stress contrast Δ u bf = 3.27 MPa does not differ too much from the experimental value Δσ u ex = 4.3 MPa. The relative difference is 24%. In accordance with the theory, the found best-fit parameter b P = 0.07 is closer to the experimental value b P = 0.05 , than to the pseudo-experimental b P = 0.106 . Similarly, the bestfit net-pressure history p bf (t) (dotted line in Fig. 6) is closer to the experimental history p ex (t) (solid line), than to the pseudo-experimental (dashed line).
The convergence of iterations was controlled by using the pseudo-experimental history in accordance with Sect. 2.3. The iterations were performed starting from the initial guess suggested in Sect. 5.4: Δ u 0 = 4.99 MPa, Δ l 0 = 5.23 MPa. Figure 7 illustrates the convergence to the exact values Δσ u ex = Δσ l ex = 4.3 MPa of the stress contrasts; Fig. 8 shows the convergence to the pseudo-experimental net-pressure history p pex (t) on iterations 1-4. The curves for two successive iterations (5, 6) practically coincide with the pseudoexperimental history.
Figures 7 and 8 demonstrate that the initial guess is quite sensible. Note that a choice by rule of thumb (without using the theoretical considerations) commonly leads either to divergence, or to much greater number of iterations. Specifically, close examination has shown that the iterations diverge, when starting values Δ u 0 and Δ l 0 are beyond the vicinity Δ u 0 = 7 MPa, Δ l 0 = 7.49 MPa. Quite similar results are obtained for the examples of high symmetric stress contrasts, given in Figs. 19 and 20 of the paper Chen et al. (2020). The intensities of stress contrasts for these examples were R u ex = R l ex = 1.17; 2.34; 3.50; 4.66. Starting from R u ex = R l ex = 2.34, the fracture propagation was practically channelized. At the final time, the ratios length-to-height were 20, 40 and 50 for R u ex = R l ex = 2.34, 3.50 and 4.66, respectively. Our calculations for the truly 3D and IP3D models gave practically the same results. Remarkably, similar to the example of the channelized propagation, considered above, the net-pressure the truly 3D model was less than that for the P3D model; consequently again the netpressure history parameter b P was less for the experimental history. Again, merely the minimal of the stress contrasts could be estimated more-less accurately, while the greatest of them corresponded to the utmost asymmetric case. We do not present graphs for these examples, because, being similar to those on Figs. 6, 7 and 8, they actually resemble the features already discussed. Note only, that for large length-to-height ratios, the time expense exceeded 3 hours.

Example of Non-channelized Fracture Propagation
Consider the example of symmetric contrasts, given by Dontsov and Peirce (2015) as "extreme" (p. 132 of the paper cited) to the applicability of the P3D model. Recall that their model in the case of symmetric contrasts, is equivalent to the IP3D model (Linkov and Markov 2020).  Figure 4 yields that the parameter b P is quite small: b P = 0.0082 ≈ 0 . This confirms that the example refers to the case when the fracture propagations becomes non-channelized, and the applicability of a P3D model becomes dubious. Still, the footprints for planar and IP3D models, while having the length-to-height ratios not exceeding three, do not differ dramatically (Dontsov and Peirce 2015;Linkov and Markov 2020). Thus the identification is still possible.
The experimental, pseudo-experimental and best-fit histories are shown in Fig. 9 by solid, dotted and dashed lines, respectively. It can be seen that, when the time exceeds the initial interval of first five minutes, the histories are quite close. The corresponding dimensionless parameters b P for these histories are, respectively, b P = 0.0082 ; 0.00712 ; and −0.00006. This suggests that the accuracy of identified stress contrasts may be quite high. The calculations confirm this suggestion: in the example considered, the identified best-fit values are Δ u bf = 0.243 MPa, Δ l bf = 0.283 MPa. They are reasonably close to the exact contrasts Δσ u ex = Δσ l ex = 0.25 MPa.
The convergence was studied by using the pseudo-experimental history with identification starting from the initial guess, provided by Fig. 4, Δ u 0 = 0.221 MPa, Δ l 0 = 0.239 MPa. Figure 10 illustrates the convergence to the exact stress contrasts. Three iterations gave practically the exact values. Figure 11 shows the convergence to the pseudo-experimental net-pressure history p pex (t) ; even the first iteration gives the history indistinguishable from the second. The figures demonstrate perfect convergence when starting from the initial guess suggested by the theory. Again, a choice by rule of thumb (without using the theoretical considerations) does not provide such a favorable identification. Commonly it leads either to divergence or to much greater number of iterations.

Example of Notable Difference in Contrasts
In general, the contrasts are non-symmetric. To see how it influences the identification, we consider the previous example retaining one of the contrasts on its level Δσ u ex = 0.25 MPa, while taking the other contrast three-fold greater Δσ l ex = 0.75 MPa. Thus R u ex = 1.194 , R l ex = 3.582 . The coefficient of asymmetry is k asym = 0.5.
The footprints, calculated by solving problems for the truly 3D and IP3D models, have been already shown in Appendix 2 (Fig. 16). It could be seen that despite the assumptions of P3D models are not met, the footprints do not differ markedly. As expected, they are strongly asymmetric about the abscissae. The calculated net-pressure histories also do not differ too much. They are shown in Fig. 12. The corresponding dimensionless parameter b P for these histories is close to zero. It is, respectively, b P = −0.0068 (experimental), and b P = −0.0234 (pseudo-experimental). This suggests that the accuracy of identification may be satisfactory. The curve for the calculated best-fit history p bf (t) , shown by dashes, confirms this suggestion. As expected, it is close to the experimental history p ex (t) ; for it, b P = −0.0126 .
The identification was performed starting from the initial guess provided by Fig. 4: Δ u 0 = 0.177 MPa, Δ l 0 = 0.192 MPa. For the pseudo-experimental history, they rapidly  converge to the exact values of the stress contrasts to the accuracy of 1.3% in 4 iterations (Fig. 13). The net-pressure history becomes indistinguishable from the pseudo-experimental on the second iteration (Fig. 14).
The identified best-fit contrasts are Δ u bf = 0.275 MPa, Δ l bf = 0.502 MPa. In accordance with the theory, the minimal of them better agrees with the exact value Δσ u ex = 0.25 MPa, than the maximal with Δσ l ex = 0.75 MPa. The example demonstrates that rough identification by using the IP3D model is possible even when the suggestion that cross sections are in plain-strain state is obviously not met.

Summary
The study of the paper expounds to which extent and how a pressure history p(t) , obtained in a HF treatment, may be used to decrease uncertainty of stress contrasts at the layers adjacent to a pay-layer.
In brief, the conclusions are as follows. Employing the normalized variables and the three dimensionless parameters proposed noticeably simplifies an overall study of options for identification of stress contrasts via a pumping history. Their using provides ranges, within which the identification is possible in principle. They also serve to obtain a sound initial guess facilitating rapid identification within the ranges distinguished. Then identification of stress contrasts via a net-pressure history may be efficiently performed by means of the IP3D model. Numerical examples for pumping histories, corresponding to data of publications on modeling planar fractures, confirm the theoretical conclusions.
The authors clearly comprehend that pumping histories, observed in field, are influenced, sometimes drastically, by factors not present in the simplest model of a planar HF in homogeneous elastic rock. Still the model accounts for some of the factors of prime significance. Since the identification by means of the IP3D model is fast and nonexpensive, the method suggested may be routinely used to get an idea on an expected HF propagation.

Appendix 1: Identification by Least-Squares Method
Consider the problem of finding m unknown parameters C j ( j = 1, … , m ) through N experimental values P i ex ( i = 1, … , N ), found at N successive time instances  1, … , N ). Besides, we have a theoretical model, which provides the theoretical values P i th ( i = 1, … , n ) at the same instances t i : P i th (C) = p th (t i , C) ( i = 1, … , N ) for any set C of the parameters C j . In vector notation, the identified components C j are assembled as the vector of unknowns C C 1 , … , C m ; the experimental values P i ex as the experimental vector P (P 1 ex , … , P N ex ) , and the theoretical values p th (t i , C) as the theoretical vector P (C) . Commonly, P (C) is continuous in the components C 1 , … , C m .
For a given set C of identified parameters, the error of its using may be estimated by the Euclidean norm: Then, the problem of identification is reduced to finding such a best-fit vector C , for which the error reaches its global minimum. This vector certainly exists, because Err(C) is a continuous function, which being non-negative is bounded by zero. The question arises how to efficiently find the best-fit vector C ? Direct search by testing various combinations of components C 1 , … , C m is very time expensive even when the number of unknowns m is two only. The time expense may be drastically reduced when having an appropriate initial guess to focus the search to its vicinity. Further improvement may employ classical methods of analysis. Specifically, the necessary conditions for a local minimum are: In general, it is a system of m nonlinear equations in m unknowns C j . Its solving by the Newton-Raphson method (e.g., Epperson 2011) yields the iterations: Herein, vectors are arranged as columns, the superscript T denotes transposition, the upper symbol −1 denotes inversion; C k is the vector of unknown parameters, obtained on the previous (k-th) iteration; L k is the sensitivity matrix composed of the sensitivity coefficients L ij . The latter are defined by the partial derivatives calculated at the point C = C k : The iterations (16) are initiated by assigning an initial vector C 0 (initial guess). They are terminated when the norm (commonly Euclidian) of the difference C k+1 − C k becomes less than an accepted tolerance . In practical applications of this and other methods of identification, iterations may E rr C j = 0 j = 1, … , m. (16) diverge. In general, convergence strongly depends on a proper initial guess (e.g., Epperson 2011).
In the problem of identification stress contrasts, the experimental vector P is composed of N known values of the experimental pressure history p ex (t i ) at reasonably chosen N time instances t i . In fact, it is sufficient to use N = 20 instances on a time interval of interest (from first to tens or hundreds minutes). The vector P may refer to the net-or total pressure. In the first case, the identified vector C includes merely two stress contrasts ( m = 2 ): C 1 = Δ u , C 2 = Δ l (equivalently, C 1 = R u , C 2 = R l ). In the second case, the in situ stress 0 is included as an additional unknown C 3 = 0 to be identified; the number of unknowns is m = 3 . Formally, additional parameters, like the pay-layer height H, and/or leak-off coefficient C leak , and/or the elasticity modulus E , may be included into identified parameters. However, as clear from the examples given, this may notably complicate the choice of a proper initial guess and deteriorate the convergence of iterations.
On each iteration (16), the theoretical vector P (C k ) is calculated by means of the IP3D model. The sensitivity coefficients L ij , defined by (17), are also found by means of the IP3D model via central differences. As mentioned, the explicit Euler integration for this model performs quite efficiently. This makes the time expense for a single iteration (16) not exceeding 3 s. Thus, for a properly chosen initial guess, when the number of needed iterations is within ten, the identification may be performed in real time of hydraulic fracturing in field. heights, which are normally within the range 1.5 < H ′′ < 15 . The fracture propagation was traced by using the IP3D model numerically from 0.5 min to two hours with a quite rough spatial mesh. Even using 10 nodal points on the fracture half-length suffices good agreement of numerical results with the analytical solutions. Specifically, the exponent b p for the net-pressure history agrees with the analytical value = 1∕5 to the accuracy of 0.003%. The exponent b x for the half-length agrees with the analytical value b x = 4∕5 to the accuracy of 0.4%. As expected, for the rough mesh employed, the agreement of the dimensional factors A p and A x is a bit worse than that for the exponents. Still, in the entire range 1.5 < H ′′ < 15 , the calculated factors A p and A x agree with their analytical values 0.700 1 Hε 1.2 and 0.602 1 Hε 0.8 to the accuracy 4% and 3%, respectively.
Both stress contrasts are zero ( Δ l = Δ u = 0 ). Figure 15 illustrates the qualitative applicability of the IP3D model in this case, for which the main assumption of the model is strongly violated. It presents the calculated footprints for the truly 3D and IP3D models in the first quadrant for the case of zero stress contrast. The footprints correspond to the normalized height H �� = 6.926 of a pay-layer and for the large time (7373 s) of a HF treatment. They are plotted in the normalized variables (2).
As could be expected, there is distinct difference between the circular and nearly linear footprints for the two models considered. Yet, the areas of the footprints differ merely 11.4%. In view of the identically met mass balance, the same difference refers also to the fracture openings averaged over the footprints. The half-length x �� * and the half-height z �� * differ from the radius r �� * merely 12.7% and 17.3%, respectively. Therefore, despite drastic difference in the form of the footprints, the IP3D model provides reasonable estimations of the length, height and average opening of a fracture. Since the zero stress contrast refers to the case quite unfavorable as regards assumptions of the P3D model, the agreement justifies its application for rough quantitative estimations of the important quantities discussed.
The comparison of net-pressure histories is complicated by the fact that axisymmetric solution for a penny-shaped fracture driven by a Newtonian fluid goes to infinity at the source point = 0 (Savitski and Detournay 2002). In contrast, for any P3D model, the net-pressure at the source is always finite. Similar drastic disagreement occurs near the fracture front, where the net-pressure defined by the axisymmetric solution, goes to minus infinity, while for a P3D model, the net-pressure at the front is always zero. Hence, the net-pressure histories at the source point certainly disagree for the planar and IP3D models when stress contrasts are zero. Merely average values of net-pressures may be compared. Still, for the both models, the mass conservation is met identically at each time instance and, consequently, on any time interval. For zero leak-off, the volume of pumped fluid identically equals to the volume of the fracture. Thus, despite the differences, it may be expected that some integral characteristics of the solutions will behave similarly and may be estimated by using the IP3D model.
As an illustration, consider the results of the IP3D model for the normalized height H �� = 6.926 of a pay-layer (this value is near the middle of the interval mentioned) at a representative time of HF treatment (7373 s). The monomial approximation of the numerical results obtained by numerical integration yields for the average net-pressure with b p = −0.332 , A p = 0.457 . The comparison with the exact theoretical values b p = −1∕3 , P av = 0.479 shows that the exponents b p agree to the accuracy of 0.4%. This again indicates that the dimensionless parameter b p is a useful general characteristic of stress contrasts. Again the agreement between the dimensional factors A p is worse than for the exponents. Still, the difference between 0.479 in (12) and A p = 0.457 in (18) is merely 4.6%. Thus, using the IP3D model for modeling net-pressure histories looks acceptable even in the extreme case of zero stress contrasts. The calculated histories may be approximated by monomial dependences of (15) type.
The lower stress contrast is infinite; the upper is zero ( Δσ l = +∞ , Δ u = 0 ). As explained in the beginning of Sect. 3.2.3, in this case, the exact monomial solution for the fracture radius remains applicable with the additional factor 3 √ 2 in front of x �� * (t) . The monomial dependence p �� av (t) = 0.479t −1∕3 of the average net-pressure does not change at all.  (2015), where under the same remaining input data, the stress contrasts were assumed symmetric ( Δσ u ex ∕Δσ l ex = 1 ), now they are asymmetric ( Δσ l ex ∕Δσ u ex = 3). The footprints are shown in Fig. 16 for the time t = 132 min. Obviously, they do not correspond to the plane-strain conditions. Still, the footprints do not differ too much, although there are notable differences in details. The areas of the footprints are actually the same; due to volume conservation, average openings are nearly the same, as well. The example demonstrates that despite violation of the main assumption, the IP3D model may serve for rough estimations of the stress contrasts.
Again we can compare the monomial theoretical netpressure with the results of calculations employing the IP3D model. For the normalized height H �� = 6.926 and the time interval of 7373 s, used in the previous case, the calculations with employing the IP3D model give the exponent b p = −0.331 . Again the dimensionless parameter b p of the netpressure history differs from the exact theoretical value −1∕3 merely 0.7%. Again it appears that this parameter may serve as a key characteristic of the net-pressure history although the main assumption of a P3D model is violated.
As above, the agreement of the exact dimensional factor A p = 0.479 with the value A p = 0.574 , found numerically for the IP3D model, is worse than that for the exponents. The difference is about 20%. However, the comparison of the analytical fracture half-lengths for the axisymmetric problem ( x �� * (t) = z �� * (t) = 3 √ 2 ⋅ 0.697t 4∕9 ) and the half-length x �� * (t) , calculated by means of the IP3D model, shows better agreement. At t = 7373 sec, the calculated half-length x �� * = 47.28 differs from the theoretical value of radius r �� * = 3 √ 2 ⋅ 36.5 = 45.99 merely 2.7%. The height 2z �� * = 70.62 differs from the diameter 2r �� * = 91.98 some greater, 23.2%. Again it appears that, despite of essential differences between the models, the IP3D model is applicable for rough quantitative estimations of the important quantities discussed.

Appendix 4: Solution to the Nordgren's Problem for Large Leak-Off in the Case of Power-Law Fluid
Consider the Nordgren's problem (Nordgren 1972), extended to an arbitrary power-law fluid. We shall use the normalized variables (2). In the case of large leak-off, similar to the case of a Newtonian fluid, studied by Nordgren, the continuity equation becomes where q ′′ is the normalized flux through a cross section of the normalized height H ′′ , * r �� * (t) = t. The ODE (24) is solved with the Cauchy condition q �� (0) = 1∕(2H �� ) , expressing that, due to symmetry with respect to r �� = 0 , merely half of the normalized (unit) influx flows through a cross section. Omitting details, integration of (24)  Substitution (25) and (27)  Integration (28) under the condition of zero opening at the fracture front yields where b p = 1∕(4n + 4) . Using (29) in (27) gives the normalized net-pressure. Its value at the source point is For a Newtonian fluid ( n = 1 ), in the physical variables, Eqs. (29), (30) become the Nordgren's (1972) results for large fluid-loss rate (Appendix 2 of the cited paper).
Recall that in the three limiting cases R l = R u = ∞ ; R l = R u = 0 ; and R l = ∞ , R u = 0 , the exact solutions have monomial forms. In these cases, results, obtained numerically by employing the IP3D model, reproduce the solutions quite accurately. By continuity, the monomial approximation (31), (32) is accurate in intermediate cases of stress contrasts, as well.
The accuracy is illustrated by using the input data of the two examples given in the paper Dontsov and Peirce (2015), for which intensities of stress contrasts are far from the extreme values. One of them refers to the experimental study of a HF in laboratory (data for Fig. 4 of Dontsov and Peirce 2015); the other to field conditions (data for Fig. 9 of Dontsov and Peirce 2015). The both refer to a Newtonian fluid and symmetric stress contrasts. In the first example, after using the normalizing (2), the three normalized input parameters are Δ l�� = Δ u�� = 0.0242 , H �� = 6.926 ; the time of reaching the pay-layer boundary is t H = 36.7s . In the second example, Δ l�� = Δ u�� = 0.229 , H �� = 4.715 ; the time of reaching the pay-layer boundary is t H = 15.5s . In terms of the intensity parameters, in the first example, R l = R u = 1.686 ; in the second, R l = R u = 1.194 . Note that the intensities R l , R u , for the both cases, are of the same order, although they refer to quite different (even in orders) values of non-normalized physical quantities.
The comparison of evaluated net-pressure histories, calculated by employing the IP3D model, and their monomial approximations p P3D (0, t) = A p t b p is presented in Fig. 17. In the both cases, the maximal relative error of

Log(t j )
A p = 10 a p the approximation occurs at the time of about 7 min, when the fracture geometry still does not correspond to planestrain conditions. The errors are, respectively, 3.2% and 1.3% for the first and second examples. At the final times, the errors are, respectively, 1.6% and 0.27%. Therefore, the monomial approximation (31) is applicable at the major part of a pressure history.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.

Fig. 17
Calculated (solid lines) and approximated (dashed lines) net-pressures. Curves 1 and 2 correspond to, respectively, data used for Fig. 4 and 9 of the paper by Dontsov and Peirce (2015)