Two-dimensional numerical investigations on the termination of bilinear flow in fractures

Bilinear flow occurs when fluid is drained from a permeable matrix by producing it through an enclosed fracture of finite conductivity intersecting a well along its axis. The terminology reflects the combination of two approximately linear flow regimes: one in the matrix with flow essentially perpendicular to the fracture, and one along the fracture itself associated with the non-negligible pressure drop in it. We investigated the characteristics, in particular the termination, of bilinear flow by numerical modeling allowing for an examination of the entire flow field without prescribing the flow geometry in the matrix. Fracture storage capacity was neglected relying on previous findings that bilinear flow is associated with a quasi-steady flow in the fracture. Numerical results were generalized by dimensionless presentation. Definition of a dimensionless time that, other than in previous approaches, does not use geometrical parameters of the fracture permitted identifying the dimensionless well pressure for the infinitely long fracture as the master curve for type curves of all fractures with finite length from the beginning of bilinear flow up to fully developed radial flow. In log–log scale the master curve’s logarithmic derivative initially follows a 1/4-slope straight line (characteristic for bilinear flow) and gradually bends into a horizontal line (characteristic for radial flow) for long times. During the bilinear flow period, isobars normalized to well pressure propagate with the fourth and second root of time in fracture and matrix, respectively. The width-to-length ratio of the pressure field increases proportional to the fourth root of time during the bilinear period, and starts to deviate from this relation close to the deviation of well pressure and its derivative from their fourth-root-oftime relations. At this time, isobars are already significantly inclined with respect to the fracture. The type curves of finite fractures all deviate counterclockwise from the master curve instead of clockwise or counterclockwise from the 1 /4-slope straight line as previously proposed. The counterclockwise deviation from the master curve was identified as the arrival of a normalized isobar reflected at the fracture tip 16 times earlier. Nevertheless, two distinct regimes were found in regard to pressure at the fracture tip when bilinear flow ends. For dimensionless fracture conductivities TD < 1, a significant pressure increase is not observed at the fracture tip until bilinear flow is succeeded by radial flow at a fixed dimensionless time. For TD > 10, the pressure at the fracture tip has reached substantial fractions of the associated change in well pressure when the flow field transforms towards intermittent formation linear flow at times that scale inversely with the fourth power of dimensionless fracture conductivity. Our results suggest that semi-log plots of normalized well pressure provide a means for the determination of hydraulic parameters of fracture and matrix after shorter test duration than for conventional analysis.


Introduction
Transient fluid flow in fractures or faults plays an important role for the production of oil and gas, as well as for fresh water supply and the production of geothermal energy, especially from artificial fracture systems, so called hot-dryrock (HDR) or enhanced geothermal systems (EGS).Flow in fractures and fracture networks may as well be important for the triggering of seismicity by precipitation (e.g., Hainzl  et al., 2006), by groundwater recharge (e.g., Saar and Manga, 2003), by hydraulic stimulation (e.g., Deichmann and Ernst, 2009;Majer et al., 2007;Shapiro and Dinske, 2009), and by water level changes in dams (e.g., Chen and Talwani, 1998).
For fractures in an impermeable rock matrix, fluid flow and pressure propagation are restricted to the fracture volume, and are thus exclusively controlled by the hydraulic diffusivity of the fractures.In contrast, fluid flow and pressure propagation in fractures is accompanied by fluid exchange with a permeable rock matrix, a rather complex problem for mathematical treatment.A first analytical solution was presented in the context of well testing (Cinco-Ley et al., 1978) that applies to the case of fluid production from boreholes subsequent to hydraulic fracturing (for the simplified geometry of a single fracture aligned with a borehole, see Fig. 1a).Cinco-Ley et al. (1978) simplified the flow field as a superposition of two fields of parallel flow: one in the fracture and one in the rock matrix, the latter perpendicular to the fracture plane (see Fig. 1b).Accounting for its peculiar geometry, this flow regime was named bilinear flow (Cinco-Ley and Samaniego-V., 1981).Evidence for bilinear flow was reported from hydraulic tests after hydraulically fracturing a low permeable matrix; for example, in tight basins that produce gas (Rushing et al., 2005;Stright and Gordon, 1983) and in sedimentary and granitic geothermal reservoirs (Häring et al., 2008;Jung and Weidler, 2000;Ortiz et al., 2011;Zimmermann, 2006).Interest in unconventional gas recovery from tight formations also triggered studies considering horizontal wells (see, for example, Du and Stewart, 1995;Jelmert and Vik, 1995;Verga and Beretta, 2001).
For constant production, a bilinear flow field is accompanied by a decrease of the wellbore pressure proportional to the fourth root of elapsed pumping time.The time window, during which this fourth-root relation can be observed, is however finite, and thus long-term predictions -of great practical importance for exploitation of liquid or gaseous resources -are erroneous when using this relationship.Therefore, constraining estimates of the end time of bilinear flow received attention in previous research (Cinco-Ley and Samaniego-V., 1981;Weir, 1999).Since radial flow dominated by the matrix properties develops when this time is exceeded, it specifically marks the end of the gain due to a stimulation operation involving hydraulic fracturing.
Until today, the physical understanding of the proposed relations for the end time of bilinear flow is incomplete.In this study, we rely on numerical simulations using a twodimensional finite element model in order to investigate the hydraulic diffusion in finite conductivity fractures.We include an analysis of the spatio-temporal characteristics of the entire pressure field in fracture and matrix in order to clarify the flow processes that lead to the termination of bilinear flow as well as to substantiate quantitative rules for the end of bilinear flow.Outlining the end time and investigating the pressure field in a dimensionless parameter space allows for us to generalize our findings obtained for specific cases to fractures with a range of dimensionless fracture conductivities.Focus is put on fractures with negligible storage capacity, and wellbore storage is also neglected.However the formulation aims at clarifying the role of fracture length.In this contribution, we first briefly give the background in terms of governing equations and non-dimensional formulation, describe the chosen modeling approach, report results, and subsequently discuss them in the light of their practical use.

