Investigations into the opening of fractures during hydraulic testing using a hybrid-dimensional flow formulation

We applied a hybrid-dimensional flow model to pressure transients recorded during pumping experiments conducted at the Reiche Zeche underground research laboratory to study the normal opening behavior of fractures due to fluid injection. Two distinct types of pressure responses to flow-rate steps were identified and numerically modelled using a radial-symmetric flow formulation for a fracture that comprises a non-linear constitutive relation for the contact mechanics governing reversible fracture surface interaction. These two groups represent radial-symmetric and plane-axisymmetric flow regimes from a conventional pressure-diffusion perspective. A comprehensive parameter study into the sensitivity of the applied hydro-mechanical model to changes in characteristic fracture parameters revealed an interrelation between fracture length and normal fracture stiffness that yield a match between field observations and numerical results. Fracture stiffness values increase with corresponding fracture length. Decomposition of the acting normal stresses into a stresses associated with the deformation state of the global fracture geometry and the contact stresses indicates that geometrically induced stresses contribute the more the lower the total effective normal stress and the shorter the fracture. Separating the contributions of the local contact mechanics and the overall fracture geometry to fracture normal stiffness indicates that the latter, the geometrical stiffness, constitutes a lower bound for total stiffness; its relevance increases with decreasing fracture length, too. Our study demonstrates that non-linear hydro-mechanical coupling can lead to vastly different hydraulic responses and thus provides an alternative to conventional pressure-diffusion analysis that requires changes in flow regime to cover the full range of observations.


Introduction
Estimation of a reservoir's effective hydraulic properties requires a consistent analysis of experimentally determined pressure and flow transients [22,23]. For individual fractures, simple analytical models for pressure diffusion have been applied when their intersection with boreholes classified them as axial or radial [27,41,40,28]. Analytical models based on solutions of the diffusion equation for fluid pressure for constant flow-rate steps document distinct differences in pressure response for the one-dimensional and the radial flow associated with axial and radial fractures, respectively. Rocks with a dense array of randomly oriented fractures may justify their treatment as porous media, leading to radial flow, too. Mathematically, the full range of responses can be addressed by regarding the dimension of the flow to be a parameter [42]. However, hydro-mechanical phenomena, such as reverse water-level fluctuations in distant monitoring wells [29,30] or insensitivity of pressure responses to increases in flow-rates, so called jacking [43] cannot be reproduced by pressure diffusion models and result in inaccurate approximation of the effective fracture characteristics [31,34,39]. Despite the growing number of treatments of hydro-mechanical coupling [34,36,37,38,6], the understanding of the influence of basic geometrical and mechanical properties of fractures on their hydraulic response to flow-rate or pressure perturbations appears still limited.
Non-Local fracture deformations triggered by perturbations of the fluid pressure along a fracture induce changes in permeability and volume of fractures with a direct impact on flow and storage characteristics and therefore on how the perturbations evolve with time and spread in space [33,31,45,44]. Accounting for hydro-mechanical interaction throughout numerical fitting of pressure p and flow-rate Q transients is a non-trivial task and requires consistent evaluation of the balance equations in an efficient manner. Evaluation of fracture opening or closing in response to a perturbation of the equilibrium state requires to consider the acting normal stresses owing to their control on the mechanical interaction between the fracture surfaces in contact. For example, large numbers of single Hertzian contacts have been invoked to characterize the mechanical interaction of two mated fracture surfaces [17,18,16]. Responses of these contacts to changes in shear and normal stress result in changes of the effective fracture aperture [19,21]. Fracture opening does not depend on local contact mechanics alone but also on the geometrical stiffness of the fracture [34]. Despite the importance of fracture stiffness for the interpretation of pumping operations, little work has been devoted to decompose these two contributions. Here, we analyze pressure transients from pumping tests conducted at the Reiche Zeche underground research laboratory, where the injection borehole penetrates a fractured rock. From a classical pressure-diffusion perspective, the hydraulic responses of the tested intervals mimic that of either radial-symmetric (positive-tangent group) or axial-symmetric (pressure-plateau group) fractures. Yet, logging and impression-packer results do not support this simple association of fracture geometry and hydraulic response. We employ a hydro-mechanical model considering radial-symmetric conditions, as applying for a radial fracture following a monolithic numerical implementation [6] to study the origin of the distinct pressure transients. Specifically, we studied the sensitivity of the hydro-mechanical model to the variation of characteristic fracture properties to identify best fits to the field data and the interdependence between the characteristic fracture properties. We highlight the influence of the two contributions to normal stiffness on the opening behaviour by separating the opening controlled by geometrical and normal contact stiffness as a function of the applied fluid pressure. We demonstrate that -with appropriate sets of parameters-distinctly different pressure responses can in principle be explained by a single model for hydro-mechanical effects, which contrasts the necessity to advocate differences in fracture orientation for pressure diffusion models.
2 Test site and experimental approach 2

