Diagnosing Hydro-Mechanical Effects in Subsurface Fluid Flow Through Fractures

Hydro-mechanically induced transient changes in fracture volume elude an analysis of pressure and flow rate transients by conventional diffusion-based models. We used a previously developed fully coupled, inherently non-linear numerical simulation model to demonstrate that harmonic hydraulic excitation of fractures leads to systematic overtones in the response spectrum that can thus be used as a diagnostic criterion for hydro-mechanical interaction. The examination of response spectra, obtained from harmonic testing at four different field sites, for the occurrence of overtones confirmed their potential for the hydro-mechanical characterization of tested reservoirs. A non-dimensional analysis identified relative aperture change as the critical system parameter.


Introduction
The flow of fluids in fractures is important for a wide range of natural processes, such as fault reactivation and seismicity in general (Guglielmi et al., 2008(Guglielmi et al., , 2015a(Guglielmi et al., , 2015bDuboeuf et al., 2017;Cappa et al., 2018), for the recovering of fluid resources (Raghavan, 2004;Maliva, 2016) or heat (Baujard et al., 2017) from the subsurface, and for the subsurface storage of energy (Matos et al., 2019) or waste (Tsang et al., 2015). Whether the focus lies on fundamental science or on geo-engineering applications, there is little alternative but to perform tests in boreholes to constrain the relevant hydraulic properties on a spatial scale of tens to hundreds of meters. Traditionally, the inverse analysis of pressure and flow transients induced by pumping operations in boreholes (Muskat, 1938;Fetter, 2001) is based on the analytical evaluation of diffusion models for isotropic and, apart from discrete boundaries, homogeneous media (Matthews, 1961;Bourdet et al., 1989;Ferroud et al., 2019). However, theoretical work (Berryman, 1992) and laboratory experiments (Schepp & Renner, 2021) highlight the significance of hydro-mechanical coupling, i.e., the interrelation between fracture or more generally conduit deformation and fluid pressure evolution, for pressure and flow transients that poses a challenge to existing characterization strategies. Some field observations, such as reverse water-level fluctuations in distant monitoring boreholes (Vinci et al., 2014a), a phenomenon related to the Noordbergum effect known for poro-elastic media (Rodrigues, 1983;Kim & Parizek, 1997), demonstrate the substantial deviation of the spatial pressure distribution from expectations rooted in the conventional diffusion approach when conduits deform, and the neglect of conduit deformation in the analyses of pressure and flow transients potentially leads to largeThely inaccurate effective hydraulic properties (Vinci et al., 2014b). Recent hydro-mechanical tests of in-situ fractures provide direct evidence for fracture deformation (Dutler et al., 2019) and document the importance of fracture deformation for flow processes in fractured reservoirs (Schmidt et al., 2021a).
Capturing the contribution in non-local changes in fracture volume to storage of injected fluid volume in numerical modeling of flow processes in high-aspect ratio fractures requires to couple the fluid-flow field with the deformation state of the surrounding rock matrix. An evaluation of the Reynolds number for high-aspect ratio fractures motivates to assume creeping flow conditions and thus a Poiseuille-type flow in the fracture in common field and laboratory test situations (Louis, 1969;Witherspoon et al., 1980). The balance equations can then be consistently extended to flow in deformable fractures, resulting in a hybrid-dimensional flow model (Vinci et al., 2015), also referred to as the lubrication equation (Batchelor, 2000). The system of governing equations for deformable fractures embedded in rock-like materials is numerically stiff, especially for weakly compressible fluids, such as water (Adachi et al., 2007;Yew & Weng, 2014;Baotang et al., 2014). The discretized system of equations requires stabilized solution strategies to achieve numerical stability during equilibrium iterations for partitioned coupling schemes (Girault et al., 2015(Girault et al., , 2016Castelletto et al., 2015) and zero-thickness interface element formulations for monolithic schemes (Segura & Carol, 2008a, 2008bSchmidt & Steeb, 2019). Different mesh-based numerical frameworks have been proposed to provide solutions in the time domain Settgast et al., 2017;Hanowski & Sander, 2016;Schmidt & Steeb, 2019).
Field scale observations are needed to verify the theoretical and numerical hydro-mechanical approaches, and to investigate the applicability of parameters determined in small-scale laboratory experiments and to thereby provide the opportunity to assess the relevance of hydro-mechanical effects not only for pumping tests but also natural processes. In the light of previously proposed test methods (Rutqvist et al., 1997;Rutqvist, 2015;Dutler et al., 2019;Schmidt et al., 2021a;Cheng & Renner, 2022), a reliable assessment approach is particularly missing for tests in a single borehole that intend to characterize conduits in their current state rather than to induce conduit modifications, such as fracture growth. Periodic pumping tests (Rasmussen et al., 2003;Renner & Messar, 2006) promise a substantial sensitivity to hydro-mechanical effects since conduit deformation causes changes in permeability and storage capacity in every period, i.e., repeatedly, as previously demonstrated for open radial fractures (Vinci et al., 2015).
We address the need for the ability to assess the significance of hydro-mechanical phenomena by simulations of harmonic pressure excitation exerted on a fluid-filled, hydraulically open but mechanically closed radial fracture. We employ our previously presented hybrid approach that assumes creeping flow of a weakly compressible, viscous fluid in a planar fracture embedded in a deformable solid (Vinci et al., 2014a;Schmidt & Steeb, 2019). In the light of the numerical results, we introduce a nondimensional description that allows us to compare the influence of hydro-mechanical coupling on the response spectrum observed for the numerical model to that of transients recorded during a range of in-situ tests.

Governing Equations and Implementation
Our numerical approach employs an implicitly coupled, non-linear, time-domain formulation of a hybrid-dimensional model (Vinci et al., 2015;Schmidt & Steeb, 2019;Kolditz et al., 2021) including fluid pressure-induced changes of local permeability and fracture volume, overcoming the limitations of purely diffusion-based models (for details see Appendix 1). The governing equations account for non-local deformation with an instantaneous impact on the fracture volume, in contrast to the timedelayed change in fluid pressure caused by diffusion, and for local changes in fracture permeability due to changes in fracture aperture, either induced by fluid pressure variations along the fracture.
The absence of material characterizing a fracture underlies its hydraulic properties, permeability and specific storage capacity, i.e., void space in a solid is needed to store and transport fluid. In contrast, the mechanics of a fracture can be addressed by the abstract concept of an interface that affects the bulk deformability of the solid body embedding it. The deformability of the interface is typically addressed by defining a stiffness, with the dimension stress per displacement, i.e., an element property differing from the material perspective expressed by, e.g., Young's modulus for a solid. In reality, fracture stiffness and permeability are both closely related to measures of fracture aperture. This relation is probably at the core of the apparently universal scaling relation between fracture stiffness and fluid flow (Pyrak-Nolte & Nolte, 2016). Yet, hydraulic and mechanical responses to changes in fluid pressure along a fracture might in detail be governed by different aspects of the fracture geometry (e.g., Vogler et al. 2018).
The aperture changes of fractures are actually controlled by two ''stiffnesses'', one related to their overall geometry and the bulk elastic properties of the medium embedding the ''void'', and one representing the contact mechanics of the touching fracture surfaces. When a fracture is mechanically open, i.e., the two surfaces loose contact, e.g., due to the fluid pressure in it exceeding the normal stress on it, solely the geometrical stiffness determines their deformation; when they are mechanically closed, the contact stiffness tends to dominate their deformation (Schmidt et al., 2021b). The deformation characteristics of individual asperities forming the contact may be approximated by that of Hertzian contacts. Such nonlinear behavior particularly predicts ''relatively large'' aperture changes and associated changes in effective hydraulic properties of fractures for modest changes in normal stresses, when the effective normal stresses on the fracture are small. The chosen constitutive model for the local fracture deformation (see Eq. (9) in Appendix 3) accounts for the strong non-linearity between changes in normal stress and fracture aperture, characteristic for mechanically closed fractures in rocks (Bandis et al., 1983;Pyrak-Nolte & Morris, 2000), but might be replaced by alternative models to match experimental observations for a specific rock type and/or fracture. We adopt a sign convention, according to which compressive stresses are positive as are decreases in aperture. The (stress dependent) specific fracture stiffness C Fr , the change in normal stress r N required for a change in fracture aperture Dd ¼ d À d 0 relative to the initial (mean) aperture d 0 at zero stress, reads with the minimum fracture aperture d min reached for r N ! 1 and E Fr the fracture Young's modulus (Gens et al., 1990;Segura & Carol, 2008b). The specific fracture stiffness C Fr increases non-linearly with increasing compressive normal stresses r N . In our modeling, we assume the constitutive relation (1) describing the ''global'' behavior of a fracture on laboratory scale to also apply ''locally'', i.e., for every material point along the fracture. Furthermore, we replace the normal stress r N by effective normal stress r N ! r N À p, calculated as the unweighted difference between normal stress and the (local, time-dependent) fluid pressure.
A numerical simulation is performed in three steps; first, the fracture is subjected to a normal stress (10 MPa \r N \100 MPa); then a uniform fluid pressure (p eq ¼ 5 MPa) is applied in the fracture domain leading to the mechanical equilibrium state assumed to be characterized by an effective stress of r eq N ¼ r N À p eq (Table 1). Finally, sinusoidal pressure perturbations are imposed at the fracture's center with an amplitude of 4 MPa. The flow into a ''rigid fracture'' is simulated using a conventional pressurediffusion model to illustrate the effects of fracture deformation by a direct comparison of results from the two models.

Dimensionless Properties Quantifying the Hydro-Mechanical Interaction Throughout Harmonic Testing
A dimensionless formulation seems mandated for a meaningful comparison of results from different experiments because of the multitude of parameters describing the hydro-mechanical system. Several important system parameters, for example fracture permeability, vary spatially and temporally once fluid pressure is harmonically perturbed. Thus, we represent them by their values for equilibrated fluid pressure at a given normal stress, as denoted by a superscript ''eq''. Specifically, we introduce the characteristic time of pressure diffusion in a fracture which constitutes an approximate measure of the time it takes for the pressure perturbation to reach a characteristic diffusion distance l Fr (Weir, 1999), here corresponding to the fracture radius or length, g fR and b f viscosity and compressibility of the fluid, and k Fr;eq ¼ ðd eq Þ 2 =12 is the intrinsic permeability of the fracture in the equilibrated yet unperturbed state, characterized by a constant aperture d eq . Furthermore, we characterize a rock domain of height h and a Young's modulus E by a nominal domain stiffness C ¼ E=h. This conversion of the material property E to a domain property C addresses the need to quantify the relation between the deformability added due to the fracture and the intrinsic deformability of the body embedding it. The employed four dimensionless parameters are (a) the dimensionless characteristic diffusion time s ¼ t d =T, a normalization of the characteristic diffusion time t d (2) by the testing period T, (b) the dimensionless domain stiffness C ¼ C=C Fr;eq that normalizes the nominal domain stiffness C by the specific equilibrium fracture stiffness C Fr;eq , i.e., (1) evaluated for r eq N , where we normalize the domain radius r by the borehole radius r b and (d) the dimensionless aperture-change X ¼ p A s, for which we multiply the amplitude of the applied fluid pressure, p A , with the relative sensitivity of aperture changes to perturbations in effective normal stress s ¼ Ào lnðdÞ=o r N ¼ o lnðdÞ=o p ¼ 1=ðdC Fr Þ (relative fracture compliance).

Hydro-Mechanical Characteristics of the Numerically Investigated Fracture Domains
We employ experimental observations on single fractures embedded in granite rock samples (Pyrak-Nolte & Morris, 2000) to determine the parameters, fracture Young's modulus E Fr and the initial fracture aperture d 0 , of the proposed constitutive relation (Eqs. 1, 9 in Appendix 3). Specifically, we fit two sets of data for fracture deformation Dd and specific fracture stiffness C Fr as a function of normal stress r N , documented in Pyrak-Nolte et al. (1987) and Gale (1987) (Fig. 2), assuming d min ¼ 0, i.e., the fracture reaches a hydraulically closed state for infinitely high normal stress. One of the two domains, B 1 , corresponds to a fracture with a ''stiff'' response, whereas the other fracture, B 2 , is compliant for low normal stresses and converges to the stiff response of B 1 only for high normal stress.
To express their stress dependence (Fig. 2), we denote properties evaluated at minimum r eq N ¼ r min N À p eq ¼ 5 MPa and maximumr eq N ¼ r max N À p eq ¼ 95 MPa effective normal stress by h andĥ, respectively (see also Appendix 5). The equilibrium fracture apertures are moderately sensitive to applied normal stress for B 1 ( d B 1 =d B 1 % 2) but highly sensitive for B 2 ( d B 2 =d B 2 % 15); the associated stiffnesses stay within half an order of magnitude for B 1 but vary non-linearly over about two orders of magnitude for B 2 and reach the high values of B 1 only at the highest imposed normal stress values of 100 MPa. The domain stiffness C B i indicates a stiffer characteristic for the smaller sample B 1 , though the same Young's modulus E sat is used, because it scales with cumulative height of the rock cylinders h B i . The normalized radius K indicates distinctly different fracture lengths for the two numerically investigated domains (K B 2 =K B 1 ¼ 10).
The stress state-dependent, dimensionless domain quantities of B 1 evolve more with stress than those of B 2 over the range of investigated effective stresses. The dimensionless domain stiffness is insensitive to changes in stress for B 1 but decreases significantly for B 2 with increasing normal stress. The fracture stiffness persistently exceeds the nominal stiffness of domain B 1 . For domain B 2 , the relation between fracture stiffness and nominal domain stiffness switches with increasing normal stress; the fracture deforms more than the rock body alone for the minimum applied normal stress but significantly less at the maximum normal stress. The relative aperture sensitivity s, determined from initial fracture aperture and fracture Young's modulus, evolves distinctly different with increasing normal stress for the two numerical domains (Fig. 5c, d); s is low and constant for B 1 , while, for domain B 2 , it is high at low normal stress and decreases to the insensitive response of domain B 1 at high normal stress.
The diffusion time (2)  of the two considered fracture scenarios are of comparable magnitude at the minimum effective equilibrium normal stress ( t d B 2 = t d B 1 ¼ 6), they differ by more than two orders of magnitude at the maximum stress state (t d B 2 =t d B 1 ¼ 360). At all stress conditions, the characteristic diffusion times are, however, significantly smaller than the selected period T of the pressure perturbation (t d =T $ 10 À7 ). The similarity of characteristic times at minimum stress and its relation to the excitation period mean that, if deformation were absent, (a) pressure diffusion would take place on a similar time scale in the two fractures and (b) pressure would be almost equilibrated along the fractures and thus close to the excitation pressure at all times. With the selection of the two samples, we thus ensure that differences in the hydraulic response of the two fractures to the periodic excitation are dominated by deformation characteristics and that the potential effect of fracture deformation on the flow process is maximized.

Figure 1
The numerical boundary conditions of the axisymmetric model comprise lateral stress r M on the outer sample surface related to radial pressure p r , axial or normal stress r N , normal stress on the injection borehole wall r pðtÞ , and the Dirichlet boundary condition associated with the fluid pressure pðtÞ at the intersection between fracture and injection borehole

Results
Evaluation of the flow-rate (q) transients in the time domain reveals distinctly different responses of the two modeled domains to the harmonic pressure (p) excitation (Fig. 3). For the ''stiff'' fracture B 1 , flow rate is close to harmonic and exhibits extrema and zero-crossings opposite to the pressure irrespective of normal stress, corresponding to a circular p À q hysteresis loop. In contrast, the flow response is visibly non-harmonic for the ''compliant'' domain B 2 , with extrema that are the closer to the pressure maximum the lower the normal stress and correspondingly the fracture's stiffness. The injection rate into the fracture peaks close before the maximum in injection pressure when the maximum fracture aperture and minimum stiffness are reached, and likewise the production rate from the fracture is the largest right after the peak in pressure.
The deviations of the flow responses from a harmonic signal observed for the compliant domain B 2 correspond to systematic overtones in the spectra at integer multiples of the excitation frequency, i.e., harmonics. The overtones exhibit decreasing amplitudes with increasing frequency at a given normal stress and with increasing normal stress at a given frequency (Fig. 5). The phases of the overtones alternate; odd harmonics (first, third,...) have a phase shift of p or 180 , i.e., are inverted, compared to the flow rate at the excited frequency, while even ones (second, fourth,...) share the phase of the flow rate signal.
For the excited period, the spectral parameters, amplitude ratio and phase shift between flow rate and pressure, exhibit systematic trends with normal stress that differ for the two samples (Fig. 4). Amplitude ratios, measures of the fracture's injectivity (e.g., Jiménez Martínez & Renner 2021), decrease with increasing normal stress but much less so for the less deformable sample B 1 . The deformable sample B 2 starts out with an order of magnitude higher flow rate than B 1 for the imposed pressure oscillation at the lowest normal stress (see also Fig. 5a, b), but its injectivity falls below that of B 1 by about an order of magnitude at the highest imposed normal stress. The amplitude ratio-normal stress relations (Fig. 4) exhibit an astonishing qualitative similarity to the relative sensitivity of aperture changes to perturbations in effective normal (Fig. 2c).
Phase shift values are all close to-but actually slightly exceed-T=4 (Fig. 4b), the limit reached for large periods for purely diffusional radial flow with a no-flow boundary (Renner & Messar, 2006). The phase shift is insensitive to normal stress for the stiff sample B 1 , but decreases with increasing normal stress for the deformable sample B 2 . The variations in phase with normal stress are minute compared to those of amplitude ratio; the proximity to T=4 is consistent with the small dimensionless diffusion times of the simulations. The scatter of phase shift values around the trends probably indicates the limit in resolution of this spectral parameter associated with the numerical simulations.

Discussion
The occurrence of systematic overtones, here actually harmonics, in the flow rate's response to the harmonic pressure boundary condition at the fracture's center may in principle result from any -or a combination-of (i) the non-constant parameter, aperture, and (ii) the related coupling term in the governing partial differential equation (6), and (iii) the material non-linearity represented by the constitutive relation for the deformation behavior of the fracture (9) employed in our mathematical model for periodic radial fluid flow in a deformable fracture. As much as a rigorous analysis of the options for a mathematical model reduction, similar to approximating (5) by (6) based on a dimensional analysis, would be beneficial for an identification of critical model features and for a reduction in computational efforts, it goes beyond the scope of the present study and we restrict here to conceptual consideration in the light of the diverse generation of harmonics.
Harmonics are well know as a resonance feature when the spatial restrictions to progressive waves match their eigenfrequencies and thereby lead to standing waves corresponding to free oscillations, e.g., the modes of musical instruments such as organ pipes. Non-constant material parameters can also lead to harmonics for spatially unrestricted progressive waves, as well known from non-linear acoustics; the solutions of the associated uncoupled hyperbolic equations with non-constant parameters are analytically well investigated and constitute series of harmonics (e.g., eq. 265 in Blackstock et al. 1998;Hargrove 1960;Garrett 2020). Likewise, the solutions of uncoupled diffusion equations, the subset of hyperbolic equations reflecting fully dissipative behavior excluding wave propagation, for forced oscillations constitute series of harmonics when the diffusivity is a non-constant parameter (Shepherd & Wiltshire, 1994).
The only intrinsic ''material'' non-linearity involved in our numerical model, the stress-dependent fracture stiffness (9), nominally affects the system's response to the imposed fluid-pressure perturbations via the coupling term and the non-constant parameter in the governing equation. Our modeling results do not allow us to assess whether a linear constitutive relation, i.e., a constant fracture stiffness, would still give rise to harmonics. The investigated stiffnesses involve a bias between absolute stiffness and non-constant stiffness (Fig. 2c). The qualitative similarity between the relation of amplitude ratio to effective normal stress (Fig. 4a)-that in field tests can be investigated by changing the mean flow rate and thus the mean fluid pressure-and the stress dependence of the fracture's relative compliance (Fig. 2c) hints at a role of the constitutive non-linearity for the system's behavior and at a related potential to constrain details of fracture mechanics in field tests.
Coupling of differential equations does not per-se yield non-linearity, a prominent counter-example in this context given by the set of coupled differential equations representing Biot's poro-elasticity theory. Yet, the type of coupling differs fundamentally between the investigated fracture problem and Biot poro-elasticity; coupling in the latter is on the level of a material point while in the former the coupling arises due to a geometrical boundary condition. Furthermore, the changes in fracture aperture due to changes in fluid pressure and/or normal stresses reflect anything but small ''strains''. Even if general statements regarding the relative size of the differential terms in (6) were gained and thus conditions justifying the neglect of the coupling term were in principle identified, one would face the challenge of accounting for the non-local effect of pressure distribution on fracture aperture by a non-constant parameter in lieu of the coupling term. Only for mechanically closed fractures experiencing fluid pressures much smaller than the acting normal stress, i.e., far away from their opening, the non-local relation expressed by the geometrical stiffness becomes subordinate compared to the solely locally operative contact stiffness, and thus a description with a nonconstant parameter but without a coupling term seems viable in principle. In any case, the non-constant parameter and the coupling term in the governing Eq. (6) are expressions of the same physical phenomenon, hydro-mechanical coupling, and thus the occurrence of harmonics can be used to assess the significance of hydro-mechanical coupling for a specific flow process.

Diagnostic Potential of Overtones for Hydro-Mechanical Interaction During Harmonic Excitation
Our numerical simulations indicate that overtones in the frequency domain of the flow response permit inferring the relevance of hydro-mechanical interaction and thus assessing the validity of a diffusion approach. The hydro-mechanical response is expected to be controlled by (a) the ''extent'' of the perturbation related to the excitation characterized by the dimensionless characteristic diffusion time t d (2) as a measure for the hydraulic penetration depth, and (b) the intrinsic ''responsiveness'' of the fracture in terms of the dimensionless aperture-change X. The relation between changes in effective aperture d and effective normal stress controls the hydro-mechanical responsiveness of the fracture and thus will be examined in detail.
In the numerical model, we assume the local mechanical equilibrium of a fluid-filled fracture under normal stress to be constrained by the stress balance i.e., the acting normal stresses r N are in equilibrium with the sum of fluid pressure p, normal contact stresses r Fr N , and stress r St related to the deformation of the surrounding bulk material, typically characterized by a structural or geometrical fracture stiffness (Murdoch & Germanovich, 2006;Schmidt et al., 2021b). We, however, neglected this geometrical stiffness, because the applied pressure amplitudes remained well below the acting normal stresses and therefore the limit of fracture-surface separation. Hence, variations in effective normal stress r N À p are exclusively compensated by contact normal stresses r Fr N and thus the numerical simulations correspond to an upper bound for the contact-related deformation of a tested fracture.
The amplitudes of overtones observed for the two domains correlate with their different dimensionless aperture-change values X (Fig. 5); overtones in the frequency domain are prominent for high dimensionless aperture-change values. Thus, we identify the amplitude of overtones as a measure for changes in flow characteristics of a tested fracture due to hydromechanical interaction.

Interpretation of Overtones in the Response Spectra of Previously Reported Field Data
We investigated the records of field studies that employed harmonic pumping (see Appendix 2 for site and test details), in the Reiche Zeche undergroundresearch laboratory (RZ), Freiberg, Germany (Renner and STIMTEC-team, 2021), at the Terrieu karstic field (TER) near Montpellier, France (Fischer et al., 2018), and at the Altona Flat Rock Site (AFR) in Clinton County, NY, USA (Guiltinan & Becker, 2015). We complement the three studies involving fractured rock masses with a study conducted at the Boise Hydrogeophysical Research Site in Boise (BHR), ID, USA (Rabinovich et al., 2015) that involved tests in a gravel aquifer, because to investigate whether our fracture-deformation based concept can be generalized to conduit deformation. Since it is technically simpler to control flow rate than pressure in field tests, for them, the roles of imposed boundary condition and response are Vol. 180, (2023) Hydro-Mechanical Effects in Fractures 2849 inverted for pressure and flow rate compared to the numerical simulations. The amplitude spectra of the recorded pressure transients can be classified into three groups (Fig. 6): overtones occur (i) at local maxima at multiples of the excitation frequency (RZ-URL, BHR), (ii) at maxima at uneven multiples of the excitation frequency (AFR), and (iii) unsystematically (TER). We compare the overtone observations with estimates of the introduced dimensionless properties to assess the role of hydro-mechanical interaction for their occurrence. Qualitatively, we expect the fields of strong and weak hydro-mechanical interaction to correspond to combinations of low and high specific fracture stiffness with low and high diffusion times. Hydromechanical coupling is expected to increase with increasing fracture deformability and applied pressure amplitudes corresponding to increasing dimensionless aperture changes and decreasing pressure diffusion times corresponding to increasing hydraulic penetration depths.
The placement of the different field studies in the ''regime map'' and the occurrence of overtones is consistent with expectations and shows a significantly higher sensitivity to the dimensionless aperture change than to the characteristic pressure diffusion time: • The dimensionless properties indicate a high dimensionless aperture-change and a moderate penetration depth for RZ. The similarity between dimensionless parameters and overtone characteristics for RZ with B2 strongly suggest that hydromechanical coupling is the cause for the non-linear response. • The AFR experiments are characterized by a low to moderate dimensionless aperture-change and a high penetration depth in agreement with the amplitude of observed overtones at uneven multiples of the excitation frequency. The absence of even multiples of the excitation frequency remains unclear and might be related to additional effects, such as non-laminar flow, not considered by the present model. • The combination of low values of dimensionless aperture change and a moderate penetration depth indicating negligible hydro-mechanical interactions is in agreement with the inconsistent occurrence of local maxima in the amplitude spectra of the TER response data. Frequencies exceeding the nominally excited one seem related to the insufficient time resolution of the data. • The amplitudes of observed overtones in the response spectrum of the BHR exceed that of overtones of the three field studies involving fractured rocks indicating significant hydro-mechanical interaction. The pronounced overtones correlate well with the high dimensionless porosity change estimated for granular media (see 4), in lieu of a dimensionless aperture change.

Implications
Our simulations were performed for a range in normal stresses mimicking the conditions, under which the parameters of the constitutive relation for fracture stiffness (1) were experimentally constrained. Yet, in field tests, the far-field stresses are fixed, and thus the effective normal stress on fractures can only be affected by changing the mean fluid pressure along a fracture by changing the mean flow rate of the periodic oscillation (see for example Cheng and Renner (2018); Jiménez Martínez and Renner (2021)). For example, our numerical results (Fig. 4a) predict a non-linear increase in injectivity with increasing mean fluid pressure for a given test period.
An account of hydro-mechanics will improve the inverse analysis of experimental field and laboratory data by separating hydro-mechanical phenomena from undesired measurement artifacts. Improved estimates of hydro-mechanical properties of a tested region then bear a significant potential to optimize, for example, the design of energy-or waste-storage and/or geothermal systems. An improved (quantitative) understanding of hydro-mechanical properties of fractures and fractured rocks on field scale is also essential for a substantiation of analyses of seismic (Brown et al., 2005) and volcanic (Aki & Koyanagi, 1981) tremor, and tidal signals in boreholes (Zeumann et al., 2009).

Conclusions
Our numerical modeling suggests that the nonlinearity associated with conduit deformation leads to overtones in the response spectrum of harmonic testing for ''common'' normal stiffness values of fractures. Overtones observed in field tests correlate with the sensitivity of hydraulically relevant conduit properties, e.g., fracture aperture or total porosity, to perturbations in fluid pressure; the correlation indicates that the deformation characteristics of conduits control the degree of non-linearity of fracture flow. Hence, overtones observed during harmonic hydraulic testing should not per-se be discarded as measurement imperfections, but are an indicator of hydro-mechanical interaction with practical relevance and bear the potential to characterize a fracture's stress state and mechanical properties, such as the specific normal stiffness. Overtones are related to the non-linearity associated with conduit deformation and focusing on them thus eliminates the ambiguity in interpretation of pressure transients regarding the role of hydro-mechanical coupling and subsurface heterogeneity.

Funding
Open Access funding enabled and organized by Projekt DEAL.Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommons.org/licenses/by/4.0/.

Data Archival
We used the modular toolbox DUNE (https://www. dune-project.org/), a free software licensed under the GPL (version 2) for solving partial differential equations (PDEs) with grid-based methods, to perform the numerical simulations. The Zenodorepository https://doi.org/10.5281/zenodo.6962627 contains our DUNE implementation of the hybrid dimensional interface elements, the files for compilation, the input and mesh files as well as the python script that generates the characteristic, dimensionless properties. The analyzed field data are compiled at https://doi.org/10.5281/zenodo.7295318.

Appendix A: Numerical Modelling
We split our numerical modelling domain into the fracture and the surrounding rock. While the rock matrix is modelled within a 3D-continuum setting, the deformable, fluid-saturated fracture is considered as a 2D-domain. The coupling of the deformable fluid-saturated fracture and the linear-elastic rock matrix is realized through compatibility conditions for the fracture aperture and the matrix deformation and through equilibrium conditions for the fluid pressure and the corresponding stress component of the rock matrix (Schmidt & Steeb, 2019). Different coupling schemes have been investigated (Schmidt & Steeb, 2019); in the case of mechanically-closed but hydraulically-open fractures, contact normal stresses and related constitutive assumptions have to been taken into account (Schmidt et al., 2021b).

A.1 Governing Equations for Deformable Fractures
Details of the derivation of the governing equations for fluid flow in a deformable fracture can be found in (Vinci, 2014). The approach incorporates the geometrical assumption of parallel plates for highaspect ratio fractures and assumes creeping flow conditions in the derivation of the balance equations. Evaluation of the balance of momentum for weakly compressible, viscous fluids in a hydraulic transmissive fracture results in a Poiseuille-type formulation (Vinci et al., 2014a;Schmidt & Steeb, 2019). Then, the relative fluid velocity w f is proportional to the pressure gradient grad p. Here, we treat a planar circular fracture such that the cubic law (Witherspoon et al., 1980) is evaluated at each material point Pðr; tÞ along the plane according to w f ¼ À dðr; tÞ 2 12 g fR op or e r ¼: À k Fr ðr; tÞ g fR where dðr; tÞ denotes the space-and time-dependent fracture aperture, k Fr ðr; tÞ ¼ d 2 ðr; tÞ=12 the spaceand time-dependent permeability, re r the position vector along the fracture, t time, and g fR the dynamic viscosity of the fluid. Using (4) in the evaluation of the balance of mass for the fluid and considering a mean velocity averaged over the fracture aperture results in a hybrid dimensional formulation (Vinci et al., 2014a(Vinci et al., , 2014b(Vinci et al., , 2015 ð5Þ where b f is the fluid compressibility and q lk accounts for leak-off, despite its facile connotation here meant to encompass all fluid exchange between the fracture and the surrounding rock, here, as commonly done, addressed as matrix. The hybrid character of the formulation lies in the reduction of the pressure field to a radial (in plane) dependence, while allowing for changes in fracture aperture (out of plane). Dimensional analysis of the governing Eq. (5) motivates to neglect the quadratic term ( II) and the convection term (III) for problems involving highaspect ratio fractures (Vinci et al., 2014a). Then, the hybrid dimensional description comprises a transient (I), a diffusion (IV), a coupling (V), and a leak-off (VI) part. In this work, we furthermore focus on the specific case of conductive fractures in poorly permeable rocks and thus leak-off driven dissipation is neglected, too, leading to the implemented approximate governing equation Fracture deformation modifies the fracture aperture dðr; tÞ and thus contributes to the fracture-flow solution via the diffusion term IV) and the volumetric coupling term V).

A.2 Governing Equations for Deformable Rock Matrix
The governing balance of linear momentum is given by Cauchy's equation where r s s denotes the solid stress tensor. In general, evolution of the deformation of the rock embedding the fracture and the pore pressure in it had to be analyzed relying on a poro-elastic theory, such as Biot's theory (Biot, 1941;Schmidt & Steeb, 2019). However, the neglect of leak-off driven dissipation in (6) allows for a single phase, linear elastic formulation for the surrounding rock (Hooke's law) where e s is the strain tensor of the solid, and volðe s Þ and devðe s Þ its volumetric and deviatoric part. In Eq. (8), we introduced the shear modulus G G sat and the saturated bulk modulus K sat of the poroelastic medium, here derived from Gassmann's lowfrequency result (Gassmann, 1951) or (Section 6.3 Mavko et al. 2009).

A.3 Model Geometry and Constitutive Assumption for Fracture
We envision a harmonic excitation in fluid pressure at the center of a circular planar interface, the fracture, between two cylindrical domains (Fig. 1). We investigate fractures in their post setting state (Bandis et al., 1983;Barton et al., 1985), i.e., fracture deformation is considered to be reversible for the investigated effective stress states. We assume a uniform initial fracture aperture and homogeneous properties of the solid cylinders to allow for a two dimensional, radial symmetric formulation of the governing Eqs. (5) and (8). The modeled fractures do not have physical ends, thus their geometrical stiffness is not considered in this study (Murdoch & Germanovich, 2006;Schmidt et al., 2021b); an assumption valid for fluid pressure amplitudes below the acting normal stresses, i.e., much smaller then the fluid pressures required for fracture-surface separation (Schmidt et al., 2021b). Yet, we imposed a noflow boundary at the end of the fractures.
The opposing fracture surfaces are locally in contact for a finite stress normal to the fracture surface, r N (Greenwood et al., 1966;Cook, 1992), leading to the contact stiffness. Shear deformations are insignificant when, as is typically the case in Vol. 180, (2023) Hydro-Mechanical Effects in Fractures 2853 hydraulic tests, the fluid pressure remains well below the jacking pressure, i.e., the pressure causing fracture surface separation, and we therefore neglect the contribution of shear stress induced fracture dilation (Schmidt et al., 2021b, a). We describe the contact mechanics on the continuum scale by an empirical, hyperbolic evolution law (Goodman, 1976;Bandis et al., 1983) Characteristic diffusion time at r eq N 2:7 Á 10 À8 1:6 Á 10 À7 (-) s Bi Characteristic diffusion time atr eq N 1:0 Á 10 À7 3:6 Á 10 À5 (-) 2854 P. Schmidt et al. Pure Appl. Geophys. where the changes in effective fracture aperture Dd relative to the initial aperture d 0 at zero stress are incorporated in the effective fracture aperture d ¼ d 0 À Dd, d min is the minimum effective fracture aperture reached for r N ! 1 and E Fr the fracture Young's modulus (Gens et al., 1990;Segura & Carol, 2008b).

A.4 Implementation
The non-linear system of partial differential Eqs. (6) and (8) is solved by an extended Dune-PDELab finite element implementation (Bastian et al., 2010). Discretization of the hybrid-dimensional flow model is realized by double-node zero-thickness elements, and the set of equations formed by monolithic assembly is solved by a direct solver in parallel. Averaging of balance of mass and momentum in the fracture domain is numerically realized by averaging the fluid pressure of facing fracture surfaces at aligning nodes of the implemented zerothickness interface elements (Segura & Carol, 2008a, 2008bSchmidt & Steeb, 2019). Numerical integration of the fracture-flow formulation is then performed on auxiliary lower dimensional elements. The auxiliary elements are used for integration purposes only, do not introduce additional degrees of freedom and therefore contribute to an increase in computational efficiency (Segura & Carol, 2008a;Schmidt & Steeb, 2019). Fracture aperture dðx; tÞ is handled on integration-point level accounting for the deformation of the zero-thickness elements. The framework uses a Newton-Raphson-based strategy to reach quadratic convergence of the equilibrium iterations.

B.1 Reiche Zeche Underground Research Laboratory
At the Reiche Zeche mine (RZ), harmonic testing was performed with a double-packer probe (Solexperts GmbH, Bochum, Germany) installed at a distance of 40.6 m from the head of 63 m long borehole BH10 drilled dipping downwards 15 to NNE from a tunnel about 130 m below the surface. The interval enclosed between the two packers had a length of 0.7 m; the injection system had a storage capacity of 1 Â 10 À12 m 3 /Pa. Stimulation preceding the harmonic testing lead to a pair of axial fracture traces at the borehole wall, complementing a preexisting circumferential fracture trace, as revealed by televiewer-logging and impression-packer testing (Adero, 2020). The data used here represent a sinusoidal variation in flow rate with a period of 65 s and an amplitude of 1.25 l/min around a constant flow rate of 3.85 l/min for a total of 5 periods. Pressure in the borehole exhibited a periodic response with an amplitude of 0.3 MPa and a mean value of 4.9 MPa.
We estimate fracture properties from results of a numerical hydro-mechanical characterization of steprate tests performed at six fractured intervals of the same borehole at different depths (Schmidt et al., 2021b). Apertures are approximated to range between 120 and 220 lm under an estimated effective normal stress of 1.1-1.6 MPa, and potential effective fracture lengths are expected to be on the order of 10 m. The fracture geometry is characterized by an estimated dimensionless domain length K RZ of 100 and is slightly larger than the fracture domain investigated for domain B 2 . The dimensionless characteristic diffusion time s RZ ranges from 2:0 Á 10 À4 to 5:1 Á 10 À4 indicating slightly slower pressure diffusion than for the numerical domains but still rapid compared to the excitation period. The dimensionless domain stiffness of 16-42 indicates a compliant fracture response, and a relative aperture-change of Classification of field tests and the performed numerical modeling regarding hydro-mechanical interaction (shading from gray to blue) based on the amplitudes of overtones in the response spectrum in relation to their dimensionless characteristic diffusion time, a measure for the expected hydraulic penetration depth, and dimensionless aperture (or porosity) change Characteristic occurrence of multiples of the excitation frequency similar to those observed for the numerical domain B 2 are present in the spectrum of the measurement data (Fig. 6 a).

B.2 Terrieu Karstic Field
At the Terrieu karstic field (TER), the harmonic pumping tested a highly conductive limestone aquifer with conduits of an estimated effective hydraulic aperture of 111-350 lm. From complementary experimental investigations, it is known that apertures can locally reach 20-50 cm (Fischer et al., 2018) and that the connected hydraulic pathways span a network extending over a total area of approximately 2500 m 2 . Monitoring wells are positioned at distances between 8 and 24 m to the pumping borehole. The dimensionless domain length of 1000 is significantly larger than for the numerical studies and the dimensionless characteristic diffusion time of 3:4 Á 10 À4 to 3:4 Á 10 À3 indicates slower pressure diffusion processes than for numerical domains B 1 and B 2 . The characteristic diffusion time is three to four orders of magnitude smaller than the applied period. We envision the voids in the karstic subsurface to differ from fractures in crystalline rocks so that changes in conduit volume are purely induced by deformation of the surrounding bulk material, i.e., we assume that changes in aperture are directly related to the domain stiffness. Overtones in the recorded data (Fig. 6b) do occur at apparently arbitrary frequencies rather than at multiples of the excitation frequency.

B.3 Altona Flat Rock Site
The harmonic pumping tests at the Altona Flat Rock Site (AFR) in Clinton County were performed near the surface (at 7.6 m depth) on a single, hydraulically dominant bedding plane fracture with a length of approximately 10 m, embedded in a lowpermeability sandstone (Guiltinan & Becker, 2015). With a dimensionless domain length of 200, the tested region is moderately larger than the ones numerically investigated. For characterization purposes, four monitoring boreholes were installed in a distance of 5 m to the excited borehole (radius 0.05 m) intersecting the fracture, too. Effective hydraulic apertures between 717 and 951 lm were estimated based on the measured fracture transmissivity (Guiltinan & Becker, 2015). The estimated characteristic diffusion time of 2:3 Á 10 À6 to 1:6 Á 10 À5 is comparable to the lower limit investigated for the numerical domains. Considering that pumping was performed near the surface, we expect a compliant fracture response based on fracture stiffnessdepth correlations previously reported for fractured granite reservoirs (Fransson, 2009;Rutqvist, 2015). The resulting dimensionless aperture-change of 0.01-0.1 indicates moderate hydro-mechanical interaction comparable to the lower range of numerical domain B 2 . The pressure response to harmonic flow excitation exhibits a dominant frequency, but shows additional amplitudes at uneven multiples of the excitation frequency (Fig. 6c).

B.4 Boise Hydrogeophysical Research Site
Harmonic tests at the Boise Hydrogeophysical Research Site (BHR) aimed to characterize a gravel aquifer, located close to the surface. The initial porosity / of the gravel aquifer ranges between 0.17 and 0.24 (Rabinovich et al., 2015). At the associated low normal stress conditions, single contacts between cobbles are expected to react compliantly to changes in effective stress induced by the pumping. The dimensionless domain length of 1000 translates to a larger tested region than numerically investigated. The high dimensionless characteristic diffusion times range between 0.6 and 2.1. We rely on the Kozeny-Carman equation, which correlates permeability and porosity for granular media (Carman, 1956), to infer deformation-induced hydraulic changes of the investigated region. Similar to the dimensionless aperturechange, we introduce the dimensionless porosity change X Ã ¼ p A o lnð/Þ=o p to quantify the hydromechanical interaction expected throughout a testing period, where we assumed that changes in effective stress are exclusively related to changes in fluid pressure. We estimated a dimensionless porosity change of 0.13-0.19 based on investigated density changes of gravel samples under small perturbations Vol. 180, (2023) Hydro-Mechanical Effects in Fractures 2857 of compression forces (Rücknagel et al., 2013). The occurrence of overtones in the spectra of the pressure response (Fig. 6d) show close similarities to that found for the responses of numerical domain B 2 (Fig. 5b). The authors stated that the occurrence of overtones scaled with the applied excitation period where longer periods resulted in higher pressure amplitudes, but did not provide any further explanation (Rabinovich et al., 2015).
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.