Governing equations for the hydraulics of a fractured well
For a well intersected by a single fracture and surrounded by a permeable matrix, two basic hydraulic equations, partial differential equations for fluid pressure p, have to be considered, namely for flow in the infinite, isotropic, and homogeneous matrix, and for flow in the fracture.The two equations are coupled by the fluid flow between matrix (subscript m) and fracture (subscript F), q F (t) (see, for example, Cinco-Ley et al., 1978).
Here The flow in the matrix is assumed to obey Darcy's law, q m = −k m ∇p/η f .Similarly, Eq. ( 2) describes the mass balance in the elements of the fracture, which has to account for an additional source term representing the fluid exchange between the matrix and the fracture (q F (t)).
Frequently, these equations are recast introducing the hydraulic diffusivities of the matrix, and the fracture, Chaudhry, 2004;Dake, 2001;Matthews and Russell, 1967).Specific solutions of the governing Eqs. ( 1) and (2) for particular initial and boundary conditions have led to the distinction of characteristic flow regimes.Radial flow, characterized by a well pressure changing proportionally to the logarithm of elapsed pumping time, results when production from (injection into) a homogeneous formation causes radial flow lines to develop normal to the production surface composed solely of the well.For homogeneous and isotropic media, diffusion of pressure perturbations obeys a linear scaling relation between the square of the characteristic propagation distance L c and the characteristic propagation time t c involving the hydraulic diffusivity of the matrix, i.e., D m ∼ L 2 c /t c (see, for example, radius of investigation or drainage in Bourdet, 2002;Chaudhry, 2004;Dake, 2001;Earlougher, 1977;Horne, 1995;Matthews and Russell, 1967).Highly permeable fractures -that is, fractures in which the pressure gradient is negligible -intercepting the well may extend the effective production surface such that flow in the subsurface is actually directed towards this extended surface rather than radial towards the well (e.g., Jenkins and Prentice, 1982).Such flow geometry is termed formation linear flow since straight flow lines are thought to result in the matrix.The well pressure changes proportionally to the square root of elapsed pumping time.The bilinear flow regime is encountered when the flow is approximately linear in both the fracture or narrow zone of high conductivity and the matrix (e.g., Butler and Liu, 1991).In this regime the finite conductivity of the fracture leads to a finite pressure gradient in the fracture (Boonstra and Boehmer, 1986;Gringarten, 1985).Fractures in an impermeable matrix result in fracture linear flow, which is per se indistinguishable from formation linear flow regarding the power law relation between well pressure and elapsed pumping time.

Modeling approach
In our study, we focus on fractures with negligible storage capacity S F → 0; that is, the fracture is considered undeformable, and the amount of fluid in the fracture is considered small enough for its compressibility to be neglected.In this approximation, Eq. ( 2) reads that is, the pressure in the fracture obeys an inhomogeneous diffusion equation with a time-dependent source term but without an intrinsic transient term.The changes of pressure with time are considered to be dominated by the fluid transfer between matrix and fracture.Previous studies revealed that only for small dimensionless times do the flow in the well and the flow in the matrix deviate from each other due to storage effects in the fracture (e.g., Weir, 1999).In the classical dimensionless form of (2) (see Eq. A1 in Cinco-Ley and Samaniego-V., 1981), the transient term on the left side appears multiplied by the diffusivity ratio κ = D m /D F .When this ratio is small, the effect of the intrinsic transient term on wellbore pressure becomes important only for extremely small values of time, and thus can be neglected for bilinear flow (see Cinco-Ley and Samaniego-V., 1981;Riley, 1991).

A. E. Ortiz R. et al.: Two-dimensional numerical investigations
We use a two-dimensional finite element model consisting of a fracture with half-length x F positioned on the x axis (at y = 0; see Fig. 1c) and intercepting a well along its axis.Flow in the fracture of finite conductivity embedded in the permeable matrix is approximated as one-dimensional and wellbore storage is neglected.Thus, (1) and (3) are implemented as and The origin of the coordinate system coincides with the well, actually represented by a point source with a flow rate q w / h determined from the true flow rate in the well q w and the height of the open well section h (Fig. 1).The fluid flow between matrix and fracture, couples the two equations, where the factor of 2 accounts for the communication via the two fracture surfaces.
In principle, numerical analysis does not require the flow geometry in the matrix to be prescribed, as do the majority of previously presented analytical solutions.The assumption of most analytical treatments that flow lines in the fracture and in the matrix remain strictly perpendicular to each other indeed cannot hold towards the end of bilinear flow.The pressure diffusion in the matrix proceeds proportional to t 1/2 , ultimately surpassing pressure diffusion in the finite conductivity fracture that scales with t 1/4 .Thus, eventually isobars have to change direction with increasing time.
We performed more than 30 simulations using COMSOL Mulitphysics®employing the linear system solver UMF-PACK.To ensure the occurrence of bilinear flow, fracture length was varied from 1.5 to 1500 m, while the further parameters remained constant (q w / h = 2 ).An effect of the model boundaries on the simulation results was avoided by locating them far from the fracture (about 600 to 1200 m).The used free-meshing technique generated an unstructured mesh with triangular elements (Fig. 1d).Close to the fracture a finer mesh was obtained by fixing the length of the fracture elements to less than 0.36 m.Thus, the total number of elements varied considerably depending on the model size determined by fracture length (with order of magnitude of about 10 4 to 10 5 ).The time step was gradually increased from 1 s to 3600 s during a simulation.A typical succession consisted of a time step of 1 s for the first 30 s, then 20 s until 600 s, 60 s until 12 000 s, 300 s until 72 000 s, and 3600 s until the end of the simulation.However, the time step was also reduced for specific phases of a simulation to study the completion of bilinear flow in detail.Initially this "zooming in" into time required repetition of a certain simulation, while the gained understanding of the systematic relations for the occurrence of bilinear flow eventually allowed for the time intervals to be set for reduced time steps right from the start of a simulation.