.1 The Reiche Zeche underground research laboratory
As part of the research program of STIMTEC, a cooperative project investigating the creation and growth of fractures in crystalline rocks to develop and optimise hydraulic STIMulation TEChniques [46,47] we performed hydraulic tests in the research mine Reiche Zeche (Rich Mine), Freiberg (Germany). The average overburden at the test site amounts to about 130 m. The foliation of the fine-to medium-grained biotite gneiss dips 5 to 15°in south-east-direction. The gneiss is penetrated by fairly randomly oriented joints with an average separation of several decimeters. Fracture counting on retrieved cores yield 4.4±2.5, but intact sections with a length of 1 to 2 m occur. In the test volume of about 40 m × 50 m × 20 m, two to three steeply dipping, east-west trending damage zones were identified with a variable width between decimeters and a few meters.
The injection borehole BH10 with a length of 63 m and a radius of 0.038 m has a strike of N31°E and a dip of 15°f rom the horizontal and thus the borehole axis intersects the foliation at an angle of 20 to 30°. Ultrasonic transmission of the test volume as well as laboratory experiments on cores revealed a pronounced anisotropy in elastic parameters. Ultrasonic waves travel almost two times faster in the direction of the foliation than perpendicular to it. Dynamic and static Young's moduli in the two directions differ by 10 to 15 %, with the low modulus observed for loading perpendicular to the foliation [48].

Experimental procedure
Experiments were performed with a double-packer probe of Solexperts GmbH, Bochum, Germany, consisting of two inflatable packers isolating an injection interval of about 0.7 m length. The probe is equipped with uphole and downhole pressure gauges, and an uphole flowmeter, all sampled with 0.2 s. Flowrates measured uphole, i.e., outside of the borehole at the pump, were corrected for the storage capacity of the injection system to derive the true flow into the rock. The storage capacity was determined in a calibration experiment, for which the probe was inserted in a hollow steel cylinder.
The uniformly applied pumping protocol comprised a sequence consisting of a) injection (with rates of 2-10 l/min) until breakdown pressure was reached, the fracking, and subsequent shut-in phase, b) three repeat injections, the refracs, with moderate rates of 3 l/min at maximum, each again followed by a shut-in phase, and c) step-rate tests involving several phases of injection with constant flow rates, successively increased from below 1 l/min to about 5 l/min. The pressure response in these step-rate tests constrains the jacking pressure, the fluid pressure at which the fracture(s) intersecting the borehole wall open. Opening is indicated by a significant increase in injectivity, the ratio between flow-rate and pressure. Impression-packer tests were performed after the entire pumping sequence to document fracture traces on the borehole wall.