Dimensionless formulation
Our numerical modeling is performed with dimensional properties, but for reporting results we use non-dimensional parameters in order to foster a fundamental understanding of bilinear flow in our conceptual study.Previous analyses of flow regimes employed a variety of non-dimensionalization approaches.Here, we use the conventional definition for dimensionless pressure (see, for example, Earlougher, 1977;Matthews and Russell, 1967): However, we use a modified definition of dimensionless time: where t D = tD m /x 2 F is the classical definition of dimensionless time for an infinite reservoir adapted for the flow in fractures by replacing radial distance r with half-fracture length x F (Cinco-Ley et al., 1978;Earlougher, 1977;Matthews and Russell, 1967).Furthermore, the dimensionless fracture conductivity is defined by The employed model parameters correspond to values of T D ranging from 0.1 to 100.Our choice of non-dimensional parameters is guided by the necessity to avoid fracture storage capacity and the request to also avoid fracture length.When fracture length is used as an explicit parameter in the conventional definition of dimensionless time, one encounters the problem that "time" becomes ill defined for very long or infinitely long fractures.The formulation should, however, be apt for fractures with a range of finite lengths, such as those created, for example, during hydraulic fracturing operations, as well as for fractures with "infinite" length, such as those encountered when length simply exceeds the influence zone of the pumping operation.The latter situation may rather be typical for stimulations in a geothermal context that create a connection between the well and either an extended network of natural fractures or a large geological fault.

Previously presented solutions for bilinear flow
The theoretical background of bilinear flow was first presented by Cinco-Ley et al. (1978) and Cinco-Ley and Samaniego-V.(1981), who studied the solution of ( 1) and (2) using a two-dimensional numerical model and Laplace transform, respectively.They demonstrated that the pressure in a vertically fractured well producing at constant flow rate is proportional to the fourth root of time in the bilinear flow regime.This result, subsequently confirmed by several approaches (e.g., Riley, 1991), reads for the non-dimensional well pressure during bilinear flow with our set of non-dimensional parameters.Thus, we get a unique relation between non-dimensional well pressure and non-dimensional time independent of any further model parameters.Some approximate analytical solutions for the pressure distribution in infinitely long fractures were derived in previous studies (see Boonstra and Boehmer, 1986;Weir, 1999).Notably, Boonstra and Boehmer (1986) demonstrated that during a certain sequence of bilinear flow, the pressure distribution is governed by a single variable combining time and distance (w in their notation).Weir (1999) subsequently emphasized the self-similarity of the pressure function with x 4 /t in contrast to Theis' solution for a homogeneous reservoir (also called line source or exponential-integral solution, Earlougher, 1977;Matthews and Russell, 1967;Theis, 1935) that admits r 2 /t as a self-similar variable.
In their seminal study, Cinco-Ley and Samaniego-V.(1981) report three expressions for the end time of bilinear flow, but did not explicitly state how these relations were derived.Essentially, two extended regimes separated by a short intermediate regime are found for end time as a function of dimensionless fracture conductivity.

Results
Presentation of the results of our numerical simulations first focuses on the evolution of well pressure.We then continue with the propagation of the pressure perturbation along the fracture and in the matrix in order to identify the mechanisms for the end of bilinear flow based on which we establish rational criteria for the end time.

Evolution of the well pressure
Following common praxis, results are presented in the form of type curves (in log-log scale) of the dimensionless well pressure and its dimensionless logarithmic derivative versus dimensionless time (Fig. 2a, b).The main features can be best explained using the derivative (Fig. 2b).In accord , where p ∞ wD is the master curve for an infinitely long fracture (see Table 1).and Samaniego-V., 1981;Weir, 1999), the derivatives follow a straight line with a slope of 1/4 over over a certain period of dimensionless time.For sufficiently high dimensionless fracture conductivities (i.e., T D 1), the derivative first turns counterclockwise into a straight line with slope ½corresponding to formation linear flow, fully developed only for T D > 50 (Fig. 2b).Ultimately, as expected (e.g., Bourdet, 2002), derivatives bend into a unique horizontal line (dp wD /d ln τ = 0.5) for all fracture conductivities, indicating that radial flow is reached.The dimensionless time to reach fully developed radial flow increases with decreasing T D , and is highest for T D = 0.The log-log plots of the well pressure and its derivative (Fig. 2a, b) suggest that type curves with T D > 1.8 and T D ≤ 1.8 bend off counterclockwise and clockwise from the 1/4-slope straight line, respectively, as previously described by Cinco-Ley and Samaniego (1981).However, introducing a normalized well pressure, p wD /2.45τ 1/4 , as a measure of the deviation from the expected bilinear behavior (10), the resulting presentation is more sensitive than the conventional type curves of the well pressure, and shows that the curve for T D = 0 actually constitutes the master curve followed by all type curves of normalized pressure for a certain time interval (Fig. 2c).The master curve (Table 1), addressed as p ∞ wD in the following, is associated with an infinitely long fracture (note that the alternative case for T D = 0, namely T F = 0, is meaningless since in this case no fluid can be injected or withdrawn via the fracture).The normalized pressure stays close to unity until dimensionless time τ > 10 −6 , when the master curve bends downwards with increasing slope (clockwise), indicating the transition from bilinear to radial flow.
All normalized type curves first deviate counterclockwise from the master curve before finally bending clockwise like the master curve itself.The latter behavior is evidenced for a few curves, but is outside the explored time range for others (Fig. 2c).Normalized type curves for T D > 10 start bending counterclockwise in the section where the normalized master curve is still at unity; those for T D < 10 only in the downward bending part of the master curve.As a consequence, normalized type curves for T D > 10 exhibit a maximum, whereas those for 1.6 < T D < 10 exhibit a succession of a minimum and a maximum.For T D 1.6 the two extrema degenerate to a single saddle point, and all normalized type curves with T D < 1.6 decrease monotonically.Normalized type curves with 1.8 < T D < 10 intersect the horizontal line corresponding to unity twice, whereas those with T D < 1.8 stay below.This behavior is probably the reason why previous authors, e.g., Cinco-Ley and Samaniego (1981), assumed a discontinuity in the behavior of the type curves near T D = 1.8.

Evolution of the pressure along the fracture and in the matrix
For a systematic analysis of the evolution of the pressure field in the fracture and in the matrix, a ratio of pressure differences, is defined, where p(x, y, t) denotes the pressure at position (x, y) in the fracture or matrix, p 0 the initial pressure (assumed identical in matrix and fracture), and p w (t) the well pressure at time t.Thus, the quantity p N compares the change in pressure at some point in the fracture or matrix to the pressure change in the well.Lines in the (x, y) plane with p N = const.are referred to as normalized isobars in the following (see examples in Fig. 3).The ratio of pressure differences notably assumes identical values when calculated using either absolute or dimensionless pressures.Normalized isobars exhibit a characteristic flame-like, almost triangular shape with a pointed tip before reaching the fracture tip, and a semi-circular front after passing it, (see Fig. 3) rather than running parallel to the fracture as assumed in previous treatments prescribing flow lines in the matrix perpendicular to the fracture (e.g., Cinco-Ley et al., 1978).The numerical simulation, however, also shows that after the start of injection or production the normalized isobars migrate with the fourth root of dimensionless time along the x D axis and with the square root of time along the y D axis for a certain time (Fig. 4).During this time period the dimen-sionless distance of the normalized isobars on the x D axis progresses according to with the value of the constant α b depending on the chosen isobar; for example, α b is about 3 and 2 for the normalized isobars p N = 0.01 and 0.05, respectively.In dimensional variables, Eq. ( 12) reads where x i (t) is the position of the isobar in the fracture for time t and here referred to as the bilinear flow diffusivity.This diffusivity combines fracture and matrix properties, and has dimensions of L 4 /T.Equations ( 13) and ( 14) are specific formulations of the self-similarity of the pressure profiles in the fracture during bilinear flow found by Weir (1999).
The dimensionless distance of the normalized isobar on the y D axis is given by during bilinear flow corresponding to in dimensional variables.Combining ( 12) and ( 15), the evolution of the ratio y iD /x iD depends on dimensionless time as The migration of the normalized isobars according to Eqs. ( 12) to ( 17) terminates for two different reasons depending on the size of the dimensionless fracture conductivity.The change in migration behavior occurs in the interval 1 < T D < 2, and we illustrate the two types of terminations by considering two examples -T D = 0.314 and T D = 5 -in Fig. 4.
For dimensionless fracture conductivities lower than 1, the normalized isobars start to slightly accelerate relative to the fourth-root-of-time migration along the x D axis long before they reach the fracture tip, i.e., x D = 1 (Fig. 4).Migration of the normalized isobars in the y direction simultaneously slows down a little bit relative to the initial square-root-oftime migration (Fig. 4a, b).The curves of x iD and y iD merge close to the interception of the two extrapolated diagnostic fourth-root and square-root relations actually occurring at (τ 1, x D = y D 1).After merging, the two curves follow the ½-slope straight line, indicating that radial flow conditions are approached.For dimensionless fracture conductivities larger than 2, migration of the normalized isobars along the x D axis decelerates and actually almost terminates for a finite time interval when reaching the fracture tip and long before the interception of the two extrapolated diagnostic fourth-root and square-root relations (Fig. 4c, d).After some time, the migration finally accelerates again and appears to approach the straight line of y iD that closely follows a square-root-of-time relation.Upon closer inspection, one notices that x iD slightly accelerates just before the prominent halt in migration associated with the arrival at the fracture tip, i.e., x iD = 1 (Fig. 4c,  d).This intermittent acceleration is caused by the reflection of the isobar at the fracture tip that one can rationalize when invoking an image fracture at the fracture tip and an image well at a distance of 2x F producing or injecting at the same rate as the real well.The "reflection" of the normalized isobar is then approximated by the superposition of the pressure fields of the two wells.Migration will accelerate when the isobars of the two wells are approaching each other from both sides of the fracture tip (Fig. 5).A more detailed view of the deviation of the normalized isobars from the fourth-root and square-root-of-time behavior is obtained by using normalized presentations, x iD /T D τ 1/4 and y iD /T D τ 1/2 , similar to the one used for the well pressure (Fig. 2c).These presentations (Fig. 6a, b) confirm that normalized x iD is constant for a certain time interval, and starts to bend upward in a similar way as the normalized well pressure bends downward (Fig. 2c).Normalized y iD is not constant even in the early stage (τ 10 −6 ), but decreases continuously with a slight increase in slope at dimensionless time τ 10 −2 (Fig. 5b).The early deviation from the square-root-of-time migration indicates that even in the direction perpendicular to the fracture, the pressure propagation is affected by the presence of the fracture at all times.The width-to-length ratio of the normalized isobars (or pressure field), y iD /x iD , initially follows a fourth-root-of-time relation (Fig. 6c), and subsequently bends clockwise from the 1/4-slope straight line simultaneously with the upward bending of normalized x iD .Within the resolution of our numerical simulation, this ratio is almost identical for all normalized isobars (Fig. 6c), suggesting that all normalized isobars have a similar shape and undergo the same evolution simultaneously in dimensionless time.The observed relationship between the ratio y iD /x iD and dimensionless time τ allows for us to determine the width-to-length ratio of all normalized isobars at any instant.During the bilinear flow period this ratio is approximately given by y iD /x iD 1.1τ 1/4 (Fig. 6c).