Intervals and selected data sets
The data used here represent part of the results of the step-rate tests performed in six intervals at depths of 24.6 m, 40.6 m, 49.7 m, 51.6 m, 55.7 m, and 56.5 m. Logging before and after fluid injection with an acoustic televiewer and impression-packer tests revealed evidence for pre-existing and induced fractures with a range of orientations (Table  2.3). We consider the circumferential fracture traces to represent radial fractures even though they do not intersect the borehole axis at a right angle as strictly required. Also, the traces classified as "axial" do not match this end-member geometry in a strict sense but the tilt to the borehole axis typical of en-echelon hydro-fractures occurring in boreholes that do not follow a principal stress axis [49]. Actually, the short traces of interval 51.6 m are not well constrained at all. Furthermore, it is impossible to associate the observed pressure transients with a specific fracture trace when intervals exhibit multiple traces. This situation is not unusual but representative of what an interpreter typically faces when tests are performed in crystalline rocks. classification of data set (see Figure 1) ‡ time it took for a pressure pulse to decay by 1/3 before the stimulation phase ⊺ time it took for pressure to decay by 1/2 during the shut-in phase after the step-rate test We selected three to five of the first low flow-rate steps for the six intervals, addressed as data sets M p a to M p c and M t a to M t c ( Table 2.3, Figure 1). The selection aimed to restrict to pressure and flow-rate couples, for which the proposed elastic model most likely applies. For some intervals, seismic activity was observed and therefore its absence during flow-rate steps could used as criterion for "elastic" response. The pressure transients induced by the step-wise increase of flow rate differ for the six intervals. We distinguish two groups of pressure evolution during a step. Pressure responses with flat tangents are evident in data sets M p a to M p c , a subset of our data that we will address as "pressureplateau group". In contrast, the data sets M t a to M t c exhibit continuously increasing pressure, the "positive-tangents group". For either group, however, the sensitivity of pressure level to flow rate diminishes with increasing flow rate, the observation interpreted as jacking.   Figure 1: The step-rate test data selected from the six intervals is divided in the pressure-plateau group labeled with M p and the positive-tangent group labeled with M t . The dark blue lines are associated to the recorded pressure data, the dark green lines show the step-wise increasing flow-rate records and the dotted grey lines represent tangents to the pressure transient at the end of a flow-rate step.

Numerical Method
When the aim is characterization rather than modification, hydraulic testing of fractures is performed below critical pressures for fracture extension, e.g., indicated by a decrease in injection pressure during constant-rate pumping, when rapid addressed as breakdown often related to tensile hydro-fracturing or when occurring over extended time periods or in a succession of small drops probably related to shearing events, either possibly accompanied by characteristic seismic activity. For tests performed at moderate injection pressures, fracture length can thus be treated constant and effects of changes in shear stress and therefore shear stiffness can be neglected. Fractures induced by hydraulic fracturing are expected to be oriented normal to the direction of the least principal (compressive) stress [34] so that they intrinsically fulfill the assumption of negligible shear stress. In the context of hydraulic characterization of fractures, their contact mechanics may thus be reduced to an account of their normal stiffness [12]. Hence, we apply a hydromechanical model considering the constitutive relation of normal fracture surface contact following a monolithic numerical implementation [6] for the numerical fitting of characteristic fracture parameters, initial width/aperture δ 0 , length l Fr , and normal stiffness parameter of a fracture E Fr . Its central aspects are recapitulated in the following, before details of the performed parameter search are presented.

Governing Equations
Flow processes of weakly-compressible, viscous fluids in high-aspect ratio fractures motivate the assumption of creeping flow conditions between two locally parallel plates, for which the balance of momentum reduces to a Poiseuilletype formulation [11,8], i.e., the relative fluid velocity w f is proportional to the pressure gradient grad p. In our continuum description, the associated cubic law is locally evaluated in the fracture domain Γ F r , i.e., on the level of a material point P(x, t), where x denotes its position vector, u(x, t) = x − X the fracture deformation relative to the reference position vector X, t time, and η fR the dynamic fluid viscosity. The locally evaluated, deformation-dependent permeability is identified as k s F r (x, t) = δ 2 /12 considering δ(u(x, t)) to be the effective fracture aperture.