End time of bilinear flow from well-pressure observations
For single-well tests, well pressure constitutes the only observable pressure, and the end time of bilinear flow is generally determined by using its deviation from the 1/4-slope straight line (in log-log plots of well pressure vs. time).According to our observations (Sect.3.1), the type curve for the infinitely long fracture rather than the 1/4-slope straight line represents the master curve, and only in its initial part up to dimensionless time τ 10 −6 is this master curve identical to the 1/4-slope straight line.For later times the master curve bends clockwise from the 1/4 -slope straight line due to the gradual transition from bilinear to radial flow.All type curves for fractures with finite length deviate counterclockwise from this master curve.Type curves for dimensionless fracture conductivity T D < 10 do this in the clockwise-bending section of the master curve, and those for T D > 10 in the straight-line section.We consequently introduce two criteria for the termination of bilinear flow: one for the clockwise deviation of the master curve from the 1/4slope straight line, and a second one for the counterclockwise deviation from the master curve.According to the underlying mechanisms, these criteria are addressed as transition criterion and reflection criterion, respectively.The time determined by the transition criterion will be addressed as transition time, and the time determined by the reflection cri-  12) and ( 13) and calculated from Fig. 4. terion as reflection time.In order to achieve higher accuracy in the determination of transition time and reflection time, the two criteria are formulated for the normalized type curves p wN = p wD /2.45τ 1/4 and p ∞ wN = p ∞ wD /2.45τ 1/4 .The transition criterion addresses the clockwise deviation of type curves and the master curve from the horizontal at unity; that is, actual pressures fall short of the bilinear relation by a relative amount of ε due to the transition to radial flow (Fig. 7).Unaffected by the fracture tip, the type curves under consideration are identical with the master curve before the transition criterion is fulfilled.Accordingly, the transition time is identical for all type curves complying with this criterion, and in particular depends not on dimensionless fracture conductivity T D but instead only on the value of ε used (Fig. 8a, c).For each ε a maximum T D , however, exists up to which the transition time can be determined.For type curves deviating from the master curve before the transition criterion is met, the transition time cannot be determined.The reflection criterion reflects the described reflection of normalized isobars at the fracture tip (Sect.3.2) that produces the counterclockwise deviation of the normalized type curves p wN from the normalized master curve p ∞ wN (Fig. 7).Our data show that for normalized type curves with T D ≥ 1, the reflection at the fracture tip is strong enough to produce intersections with the curved deviation lines ( 19) up to at least ε = 0.05, and thus a reflection time τ r can be determined.Type curves with T D < 1 also show a reflection, but the associated deviation from the normalized master curve remains quite small, and the reflection criterion may not be met for any ε of practical significance.The reflection time is proportional to T −4 D for all type curves with T D > 1 (Fig. 8).This relation is intuitively understandable when recalling our observation that the normalized isobars in the fracture migrate proportional toτ 1/4 .The time it takes for a normalized isobar to propagate from the well to the fracture tip is therefore proportional to x 4 F , and since T D is inversely proportional to x F , the observed relation between reflection time and dimensionless fracture conductivity results.For dimensionless fracture conductivities T D < 1, this relation may no longer be valid since in these cases migration of the normalized isobars starts to accelerate relative to the fourth-root-of-time migration long before the reflection criterion is fulfilled.
The arrival time of a normalized isobar p N at the fracture tip is smaller than the reflection time by a factor of 16 (Fig. 8a, c) when the same value is used for ε and p N (e.g., ε = p N = 0.05).This numerical relation can be explained by turning again to the above-introduced concept of an image fracture and an image well.In this concept, the reflection of the normalized isobar at the fracture tip is approximated by the superposition of the normalized pressure profiles of the two wells.Inserting 2x FD instead of x FD in Eq. ( 12) increases the time by the observed factor of 16.
During the bilinear flow period the dimensionless well pressure is proportional to τ 1/4 , and thus its values are constant for the transition time and proportional to T −1 D for the reflection time (Fig. 8b, d).Dimensionless well pressures for the reflection time and for the time of arrival of the normalized isobar at the fracture tip differ by a factor of 16 1/4  =2.The end time of bilinear flow reported by Cinco-Ley and Samaniego-V.(1981) differs significantly from our time estimates for dimensionless fracture conductivities T D < 5 (Fig. 8).When reporting their three regimes, Cinco-Ley and Samaniego-V.(1981) were not specific on their criterion for quantitative determination.Their data likely reflect the shortcomings encountered when relying on a deviation from the 1/4-power relation without investigating the deviation in detail, especially for type curves with 1.6 < T D < 2.5; end time apparently becomes a discontinuous function of T D as reported in Cinco-Ley and Samaniego-V.(1981).

Discussion
Introducing the dimensionless time τ according to (8) as an alternative to the approach by Cinco-Ley and Samaniego-V.(1981), while leaving all other dimensionless parameters consistent with this previous analysis, proved to be favorable for a better understanding of bilinear flow.The new dimensionless time permits for identification of a unique function of the dimensionless well pressure for an infinitely long fracture (T D = 0) applicable from the beginning of bilinear flow  18) and ( 19), respectively, and dimensionless arrival time τ a of the indicated isobar at the fracture tip for (a) ε = p N = 0.01 and (c) ε = p N = 0.05.The dimensionless end time reported by Cinco-Ley and Samaniego (1981) is represented by thin grey lines.Fit equations for the calculated end and arrival times are represented with blue lines and blue dashed lines, respectively.Note that τ r and τ a differ by a factor of 16.Dimensionless well pressure at τ e and τ a .Note that the dimensionless well pressures for τ r and τ a differ by a factor of 2. Fit equations for the calculated dimensionless well pressure are represented with blue lines and blue dashed lines by τ = τ e and τ = τ a , respectively.up to fully developed radial flow.A normalized presentation of computed type curves for various dimensionless fracture conductivities showed that the type curve for the infinitely long fracture rather than the 1/4-slope straight line, as assumed by Cinco-Ley and Samaniego (1981), constitutes the master curve for the type curves of all fractures with finite length.The latter all deviate counterclockwise from the master curve instead of clockwise (for low fracture conductivity) or counterclockwise (for high fracture conductivity) from the 1/4-slope straight line as assumed by these authors.
The master curve itself starts its clockwise deviation from the 1/4-slope straight line at a dimensionless time τ 10 −6 .However, clockwise deviation builds up very slowly and becomes noticeable in the commonly used log-log presentation only at τ 10 −2 .Furthermore, this clockwise deviation is counteracted by a counterclockwise bending for finite fractures.These counteracting effects are balanced best for fracture conductivities close to 2 so that these type curves stay close to the 1/4-slope straight line the longest, explaining why the end times of bilinear flow given by Cinco-Ley and Samaniego (1981) strikingly peak for T D ∼ 1.8.
In order to distinguish the two types of processes that lead to a termination of bilinear flow, transition to radial flow, and isobar reflection at the fracture tip preceding the transition to intermittent formation linear flow, we replaced the term "end time" by the specifications "transition time" and "reflection time" and established corresponding criteria.Application of these criteria to normalized type curves is especially advantageous for fractures with dimensionless conductivity between 1 and 2 since in this interval both transition time and reflection time can be determined.The transition time that is independent of dimensionless fracture conductivity can be determined for dimensionless fracture conductivities up to about 2; reflection time, inversely proportional to the fourth power of dimensionless fracture conductivity, can be determined for dimensionless fracture conductivities down to about 1.
Investigating the migration of isobars for (effectively) infinite fractures, we found that isobars normalized with respect to the well pressure migrate proportional to τ 1/4 along the x axis (i.e., in the fracture) and approximately proportional to τ 1/2 along the y axis (i.e., in the matrix perpendicular to the fracture) up to a dimensionless time τ 10 −2 .Their widthto-length ratio is close to 1.1τ 1/4 during this time period, and is independent of dimensionless fracture conductivity.This relation approximately holds for all normalized isobars with p N < 1 that therefore have similar shape at any instant.When τ is known, this simple relation allows for determination of the width-to-length ratio of the isobars at any instant.The ratio is 0.035 for τ 10 −6 , when the first sign of deviation from the fourth-root-of-time behavior of the well pressure is noticeable relying on the semi-log plots of normalized well pressure (Fig. 2c, Fig. 7).For τ 10 −2 , when the deviation becomes noticeable in conventional log-log presentations, the ratio amounts to 0.35.Thus, in a strict sense, bilinear flow ends at a very early state of the shape evolution, but the termination becomes noticeable in conventional type curves only at a much later state, i.e., at a time when the growth of the width-to-length ratio also starts to deviate noticeably from the fourth-root-of-time relationship.Shape evolution of normalized isobars was computed up to dimensionless time τ 1, where the width-to-length ratio is about 0.7 and apparently needs about two orders of magnitude more time to approach 1 for true radial flow in the case of the infinite fracture, a suggestion that has to be checked by further studies.
For finite fractures the shape evolution is disturbed when the isobars approach the fracture tip by a process here described as "reflection".This reflection is noticed by a clockwise deviation of the well pressure from the master curve at a time 16 times later than its actual occurrence at the fracture tip.Interestingly, the disturbance in the shape evolution shortens the time it takes to approach radial flow conditions characterized by a width-to-length ratio of 1.This shortening can be explained by the fact that migration of the isobars along the x axis is retarded after passing the fracture tip.
Our results substantiate and quantify the analytical finding of the self-similarity of the pressure functions during bilinear flow with the variable x 4 /t (see Boonstra and Boehmer, 1986;Weir, 1999).Thus, scaling considerations have to use D b ∼ l 4 c /t c with the hydraulic diffusivity for bilinear flow, D b , given by Eq. ( 14) rather than the traditionally used scaling relation between hydraulic diffusivity D hyd ∼ l 2 c /t c (e.g., Matthews and Russell, 1967) and characteristic length (l c ) and timescale (t c ) of pressure diffusion when a prominent two-dimensional hydraulic conduit (e.g., fracture or fault) is present.Due to these differences in scaling behavior, the characteristic propagation distance of the pressure perturbation in the vicinity of a fracture-like heterogeneity will exceed the classical penetration estimate for short dimensionless times, while it will increase much less than the classical penetration estimate for long dimensionless times.Some aspects of the conclusions based on hydraulic scaling relations as in analyses of hydraulic stimulation (e.g., Shapiro and Dinske, 2009) or aftershock (e.g., Mukhopadhyay et al., 2011) and swarm activity (e.g., Hainzl, 2004) may therefore have to be reconsidered in cases.