4
The hybrid-dimensional formulation is obtained by inserting the balance of momentum into the balance of mass, derived for a deformable fracture. The outcome of a dimensional analysis of the resulting partial differential equations suggests that quadratic and convective terms are negligible [8,9,32] and the accordingly reduced hydro-mechanical governing equation reads ∂p ∂t comprising a transient I), a diffusion II), a coupling III), and a leak-off term IV ), where β f denotes the fluid compressibility, and q lk leak-off, i.e., the flow-rate from the fracture into the surrounding rock mass. The deformation dependent effective fracture aperture δ(u(x, t)) contributes to the characteristic diffusion process by term II) and to volume changes of the fluid domain by term III), which strongly couples the solution of the fracture-flow domain to the deformation state of the surrounding matrix.
The rock matrix surrounding the fracture might be treated by purely elastic or by biphasic poro-elastic (e.g., Biot's theory [1]) formulations depending on the application in mind. For the typically substantial difference in the characteristic times of pressure diffusion in the fracture and in a surrounding crystalline rock, a biphasic description results in oscillations of the pore-pressure solution, when time discretization and material properties are chosen in the relevant range to model the conducted field experiments. Hence, this work refrains from treating the matrix by Biot's full theory but approximates the material behavior with Gassmann's low frequency result [2,3].
The intact gneiss exhibits a permeability < 10 −20 m 2 [48]. Thus, leak-off from a fracture, into which fluid is injected from a borehole, into the "surrounding" is controlled by its intersection with other fractures. The hydraulic testing in BH10 revealed that the pre-existing fractures in the gneiss exhibit vastly variable hydraulic properties, as for example evidenced by the results of the pressure-pulse tests ( Table 2.3). We thus face a range of possible scenarios for the induced or pre-existing fractures intersecting the borehole. They may intersect only poorly permeable pre-existing fractures or linking up with a highly permeable pre-existing fracture. We consider either scenario to be suitable for an approximate description that neglects leak-off. For the second scenario, our modeling will simply gain the equivalent properties of a single fracture, since a variation of properties along a fracture is not tackled and thus a distinction of "individual" fractures composing a conduits is not possible. Neglecting leak-off likely overestimates "effective" length because all of the injected fluid volume has to be stored in the fracture. Applying a single fracture model with fixed geometry to this mix of fractures intends to test the versatility of the model and to determine equivalent fracture properties in a consistent way.

Constitutive Relations
Traditionally, normal contact models are expressed in terms of fracture deformation relative to the position corresponding to the first, stress-free contact of the two fracture surfaces, where fracture closing is described by positive deformation values [21,15]. For the response of the fracture to changes in normal stress, we use a modified non-linear elastic constitutive relation based on the model proposed by [15,14] that characterizes the normal deformation U e of an interface with two parameters, the initial stiffness at vanishing normal stress, E Fr eq , and the maximum displacement for infinite stress U max . To be consistent with the governing flow eq. (2), we formulate (3) in terms of relative aperture changes where the fracture deformation U e is defined as changes of the hydraulic aperture δ relative to the initial mechanical aperture δ mech 0 and the maximal deformation U max is defined with respect to the difference between the minimal mechanical fracture aperture δ min , approached for infinite normal stress, and the initial aperture, obviously obeying δ mech 0 > δ min . Neither absolute values of nor changes in mechanical and hydraulic apertures of fractures do have to coincide; particularly true once contact is established and the effective values of these aperture measures strongly depend on contact details and percolation characteristics in the fracture plane [13]. The resulting hydraulic, respectively mechanical fracture aperture might then be expressed in terms of where we simplistically relate initial mechanical and initial hydraulic fracture apertures by δ mech 0 = δ 0 /s 0 . Introducing the constant, dimensionless parameter s 0 ≥ 1 intends to distinguish pore space accessible for fluid flow characterized by the hydraulic aperture and the mechanical response of the contact surface characterized by the mechanical aperture [13,12,20]. Inserting (4) into (3) gives In principle, coupled hydro-mechanical simulations of deformable fractures require to numerically determine the equilibrium state of a fracture before the perturbation of its mechanical state, here associated with the pumping operations. Instead, we reformulate (6) using the aperture δ eq that reflects the unperturbed in-situ normal stresses σ Fr N,eq : i.e., we shift the reference state of the fracture to the in-situ stress level. Simple manipulations yield the relation in its implemented form where ∆δ eq = δ − δ eq defines the relative aperture change regarding the equilibrium fracture aperture and E Fr eq is introduced as the normal stiffness parameter of the equilibrium state.
The reduction of the numerical model to a single fracture interacting with the low permeable surrounding gneiss resulting in negligible exchange between fracture and solid domain similar to undrained conditions, negligible contribution of shear forces due to the low viscosity of the pore fluid, water, and by the low frequency range (≪ 100 Hz) of the perturbations induced by the step-rate tests motivates the treatment of the surrounding matrix by a single phase formulation and neglect of the leak-off term IV ). Hence, the deformation state of the linear-elastic rock matrix, embedding the fracture, is characterized by effective bulk K ef f and shear µ ef f modulus representing Gassmann's low frequency result [2,3]. In eq. (9), φ 0 denotes the initial porosity of the porous matrix, K s the (average) modulus of the compressible grains composing the matrix, K and µ the bulk and the shear modulus of the dry skeleton, and K f the bulk modulus of the fluid.