Implications for well-test analysis
Our findings may be used to determine geometric and hydraulic properties of fracture and matrix from short injection and production tests for which conventional well-test analysis would not be applicable.However, this evaluation requires excellent test conditions and high-quality pressure data.We recommend analyzing plots of normalized well pressure (Fig. 6) in addition to the conventional log-log plots and derivative (exemplified in Fig. 2 a,b).Using these normalized plots one may obtain the desired information on fracture transmissibility, matrix permeability, and fracture length at a much earlier time than with conventional procedures as explicitly outlined below.
The first step of any analysis is to determine the slope M characterizing the bilinear flow section of a diagram of the change in well pressure vs. fourth root of time ( p w = Mt 1/4 ).The slope is then used to construct a p w /Mt 1/4 vs. time diagram.In case bilinear flow ended during the pumping operation (indicated by either a clockwise or an counterclockwise deviation from a horizontal line at 1), transition time (t t ) and/or reflection time (t r ) and the corresponding well pressure (p wt and/or p wr ) can be read from this diagram for a considered relative deviation, i.e., a specific value of ε in ( 18) or ( 19).The following cases and their potential for determining fracture and matrix characteristics have to be distinguished: Case 1: When the well-pressure record constrains only the slope M, then the product T F (k m s m ) 1/2 can be determined using Case 2: When the well-pressure record exhibits a clockwise deviation of relative magnitude ε (that is, the transition time is known), then the permeability of the matrix k m can be determined from where p wtD = 0.38, 0.68, and 0.92 for ε = 0.01, 0.03, and 0.05, respectively (Fig. 8b, d).In case the storage coefficient the storage coefficient of the matrix, s m , is known or can be reasonably estimated, the fracture conductivity can be derived by Then, also the two diffusivities of the system, the one for bilinear flow D b ( 14) and the one for the matrix D m , are constrained.Furthermore, the ratio y iD /x iD (or y i /x i ) can be determined from The dimensionless fracture conductivity obeys the relation T D < T D max with T D max = 2.5, 1.6, and 1.0 for ε = 0.01, 0.03, and 0.05, respectively, and thus one also has a constraint on fracture length, i.e., x F > T F /T D max k m .Case 3: When the pressure record contains a counterclockwise deviation of relative magnitude ε but no clockwise deviation, then one knows the reflection time and has the relation T D > 10.The dimensionless fracture conductivity T D cannot be further constrained, however, since all type curves rapidly rise in a similar way.Thus, only the product T F (k m s m ) 1/2 can be determined (as in case 1).In addition, if the matrix properties (k m and s m ) are known or can be reasonably estimated, one can infer from the reflection time, where C = 1.73, 1.41, and 1.25 for ε = 0.01, 0.03, and 0.05, respectively (Fig. 8a, c).
Case 4: When the pressure record exhibits a clockwise deviation from the 1/4-slope straight line succeeded by a counterclockwise deviation from the master curve p ∞ wN , then transition time as well as reflection time are known.Such data allow for determination of matrix permeability, fracture transmissibility, and the ratio y i /x i (as in case 2).In addition, the dimensionless fracture conductivity T D can be quantified by looking for a pair of matching transition times and reflection times in Fig. 8, and with that the fracture length x F can be determined according to

Conclusions
Using two-dimensional numerical modeling, we investigated the evolution of the pressure field in and around a fracture imbedded in a permeable matrix during injection or production tests at a constant rate in a borehole aligned with the fracture.The understanding of the well-pressure evolution has gained significantly from introducing a new dimensionless time containing only the transport parameters of fracture and matrix as well as the storage coefficient of the matrix but no geometrical or storage parameters of the fracture.In this presentation, type curves of dimensionless well pressure for fractures with finite length evolve from a single master curve when dimensionless time progresses.The unique master curve corresponds to an infinitely long fracture and comprises two stages with an extended transition in between.The early and the late stage are characterized by pressure in the well increasing with time to a power of 1/4 (bilinear flow) and the logarithm of time (radial flow), respectively.For fractures of finite length, well pressure always deviates from the master curve towards higher pressures; that is, all type curves branch off counterclockwise from the master curve instead of clockwise or counterclockwise from the 1/4-slope straight line as considered by Cinco-Ley and Samaniego-V.(1981).
Nevertheless, two mechanisms have to be distinguished for the termination of bilinear flow depending on fracture and matrix properties.
For any fracture of finite length, the propagation of the pressure front in the fracture will eventually be affected by the fracture tip.Fractures with a dimensionless conductivity T D > 10 qualify as fractures with high conductivity since for these the reflection of the pressure front at the fracture tip happens long before substantial migration of isobars in the matrix.The reflection leads to a reduction of the pressure gradient in the fracture, and thus signals the transition to formation linear flow.Termination of bilinear flow is noticed by an increase of well pressure relative to the horizontal section of the normalized master curve that occurs, however, only 16 times later than the actual reflection at the fracture tip.In contrast, for fractures with low conductivity (T D < 1), migration of isobars in the matrix becomes significant long before the pressure front in the fracture approaches the fracture tip due to the difference in the power in the relation with time, i.e., square root and fourth root for matrix and fracture, respectively.The gaining of pressure propagation in the matrix on that in the fracture ultimately results in radial flow, indicated by normalized well pressures falling below the horizontal at unity and width-to-length ratios of isobars deviating significantly from an initial relation to the fourth root of dimensionless time and indicating substantial inclination of isobars with respect to the fracture.For an intermediate range of fracture conductivities (1 < T D < 10), reflection at the fracture tip interferes with the transition to radial flow and normalized well pressure exhibits a peculiar succession of decrease, increase, and decrease in cases.
The two criteria introduced for the deviation of the mastercurve from the fourth-root-of-time behavior (transition criterion) and for the deviation of the type curves for finite fractures from the master curve (reflection criterion) revealed that the transition time is independent of the dimensionless fracture conductivity and applies to the infinite fracture as well as to all finite fractures whose type curves do not branch off from the master curve before this end time is reached.The reflection time is inversely proportional to dimensionless fracture conductivity to a power of 4 corresponding to the fourth-root-of-time migration of the normalized isobars in the fracture expressed by a scaling relation that includes a bilinear diffusivity with dimensions of L 4 /t.
The gained insight into the relation between the entire flow field and the peculiarities of the recorded wellbore pressure permits for constraining hydraulic and geometrical parameters of the subsurface in practice.Using semi-log plots of   normalized well pressure in addition to the common log-log diagrams improves the sensitivity of analyses in particular for dimensionless fracture conductivities smaller than 3, and hydraulic parameters of matrix and fracture may be determined after shorter test duration than necessary for conventional analysis.
Our results substantiate and quantify the previously discussed self-similarity of the pressure functions during bilinear flow.In the presence of planar conduits with higher conductivity than the enclosing matrix pressure, diffusion obeys a scaling between characteristic propagation distance and characteristic time leading to faster and slower pressure propagation for small and large dimensionless times, respectively, in comparison to the classical scaling for hydraulic diffusivity in homogeneous media.Therefore, care has to be taken when using scaling arguments for situations comprising preferred hydraulic pathways, as, for example, generated during stimulation, or likely associated with earthquakes.
, k m [m²] denotes the matrix permeability, s m [Pa −1 ] the specific storage capacity of the matrix, η f [Pa s] the fluid viscosity, T F = b F k F [m³] the fracture conductivity with fracture width b F [m] and fracture permeability k F [m 2 ], and S F = b F s F [m Pa −1 ] the fracture storativity with specific storage capacity of the fracture s F [Pa −1 ] (see also the nomenclature at the end of the paper).Equation (1) expresses the mass balance in the volume elements of the matrix.Storage of fluid in a matrix volume associated with pressure transients, i.e., −s m ∂p/∂t, constitutes the sole source term, and thus equals the divergence of the flow through a volume element, ∇ • q m .