Numerical Model
The flow model requires input values for the elastic parameters of the matrix K and µ, the fluid bulk modulus K f , the initial fracture opening δ 0 , the equilibrium-normal stiffness parameter E Fr eq , the dimensionless contact parameter s 0 , and flow-boundary conditions for the intersection of the fracture plane with the borehole, as prescribed by the individual experimental protocols followed for the tests in the six intervals. The chosen parameters are listed in Table  4; while Freiberg gneiss is anisotropic (see [48]), for simplicity, we rely on representative isotropic material parameters. Freiberg gneiss appears homogeneous in the tested rock volume and we therefore employ a constant value of 4 for the parameter s 0 , determined from exploratory calculations. The chosen value represents a rather distinct difference between initial hydraulic and mechanical fracture apertures. The remaining model parameters, the initial hydraulic fracture aperture δ 0 and the equilibrium fracture normal stiffness E Fr eq along with the fracture length l Fr resulting from the discretization of the modelled domain, determine the effective hydraulic conductivity and the storage capacity of the tested fractures and the initial normal stress acting on the fracture surfaces. We seek optimized values for these by analysing the misfit between numerical pressure transients and observed pressure transients.
The assumption of a radial-symmetric fracture geometry and a linear-elastic response of the poro-elastic matrix reduce the total number of degrees of freedom (DoF) due to a reduction of the dimension and the neglect of pore-pressure effects in the surrounding domain, respectively. The reduction of DoF results in a high efficiency of the method; simulations of transients require just several minutes on a standard desktop PC with the used numerical discretization corresponding to around 20.000 DoF for the whole set of modelled fractures. 6

Quantification of misfit
Identifying "the" best numerical fit requires examination of the evolution of the misfit between observed transients p exp and the modeled transients p num . We quantify misfit by a normalized L 2 -error measure i.e., misfit is reduced to a single scalar value for each considered parameter combination. The variation of the misfit with combinations of the model parameters, i.e., the existence of global and/or local minima, is not known a priori, but is required to determine the quality and variance of best numerical fits. Iso-surfaces of misfit in the three-dimensional space of the model parameters {δ 0 , E Fr eq , l Fr } were gained from interpolation between the discrete values of actually performed calculations.

Strategy of parameter search
We separately investigated the misfit spaces for the pressure-plateau group and the positive-tangent group. The sensitivity of the model to its parameters was studied in a total of 1144 and 735 simulations for the pressure-plateau and the positive-tangent group, respectively. Numerical fits suitable for characterization of tested fractures should possess a normalized error of approximately e L2 ≤ 0.05 corresponding to absolute pressure deviations continuously below 0.2 MPa, subsuming the effect of flow-rate fluctuations due to irregularities of the pump, and intrinsic accuracy of sensors and the correction of flow rate for storage capacity of the injection system. The investigated ranges of the individual parameters (Table 3) were defined based on an exploratory analysis starting from educated guesses for the in-situ tests. This exploratory search indicated the existence of a potential misfit minimum below the defined error threshold whose location we then investigated further. For the subsequent analysis of the remaining sets of measurement data, knowledge about the existence of a local minimum motivated a strategy based on biased guesses to numerically determine parameter combinations resulting in the investigated error ranges corresponding to high quality fits. On average, this resulted in approximately 20-25 simulations per transient.

Results
The parameter study aimed at the identification of parameter combinations yielding good fit between field data and numerical simulations to understand the characteristics of fractures responsible for the two distinct groups of pressure transients. In a first step, we focused on one data set of each group, the transient pressure response M p a obtained from hydraulic tests at 51.6 m borehole depth, representing the pressure-plateau group, and data set M t a corresponding to tests conducted at a borehole depth of 24.6 m, representative for the pressure-tangent group before we used the gained knowledge about the existence and magnitude of a local error minimum to reduce the computational costs of the numerical fitting procedure of the remaining data sets by focusing on parameter combinations resulting in low error values.

Pressure-plateau group
The error surfaces of the pressure-plateau group possess an ellipsoid-like shape (Figure 4.1.1); the closed surface with accurate numerical solutions of e L2 ≤ 0.0275 indicates the existence of an error minimum. Examination of the error evolution with fracture length (A I to H I in Figure 4.1.1) shows consistency with the introduced iso-surface plot since the error reaches a minimum at a fracture length of 50.0 m. The error evolution is asymmetric around this minimum, it increases less steep for fractures with an increasing than for fractures with a decreasing length. The corresponding unique error minimum in the parameter space occurs for the parameter combination D I , which consists of the initial hydraulic fracture aperture δ 0 =36.0 µm, an equilibrium fracture normal stiffness parameter of E Fr eq = 5.6 MPa, and a fracture length of 50.0 m. The model shows higher sensitivity to the equilibrium-fracture normal stiffness parameter and the initial hydraulic aperture than to fracture length. Individual variations of the fracture stiffness and initial hydraulic aperture relative to the values obtained for the identified minimum (exemplified by D a I to D d I in Figure  4.1.1) result in pronounced under-and overestimation of the measured transients, respectively.

Positive-tangents group
For the positive-tangent group, the misfit surface identifying relevant parameter combinations with absolute differences to the measured data consistently below 0.2 MPa, defined by e L2 ≤ 0.055, consists of two connected ellipsoidal shapes with different orientations of their axes. A single potential minimum is indicated by the closed iso-surface with e L2 ≤ 0.03. The error evolution with fracture length (A II to G II in Figure 4.1.2) confirms the existence of a minimum for the parameter set B II , consisting of an equilibrium-fracture normal stiffness parameter E Fr eq = 2.8 MPa, an initial aperture of 42.0 µm, and a fracture length of l Fr = 4.75 m. Large misfits result when fracture length decreases below 4.75 m; however, misfit is less sensitive to variations in fracture length above this value. Variation of the equilibrium-fracture normal stiffness parameter E Fr eq and initial hydraulic aperture δ 0 relative to the parameter set B II (i.e., parameter sets B a II to B d II in Figure 4.1.2) reveals a higher sensitivity of the model to changes of the stiffness parameter than the initial hydraulic aperture.

Characteristic fracture properties
For each of the remaining four data sets, parameter sets resulting in fits with error values close to those obtained for the local minima in the two examples above could be determined ( Figure 4, Table 4). The error for data set N p c is exceptionally high compared to that of other sets when we do not neglect the first pumping step ( Table 4) that involves a delayed pressure increase (Figure 1). This interval looses water when isolated and has to be refilled after an extended shut-in period. The identified error minima correspond to distinctly different parameter ranges for the two groups.
The elongated enclosing hull of the determined parameter combinations for the pressure-plateau group visualizes the recognized scaling of fracture stiffness and initial fracture aperture with fracture length. For the positive-tangent group, we observed scaling of the fracture stiffness with fracture length, but no correlation between initial aperture and fracture length does exist (Table 4, Figure 4). Optimal parameters of the two groups occupy distinctly different volumes of the misfit space.

Discussion
Summarizing, throughout the numerical fitting of transient measurement data with pronounced pressure-plateaus a unique minimum in terms of a corresponding parameter combination could be identified in the range of investigated material parameters. The sensitivity analysis proofed continuously increasing errors for relative changes of parameters compared to the set D I which indicates that no further minima exist in experimentally justifiable limits of the parameters. The model exhibits a high sensitivity for the equilibrium fracture normal stiffness and the initial fracture aperture, whereas simulated pressures are fairly insensitive to changes in fracture length.   [52,50,51]. The fracture lengths of meter-scale derived for the positive-slope group of pressure transients are consistent with the spatial scale of the test volume and the dimensions of seismicity clouds observed during the corresponding stimulations. The decameter-scale fracture lengths modelled for the pressure-plateau group appear long at first glance. Yet, considering the shape of the misfit iso-surface of this group that documents an insensitivity of the model to changes in fracture length beyond a critical lower bound, fracture lengths barely exceeding 10 m cannot be excluded per-se. Furthermore, the model involves only a single fracture, neglecting leak-off into intersecting fracture systems and thus its application to data determines properties of an equivalent fracture potentially subsuming pre-existing fractures with a comparable or higher conductivity than that of the fracture intersecting the borehole.