Fig. 2 .
Fig. 2. Type curves of the dimensionless well pressure (a), its logarithmic derivative (b), and of the normalized well pressure (c) as functions of dimensionless time.Note that (c) is actually a restricted zoom of the data presented in (a) and (b).For example, the clockwise bending of the curve for T D = 100 prominent in (a) is outside of the chosen scale.The red line in (c) represents the normalized curve p ∞ wD /2.45τ 1/4, where p ∞ wD is the master curve for an infinitely long fracture (see Table1).

Fig. 3 .
Fig. 3. Snapshots of the normalized pressure field p N for fractures with dimensionless fracture conductivity of T D = 0.314 (top) and T D = 5 (bottom).

Fig. 4 .
Fig. 4. Dimensionless distances x iD and y iD of isobars of normalized pressure p N along the x D and y D axis as a function of dimensionless time for (a) T D = 0.314, p N = 0.01; (b) T D = 0.314, p N = 0.05; (c) T D = 5, p N = 0.01; and (d) T D = 5, p N = 0.05.

Fig. 5 .
Fig. 5. Evolution of normalized distances x iD /T D τ 1/4 (a) and y iD /T D τ 1/4 (b), and of the ratio y iD /x iD (c) as functions of dimensionless time τ for normalized isobars p N = 0.01 and p N = 0.05.Dimensionless fracture conductivity is T D = 0.314 in all cases.

Fig. 6 .
Fig. 6.Evolution of normalized distances x iD /T D τ 1/4 (a) and y iD /T D τ 1/4 (b), and of the ratio y iD /x iD (c) as functions of dimensionless time τ for normalized isobars p N = 0.01 and p N = 0.05.Dimensionless fracture conductivity is T D = 0.314 in all cases.For τ < 10 −2 the normalized distance x iD /T D τ 1/4 is approximately constant and equal to α b introduced in Eqs.(12) and (13) and calculated from Fig. 4.

Fig. 7 .
Fig. 7. Type curves of the normalized well pressure for indicated values of dimensionless fracture conductivity T D .The dashed grey lines are the deviation lines representing the end time using the transition criterion (18); the solid grey lines are the deviation lines representing the reflection criterion (19).

Fig. 8 .
Fig. 8. Dimensionless end times τ e according to the transition mechanism (τ t ) and reflection mechanism (τ r ), Eqs.(18) and (19), respectively, and dimensionless arrival time τ a of the indicated isobar at the fracture tip for (a) ε = p N = 0.01 and (c) ε = p N = 0.05.The dimensionless end time reported by Cinco-Ley and Samaniego (1981) is represented by thin grey lines.Fit equations for the calculated end and arrival times are represented with blue lines and blue dashed lines, respectively.Note that τ r and τ a differ by a factor of 16.Dimensionless well pressure at τ e and τ a .Note that the dimensionless well pressures for τ r and τ a differ by a factor of 2. Fit equations for the calculated dimensionless well pressure are represented with blue lines and blue dashed lines by τ = τ e and τ = τ a , respectively.
of fracture during bilinear flow, Eq. (14), [m 4 s −1 ] D F hydraulic diffusivity of (isolated) fracture,D F = T F /η f S F [m 2 s −1 ] D m hydraulic diffusivity of matrix, D m = k m /η f s m [m 2 s −1 ] h heightof the open well section, fracture height, [m] dimensionless well pressure (Figs.2c and 6) [-] p weD dimensionless well pressure at end time of bilinear flow, [-] p wN normalized well pressure, i.e., normalized by the 1/4-relation for bilinear flow, [-] q w flow rate in the well [m 3 s −1 ] s m specific storage capacity of matrix [Pa −1 ] S F storativity of the fracture [mPa −1 ] t time, [s] T F fracture conductivity (transmissibility), [m 3 ] T Ddimensionless fracture conductivity, Eq. (9), [-] x, y spatial coordinates along, normal to the fracture with origin at the well,[m] x F fracture half length, [m] x D , y D dimensionless coordinates (x D = x/x F , y D = y/x F ) [-] x iD ,y iD dimensionless distances of normalized isobars from the well (along the x D -and y D axis, respectively) [-] Greek symbols α b constant for pressure diffusion in fracture during bilinear flow, Eq. (12), [-] α m constant for pressure diffusion in matrix, Eq. (15), [-] gamma function [-] η f fluid viscosity, [Pa s] τ dimensionless time, Eq. (8), [-] τ a dimensionless arrival time (of the normalized isobar at the fracture tip) [-] τ e dimensionless end time of bilinear flow [-] τ r dimensionless reflection time (arrival time of the reflected normalized isobar at the well) [-] τ t dimensionless transition time [-]