Distribution of pressure along the fracture
Knowing the distribution of fluid pressure along the fracture is a crucial pre-requisite for substantial stimulation modeling. Since our hydro-mechanical model includes the entire fracture, we can use the determined parameter sets to investigate the pressure distribution in the fracture at any point during the step-rate tests. We focus on the pressure states at the end of each applied flow-rate step for data sets N p a and N t b that are representative for their corresponding group to examine whether the characteristics of the transients observed in the borehole bear information on the pressure distribution along the fracture.
For N p a , the representative of the pressure-plateau group, pressure gradients along the fracture are higher than for N t b of the positive-tangent group, for which the pressure profile is almost flat, i.e., the pressure in the fracture is equilibrated during every state of the pumping ( Figure 5). This significant difference in pressure distribution reflects the critical interrelation between local deformation and its consequences for local flow and storage. For the long fractures of the pressure-plateau group, the local deformation and thus permeability decrease with distance from the injection point, but storage of fluid becomes easier close to the borehole where, due to the increased fluid pressure the fracture is already less stiff than at its end. Pressure in the relatively short fractures of the positive-tangent group increases during each flow-rate step, as documented by the continuous increase in pressure in the borehole, and from step to step. The close to constant pressures document that pressure is not controlled by transport restrictions in them but by their storage capacity, which is limited owing to the direct effect of fracture length on fracture volume and also on geometrical fracture stiffness, as detailed in the next section.

Specific normal stiffness
The contact mechanics of the six investigated fractures is uniquely determined by the parameters constrained by the modeling. Evaluating the constitutive relation (8) with the found equilibrium-fracture normal stiffness parameter E Fr eq and the initial hydraulic fracture width δ 0 yields their opening and closure behavior when subjected to normal stresses deviating from the equilibrium stress ( Figure 6). The corresponding specific contact stiffnesses reflect the strong non-linearity of the constitutive relation; close to the equilibrium stress specific stiffness vary between 10 2 and 10 3 MPa/mm and thus fall well within the range of previous in-situ observations and laboratory studies [50,51,12]. 9 When extended towards fracture closing, i.e., when effective normal stresses exceed the equilibrium in-situ stresses, the stiffness of all tested fractures converge to a rather limited range of 2 · 10 3 MPa/mm to 4 · 10 3 MPa/mm. The equilibrium-fracture normal stiffness parameter E Fr eq defines the tension limit, at which the two fracture halves separate, and thus provides a constraint on the equilibrium stress that the fractures experience in-situ. Since the obtained equilibrium-fracture normal stiffness parameters of the pressure-plateau group are higher than the ones of the positive-tangent group, the predicted equilibrium normal stresses for fractures of the plateau group, ranging between 3.5 MPa and 5.6 MPa, are larger than the ones for fractures of the positive-tangent group, ranging from 2.1 MPa and 2.8 MPa. Magnitude and range of these predictions are consistent with the stress state inferred for the test volume at Reiche Zeche [48].
The contribution of the elastic medium, in which the fractures are embedded, to the stiffness of the entire system is conventionally addressed as geometrical stiffness [34,8]. We evaluated the balance between contact stiffness and geometrical stiffness by prescribing constant fluid pressures in the range of the experimental pressure levels with an increment of 0.5 MPa for the two equivalent fractures found for intervals 51.6 m (M p a ) and 40.6 m (M t b ), representing the two pressure-transient groups and constituting the upper and lower bounds of the parameter space of optimal fits in terms of fracture length and normal stiffness parameter, respectively.
The prescribed levels of fluid pressure lead to local deformations according to the constitutive relation (8) and associated local normal contact stresses σ Fr N , which we integrate over the fractures' lengths. Mechanical equilibrium across the fracture requires changes in fluid pressure and total normal stress to hold ∆p = ∆σ Tot N . Thus, the mismatch between the applied fluid pressure and the integrated normal contact stresses corresponds to the normal stress exerted on the fracture by the deformation of the surrounding material, here addressed as geometrical normal stress σ G N . The decomposition of the changes in total acting normal stress ∆σ Tot N = ∆σ Fr N + ∆σ G N (11) gives changes in the geometrical stress as The stress balance differs for the two investigated fractures and varies with fluid pressure for an individual fracture (Figure 7). For the long (50 m) fracture of the pressure-plateau group, force balance across the fracture is dominated by contact stresses, while the contribution of geometrical normal stress is significant for the short (3.7 m) fracture of the positive-tangent group, the more the higher the fluid pressure. The changing relative contributions result from the non-linearity of the normal-contact stress formulation introduced by eq. (8).
We transfer the findings for the normal-stress decomposition (Figure 7) to a corresponding decomposition of specific normal stiffness. Geometrical stiffness of the two investigated fractures is studied by conducting a numerical analysis of their opening behaviour under the assumption of negligible contact normal stresses, which results in p = σ G N . The discretized geometrical stiffness is then evaluated by means of an averaged aperture and the discretized normal stress change K G = ∆σ G N /∆δ G ; fracture contact normal stiffness is calculated by inserting the fitted parameters into the analytical derivative K Fr = ∂σ Fr N /∂δ Fr . The combined stiffness is numerically determined by K Com = ∆σ Tot N /∆δ Com . Since the relation between σ Tot N = p, σ G N , and σ Fr N is known from the stress decomposition (Figure 7), the resulting stiffness components can be transformed to express the stiffness of the combined model for a uniform fluid pressure.
The sum of transformed geometrical and contact normal stiffness agrees with the combined stiffness (Figure 8) lending support to the assumptions made regarding the transformation of different stress states to the acting effective normal stress. Geometrical stiffness is found to be negligible for data set N a p of the pressure-plateau group, for which the combined specific stiffness is well approximated assuming equivalence of fluid pressure and acting normal contact stresses (see Figure 6). In contrast, the combined specific stiffness for data set N t b , the representative of the positive-tangent group, is a superposition of both stiffness components converging towards the geometrical stiffness with increasing fluid pressure. In fact, the geometrical stiffness forms the lower bound for combined stiffness. The contribution of the geometrical stiffness to the total fracture stiffness is highest for fluid pressures above the identified asymptotic stress values of the contact model, i.e., at the onset of separation of the fracture halves.

Conclusions
We numerically modeled the opening characteristics of fractures during in-situ hydraulic tests using a hydromechanical flow model, here implemented for radial fractures with non-linear contact mechanics and without leakoff. Systematic variations of experimentally observed pressure transients lead us to distinguish two groups, one with continuously (positive-tangent) and step-wise rising (pressure-plateau) pressure response to stepwise increases in flow rate. Our comprehensive parameter study unveiled the effect of fracture length, equilibrium-normal stiffness parameter, and initial fracture aperture on the the misfit between observed and calculated pressure transients. The proposed hydro-mechanical coupling can explain the strikingly different pressure transients within experimental uncertainty and thus provides a perspective to the response of fractures to pumping operations alternative to the traditional pressurediffusion analyses, which relate the distinct pressure groups to differences in flow regime related to differences in the orientation of the fracture relative to the borehole The identified minima in mismatch between observed and calculated pressure transients correspond to different fracture properties for the two groups, though we noticed a close correlation between fracture length and fracture normal stiffness resulting in a specific mismatch for the positive-tangent group. Pressure plateaus are characteristic of relatively long and stiff fractures, while relatively short and compliant fractures lead to continuously increasing injection pressures. The pressure distribution along the fractures differs significantly for the two groups of pressure transients; pronounced non-linear pressure gradients develop in the long fractures of the pressure-plateau group during injection, while the pressure in the short fractures of the positive-tangent group remains close to the injection pressure along their entire length. Our observations on pressure distribution motivate to investigate the validity of the common practices of normal stress estimation from shut-in and jacking pressures.
Throughout the performed step-rate tests, fluid injection results predominantly in opening of fractures. Yet, the implemented non-linear constitutive relation for the contact mechanics of the fractures covers their opening and closing relative to the initial in-situ state associated with the in-situ stress state. Thus, evaluation of the constitutive relation with the determined fracture parameters allowed us to investigate the contributions of local contacts and overall fracture geometry to stress balance and thus bulk stiffness. With decreasing effective stress, the role of the contacts diminishes and total stiffness approaches the lower bound constituted by the geometrical stiffness.
The proposed hydro-mechanical model exhibits diminished sensitivity to fracture length when a flow-rate step results in a constant injection pressure. Thus, extending the pumping duration will not only help to discriminate between the alternatives of flow regime vs. hydro-mechanical effects, but may also reduce uncertainty of model parameters in case pressure ultimately deviates from the apparent early plateau. Future numerical work should explore different scenarios for the relation between mechanical and hydraulic apertures and its evolution with changes in effective normal stresses.           (Table 4). Bottom: Semi-logarithmic representation of the specific contact stiffness K Fr N as a function of the acting normal contact stresses. Positive values indicate compressive stresses (relative to the equilibrium stress). The legend applies to top and bottom, specifically dark blue lines represent data sets of the pressure-plateau group and grey lines that of the positive-tangent group.