Exploratory use of periodic pumping tests for hydraulic characterization of faults

Exploratory use of periodic pumping tests for hydraulic characterization of faults SUMMARY Periodic pumping tests were conducted using a double-packer probe placed at four different depth levels in borehole GDP-1 at Grimselpass, Central Swiss Alps, penetrating a hydrothermally active fault. The tests had the general objective to explore the potential of periodic testing for hydraulic characterization of faults, representing inherently complex heterogeneous hydraulic features that pose problems for conventional approaches. Site selection reflects the specific question regarding the value of this test type for quality control of hydraulic stimulations of potential geothermal reservoirs. The performed evaluation of amplitude ratio and phase shift between pressure and flow rate in the pumping interval employed analytical solutions for various flow regimes. In addition to the previously presented 1-D and radial-flow models, we extended the one for radial flow in a system of concentric shells with varying hydraulic properties and newly developed one for bilinear flow. In addition to these injectivity analyses, we pursued a vertical-interference analysis resting on observed amplitude ratio and phase shift between the periodic pressure signals above or below packers and in the interval by numerical modeling of the non-radial-flow situation. When relying on the same model the order of magnitude of transmissivity values derived from the analyses of periodic tests agrees with that gained from conventional hydraulic tests. The field campaign confirmed several advantages of the periodic testing, for example, reduced constraints on testing time relative to conventional tests since a periodic signal can easily be separated from changing background pressure by detrending and Fourier transformation. The discrepancies between aspects of the results from the periodic tests and the predictions of the considered simplified models indicate a hydraulically complex subsurface at the drill site that exhibits also hydromechanical features in accord with structural information gained from logging. The exploratory modeling of vertical injectivity shows its potential for analysing hydraulic anisotropy. Yet, more comprehensive modeling will be required to take full advantage of all the pressure records typically acquired when using a double-packer probe for periodic tests. Fracture and flow, Hydrothermal systems, Europe, Fourier analysis, Numerical modelling, Mechanics, theory, and modelling 1 INTRODUCTION Analyses of pressure and/or flow transients associated with pumping operations in wells, that is, injection or production, constitute the primary means for hydraulic characterization of the subsurface (e.g. Fetter 2001). Effective hydraulic properties are determined relying on models considering specific flow geometries and boundary conditions. Specific aspects of observations may guide the selection of one out of a number of alternative models (e.g. Matthews & Russell 1967; Bourdet et al.1989). Periodic pumping tests involve consecutive periods of injection and/or production resulting in alternating flow and pressure and promise several practical and analytical benefits. They may be seen as either an alternative or a complement to conventional testing for constraining effective hydraulic properties (e.g. Rasmussen et al.2003; Renner & Messar 2006), diagnosing the flow regime (e.g. Hollaender et al.2002) and providing valuable information on subsurface heterogeneity (Renner & Messar 2006; Ahn & Horne 2010; Fokker et al.2012; Cardiff et al.2013). The periodic signals can be easily separated from noise (Renner & Messar 2006; Bakhos et al.2014) or background variations associated with long-term operation (e.g. Guiltinan & Becker 2015). In addition, numerical modeling of periodic tests in the frequency domain is much faster than that of typical transients in the time domain (Cardiff et al.2013; Fokker et al.2013). As is the case for conventional testing, effective (or equivalent) hydraulic parameters may be derived from periodic pressure and flow-rate observations by relying on a simple diffusion equation for a homogeneous and isotropic medium in the first step of analyses. Yet, their variation with pumping period can be attributed to—among other possible reasons—spatial heterogeneity (e.g. Renner & Messar 2006; Rabinovich et al.2015). Renner & Messar (2006) performed comprehensive field tests and found that diffusivity values estimated from the analysis of flow-rate and pressure responses at the pumping well increased with increasing period. In contrast, the diffusivity gained from analysis of pressure responses at the pumping and a monitoring well decreased with increasing period, as also found in field tests by Becker & Guiltinan (2010), Rabinovich et al. (2015) and Guiltinan & Becker (2015). The estimated diffusivity values of the two analyses agree at a critical period corresponding to a penetration depth on the order of the well distances. A similar phenomenon was observed in oscillatory hydraulic tests on rock samples (Song & Renner 2007). Fokker et al. (2013) successfully modeled the field observations using a heterogeneous subsurface suggesting that the disparity of hydraulic properties below the critical period reflects effective hydraulic properties of different parts of the subsurface probed depending on period. The work described here aimed at assessing the use of periodic testing for in situ hydraulic characterization of fault structures on the scale of metres to decametres which is of considerable interest in fundamental science and industrial applications. Previous fundamental research covered the relation between fault zone architecture as constrained by structure investigations (e.g. Wibberley & Shimamoto 2003), hydraulic properties (Caine et al.1996; López & Smith 1996) in cases derived from laboratory experiments (e.g. Mitchell & Faulkner 2009), mineralization (Sheldon & Micklethwaite 2007) and earthquake-source mechanics (Sibson et al.1975; Bizzarri 2012). From a perspective of industrial applications, the hydraulic characteristics of faults play an important role for assessing the performance of a reservoir in general (e.g. Fisher & Knipe 2001) but also for stimulations. The term hydraulic stimulation refers to an operation in wells initially developed by the hydrocarbon industry (Economides & Nolte 2000; Shaoul et al.2007) but now also used for developing deep geothermal reservoirs (Baria et al.1999; Schindler et al.2008; Genter et al.2010; Houwers et al.2015; Garcia et al.2016). We present and analyse results from periodic pumping tests in a well penetrating a strike-slip fault considered an analogue for a heat exchanger as needed for economic recovery of petrothermal energy. Our work addresses the need to develop testing protocols and evaluation methods suitable for heterogeneous reservoirs, such as faulted rocks. Equivalent hydraulic parameters are estimated by comparing the relation between flow rate and pressure with a range of previously presented and newly derived analytical solutions. The pressure record in the borehole section separated by packers is compared to that above or below the probe and the potential of the found amplitude ratios and phase shifts for the derivation of hydraulic parameters is explored relying on a numerical model. 2 EXPERIMENTAL DETAILS AND DATA ANALYSIS We use the term injectivity analysis when analysis involves flow rate and pressure in the actively tested section of the pumping well, addressed as interval. Advanced double-packer tests provide the pressure responses above and below the interval, too. We refer to the comparison of these various pressures measured in a single well as vertical-interference test following the convention of addressing the analyses of pressure responses in monitoring wells as interference tests (Matthews & Russell 1967). 2.1 Pumping tests 2.1.1 Location The close to E–W striking Grimsel Breccia Fault zone (GBFZ, Fig. 1) crops out at Grimselpass, Central Swiss Alps, with an altitude of 2160 m above sea level. To the north and the south, two shear zones subparallel to GBFZ, namely Grimsel Shear Zone (GSZ 2) and the Southern Boundary Shear Zone, crop out at distances of about 80 and 500 m, respectively (Fig. 1). The fault zone is hydrothermally active as evidenced by a 250 m wide zone of hydrothermal discharge in the Transitgas AG Tunnel that runs in NS direction at 200 m depth. The origin of the warm fluid is estimated at 4 ± 1 km depth with a reservoir temperature of at least 100 °C (Belgrano et al.2016). Borehole GDP-1 (coordinates of 669΄464/157΄025) with a radius of 0.048 m was drilled about 100 m west of the freshwater reservoir Totensee with a deviation of 24° from the vertical and an azimuth of 164° to intersect the GBFZ approximately normal to its strike. The Transitgas AG tunnel passes by the borehole at about 180 m distance to the west. The deviated borehole was drilled to a final depth of 125.3 m using a thixotropic polymer fluid. The true vertical depth is about 114.5 m. The borehole is open except for the casing stabilizing the first 2.1 m at the top of the well (Fig. 2). Figure 1. View largeDownload slide Sketch of the borehole site. The trajectory of borehole GDP-1 is indicated by the orange line. The depth of the Transitgas AG Tunnel (trajectory in light blue) is 200 m; the zone of inflow is indicated in dark blue. Estimated surface traces of shear zones and faults are represented by grey and dashed black lines, respectively. GBFZ: Grimsel Breccia Fault zone. GSZ 2: Grimsel Shear Zone. SBSZ: Southern Boundary Shear Zone. Figure 1. View largeDownload slide Sketch of the borehole site. The trajectory of borehole GDP-1 is indicated by the orange line. The depth of the Transitgas AG Tunnel (trajectory in light blue) is 200 m; the zone of inflow is indicated in dark blue. Estimated surface traces of shear zones and faults are represented by grey and dashed black lines, respectively. GBFZ: Grimsel Breccia Fault zone. GSZ 2: Grimsel Shear Zone. SBSZ: Southern Boundary Shear Zone. Figure 2. View largeDownload slide Sketch of the measurement setup. The flow rate recorded by the flow meter at the surface is denoted as Qsurface. Pressure recorded above the upper packer, in the interval and below the lower packer are represented as pabo, $${p_{{\mathop{\rm int}} }}$$ and pbel, respectively. Figure 2. View largeDownload slide Sketch of the measurement setup. The flow rate recorded by the flow meter at the surface is denoted as Qsurface. Pressure recorded above the upper packer, in the interval and below the lower packer are represented as pabo, $${p_{{\mathop{\rm int}} }}$$ and pbel, respectively. 2.1.2 Equipment and procedures The test system (Fig. 2) consisted of a reservoir filled with water from Totensee, a rotary pump and several magnetic-inductive flow meters covering a range from 0.35 to 400 L min−1 (accuracy better than 0.5 per cent of full scale) at the surface. A variable-length double-packer probe with downhole pressure gauges (full scale 2 MPa, accuracy 0.1 per cent of full scale) was lowered into the borehole. The packers with 1 m length were inflated from their initial diameter of 72 mm using nitrogen with pressures of up to 2.5 MPa. Data were digitally recorded with a resolution of 16 bit and a time step of 1 or 2 s. Four depth levels were selected for the testing relying on the preliminary results of borehole logging and inspection of drilled cores (Table 1 and Appendix  A). Tests with an interval length of 9.4 m at all four depth levels (i1–i4) were complemented by one test with a short interval of 1.7 m (i5) as part of i2. The latter isolated the prominent fault at a depth of ∼105 m with an orientation of 150/85 estimated from an optical log and core analyses (Appendix  A). The orientation of this fracture is actually also representative of the rather steeply dipping fractures striking subparallel to GBFZ that dominate in interval i1 and all of i2 (Table 1). Structural characterization of the main fault (i3) was not possible owing to extensive borehole breakout and core loss. The shallowest interval at ∼45 m (i4) exhibits a slightly higher fracture density than i2 and i3 mainly resulting from two perpendicular sets of E–W striking fractures, of which one is steeply dipping. Conventional tests (e.g. constant head, constant flow, slug, pulse and recovery tests) were performed at all intervals, while periodic pumping tests were conducted at all depth sections but i4. Table 1. Specifics of the test intervals. Int.  Packer position  Depth at centre of probe  Fracture density  Comment    Lower (m)  Upper (m)  (m)  (1 m−1)    i1  121.8  112.4  117.1  7.6  Several small breccia  i2  111.4  102.0  106.7  5.5  Prominent, isolated fault  i3  86.6  77.2  81.9  –*  Main fault with breccia  i4  49.6  40.2  44.9  8.8  Several fault strands  i5  105.8  104.1  105.0  6.4  Prominent, isolated fault  Int.  Packer position  Depth at centre of probe  Fracture density  Comment    Lower (m)  Upper (m)  (m)  (1 m−1)    i1  121.8  112.4  117.1  7.6  Several small breccia  i2  111.4  102.0  106.7  5.5  Prominent, isolated fault  i3  86.6  77.2  81.9  –*  Main fault with breccia  i4  49.6  40.2  44.9  8.8  Several fault strands  i5  105.8  104.1  105.0  6.4  Prominent, isolated fault  *Severe breakouts and core loss, neither logging (see Appendix  A) nor core analysis revealed structural information. View Large Initial selection of oscillation period and flow rate was based on forward calculations assuming radial flow in a homogeneous medium with diffusivity values considered typical for a fractured formation. The programme was then adjusted on the basis of the first results such that ultimately tests were performed at periods τ between 60 and 1080 s (Table 2). A single period consisted of two injection phases, each with a duration of half of the period, with nominal flow rates that differed by about 3 L min−1. By varying the lower rate in a period between 0 and 7 L min−1, we also modified the average pressure at which periodic testing took place (Fig. 3). According to the simple scaling relation $${z_p} \sim \sqrt {D\tau } $$, the chosen periods correspond to penetration depths zp of several decimetres to tens of metres for assumed values of hydraulic diffusivityD between 1E–5 and 1E–2 m2 s−1. Thus, neither Totensee nor the tunnel are expected to constitute relevant boundary conditions for the exerted flow. Figure 3. View largeDownload slide Pressures (left y-axis) and flow rate (green, right y-axis) recorded during the periodic testing at a depth of about 105 m (i5) for oscillation periods of (a) 60 s (i5–1, i5–2, see Table 2) and (b) 180 and 540 s (i5–3, i5–4). Figure 3. View largeDownload slide Pressures (left y-axis) and flow rate (green, right y-axis) recorded during the periodic testing at a depth of about 105 m (i5) for oscillation periods of (a) 60 s (i5–1, i5–2, see Table 2) and (b) 180 and 540 s (i5–3, i5–4). Table 2. Results of the Fourier analyses of recorded pressures and flow rate. Amplitude ratios and phase shifts are only reported when the amplitude spectrum confirmed that the imposed frequency is the dominant frequency in the signal. The results for tests that included sequences of zero-flow rate and therefore were manually shifted along the time axis are indicated by italic numbers. Int.  Period  Nominal flow rate  Injectivity analysis  Vertical-interference analysis        Qcor–pint  Qsurface–pint  pbel–pint  pabo–pint        Phase shift  Amplitude ratio  Amplitude ratio  Phase shift  Amplitude ratio  Phase shift  Amplitude ratio    (s)  (L min−1)  (cycles)  (m3 s−1 Pa−1)  (m3 s−1 Pa−1)  (cycles)  (–)  (cycles)  (–)  i1  180  3/0  3.65 ± 0.03E−01  4.54 ± 0.26E−10  2.05 ± 0.02E−09  4.43 ± 0.31E−01  6.90 ± 0.59E−03  6.14 ± 0.35E−01  1.40 ± 0.23E−03        1.08 ± 0.05E−01  5.21 ± 0.12E−10            i2–1  180  10/7  9.89 ± 0.22E−02  4.06 ± 0.09E−09  5.23 ± 0.08E−09  7.71 ± 0.25E−01  2.01 ± 0.30E−01  –  –  i2–2  180  3/0  1.70 ± 0.03E−01  3.76 ± 0.05E−09  5.26 ± 0.01E−09  8.74 ± 0.05E−01  3.44 ± 0.09E−01  –  –  i3–1  180  7/4  2.03 ± 0.05E−01  2.86 ± 0.30E−09  4.23 ± 0.09E−09  8.62 ± 0.03E−01  1.22 ± 0.04E−02  –  –  i3–2  180  3/0  3.24 ± 0.05E−01  1.06 ± 0.03E−09  2.67 ± 0.02E−09  8.91 ± 0.03E−01  1.51 ± 0.03E−02  –  –        2.65 ± 0.05E−01  9.88 ± 0.31E−10            i3–3  480  3/0  1.19 ± 0.02E−01  9.64 ± 0.10E−10  1.47 ± 0.02E−09  8.97 ± 0.04E−01  2.74 ± 0.07E−02  –  –  i3–5  1080  7/4  2.39 ± 0.08E−02  6.19 ± 0.20E−10  7.19 ± 0.23E−10  8.80 ± 0.10E−01  3.85 ± 0.17E−02  –  –  i5–1  60  3/0  3.80 ± 0.08E−01  5.19 ± 0.13E−09  9.43 ± 0.11E−09  9.46 ± 0.06E−01  1.53 ± 0.05E−01  –  –        1.40 ± 0.01E−01  4.90 ± 0.03E−09            i5–2  60  10/7  1.02 ± 0.03E−01  6.60 ± 0.13E−09  1.04 ± 0.01E−08  9.16 ± 0.17E−01  2.24 ± 0.20E−01  –  –  i5–3  180  10/7  7.40 ± 0.17E−02  4.65 ± 0.05E−09  5.63 ± 0.04E−09  8.94 ± 0.28E−01  2.85 ± 0.46E−01  7.04 ± 0.77E−01  7.60 ± 3.40E−03  i5–4  540  10/7  8.40 ± 0.11E−02  3.21 ± 0.03E−09  3.50 ± 0.04E−09  2.62 ± 1.67E−02  2.70 ± 0.15E−01  7.15 ± 0.51E−01  3.50 ± 0.59E−03  Int.  Period  Nominal flow rate  Injectivity analysis  Vertical-interference analysis        Qcor–pint  Qsurface–pint  pbel–pint  pabo–pint        Phase shift  Amplitude ratio  Amplitude ratio  Phase shift  Amplitude ratio  Phase shift  Amplitude ratio    (s)  (L min−1)  (cycles)  (m3 s−1 Pa−1)  (m3 s−1 Pa−1)  (cycles)  (–)  (cycles)  (–)  i1  180  3/0  3.65 ± 0.03E−01  4.54 ± 0.26E−10  2.05 ± 0.02E−09  4.43 ± 0.31E−01  6.90 ± 0.59E−03  6.14 ± 0.35E−01  1.40 ± 0.23E−03        1.08 ± 0.05E−01  5.21 ± 0.12E−10            i2–1  180  10/7  9.89 ± 0.22E−02  4.06 ± 0.09E−09  5.23 ± 0.08E−09  7.71 ± 0.25E−01  2.01 ± 0.30E−01  –  –  i2–2  180  3/0  1.70 ± 0.03E−01  3.76 ± 0.05E−09  5.26 ± 0.01E−09  8.74 ± 0.05E−01  3.44 ± 0.09E−01  –  –  i3–1  180  7/4  2.03 ± 0.05E−01  2.86 ± 0.30E−09  4.23 ± 0.09E−09  8.62 ± 0.03E−01  1.22 ± 0.04E−02  –  –  i3–2  180  3/0  3.24 ± 0.05E−01  1.06 ± 0.03E−09  2.67 ± 0.02E−09  8.91 ± 0.03E−01  1.51 ± 0.03E−02  –  –        2.65 ± 0.05E−01  9.88 ± 0.31E−10            i3–3  480  3/0  1.19 ± 0.02E−01  9.64 ± 0.10E−10  1.47 ± 0.02E−09  8.97 ± 0.04E−01  2.74 ± 0.07E−02  –  –  i3–5  1080  7/4  2.39 ± 0.08E−02  6.19 ± 0.20E−10  7.19 ± 0.23E−10  8.80 ± 0.10E−01  3.85 ± 0.17E−02  –  –  i5–1  60  3/0  3.80 ± 0.08E−01  5.19 ± 0.13E−09  9.43 ± 0.11E−09  9.46 ± 0.06E−01  1.53 ± 0.05E−01  –  –        1.40 ± 0.01E−01  4.90 ± 0.03E−09            i5–2  60  10/7  1.02 ± 0.03E−01  6.60 ± 0.13E−09  1.04 ± 0.01E−08  9.16 ± 0.17E−01  2.24 ± 0.20E−01  –  –  i5–3  180  10/7  7.40 ± 0.17E−02  4.65 ± 0.05E−09  5.63 ± 0.04E−09  8.94 ± 0.28E−01  2.85 ± 0.46E−01  7.04 ± 0.77E−01  7.60 ± 3.40E−03  i5–4  540  10/7  8.40 ± 0.11E−02  3.21 ± 0.03E−09  3.50 ± 0.04E−09  2.62 ± 1.67E−02  2.70 ± 0.15E−01  7.15 ± 0.51E−01  3.50 ± 0.59E−03  View Large 2.1.3 Undisturbed hydraulic heads We constrained the natural hydraulic heads from the pressure levels above, in, and below the interval asymptotically approached with time after packer inflation (Fig. 4). The recorded pressures suggest that the unperturbed water level of about 32 m below the surface, corresponding to about 700 L of water in the well, is actually the consequence of a dynamic situation with a prominent outflow at the borehole bottom, that is, at a true vertical depth of more than 110 m, and an inflow between about 73 and 110 m. Relative to the unperturbed pressure the pumping operations employed overpressures of at most 300 kPa (Fig. 4). Pumping pressures reached the level corresponding to a water column filling the entire borehole only for the interval i5. Figure 4. View largeDownload slide Equilibrium pressure of intervals (pint) and below probe (pbel) compared to pumping pressure in the interval and two hydrostatic water pressures for water levels at the surface (0 m) or at 32 m depth in the well. Figure 4. View largeDownload slide Equilibrium pressure of intervals (pint) and below probe (pbel) compared to pumping pressure in the interval and two hydrostatic water pressures for water levels at the surface (0 m) or at 32 m depth in the well. 2.2 Fourier analysis We performed two basic steps of data processing before Fourier transformation. Recorded time-series were cut to lengths corresponding to integer multiples of the pumping period and pressure records were automatically detrended to avoid artefacts due to differences between pressure levels at start and end of a pumping sequence (Fig. 5). Transformation yield amplitude and phase spectra from which we determined amplitude ratios and phase shifts between corrected flow rate (see Section 2.3.1) and interval pressure but also between interval pressure and pressure below and above the interval for the imposed period (see Table 2). The ratio between the amplitudes of flow rate and interval pressure corresponds to the injectivity index determined from conventional tests (e.g. Lyons 2010). Figure 5. View largeDownload slide Comparison of (a) recorded and (b) detrended pressure during the periodic testing at a depth of about 105 m (i5–3). For detrending, the average of the upper and lower envelopes is subtracted from the signal. Figure 5. View largeDownload slide Comparison of (a) recorded and (b) detrended pressure during the periodic testing at a depth of about 105 m (i5–3). For detrending, the average of the upper and lower envelopes is subtracted from the signal. To constrain the uncertainty of phase shift and amplitude ratio, we performed a sliding-window analysis, that is, a window with a length of three nominal periods was successively moved over the entire data set using a step size of 10–20 per cent of a period. For each window position, phase shift and amplitude ratio was determined and the standard deviation of all values is here reported as uncertainty. 2.3 Injectivity analysis 2.3.1 Correction for the storage capacity of the tubing Injectivity analysis requires the true flow rate into the rock that differs from the flow rate measured at the surface owing to the storage of fluid in the hydraulic system, composed of the tubing and the borehole section enclosed by the packers. For a closed system entirely filled with fluid, its capacity to store fluid stems from the finite fluid compressibility, the deformability of the tubing and of the pressurized rock section in the interval, and changes in packer seating due to changes in interval pressure. For an open system, one has to additionally account for a geometrical storage capacity that specifies how much fluid has to be added to/removed from the water column for a unit change in its height and thus in pressure it exerts. The used tubing with an inner diameter of 23.7 mm exhibits a storage capacity dominated by its geometrical component of $${S_{{\rm{tube}}}} = {{{A_{{\rm{tube}}}}} / {\rho g}} = 4.8 \times {10^{ - 8}}\,{{\rm{m}}^3}\,{\rm{P}}{{\rm{a}}^{ - 1}}$$ where Atube, ρ and g denote cross-section of the tubes (m2), fluid density (kg m−3) and gravitational acceleration (m s−2), respectively. The contribution from the compressibility of water is at least three orders of magnitude smaller, as is the contribution from the deformation of the steel tubes with a thickness of almost 5 mm. We constrained the combined contribution of rock deformation and changes in packer seating to interval-storage capacity by a rapid pulse test for the long interval of 9.4 m length (corresponding to an interval volume of about 60 L). The valve of the probe was briefly opened to pressurize the interval by the water column in the tubing. The change in volume was determined from the change in water-column height in the tubing and the change in pressure was recorded by the interval sensor. The ratio between these changes yields an interval-storage capacity of $${S_{{\mathop{\rm int}} }}( {9.4\ {\rm{m}}} ) = 6.0 \times {10^{ - 10}}\,{{\rm{m}}^3}\,{\rm{P}}{{\rm{a}}^{ - 1}}$$, that is, almost two orders of magnitude smaller than the geometrical storage capacity of the tubing (or 1.3 per cent). The shorter interval length of 1.7 m (corresponding to an interval volume of about 10 L) was not explicitly tested, but the reduced length should lead to a further reduction of the storage capacity roughly proportional to interval length, that is, Sint(1.7 m) ≃ 1.1 × 10 − 10m3Pa − 1. Thus, we neglected all contributions to storage capacity but the geometrical one of the tubing when determining true flow rates into or out of the rock formation from the flow rate recorded at the surface, Qsurface(L min−1), according to   \begin{equation}{Q_{\rm cor}} = {Q_{{\rm{surface}}}} - {S_{{\rm{tube}}}}{\dot{p}_{{\rm{int}}}}.\end{equation} (1) The time derivative of the interval pressure, $${\dot{p}_{{\rm{int}}}}$$, was calculated using a Fourier transformation of the recorded pressure. 2.3.2 Flow delay All tests with a water level in the tube below surface level are potentially affected by the conventional approach of remotely measuring flow rate at the surface and the associated time delay between flow in the gauge and actual addition of fluid to the water column loading the interval. Yet, when correcting flow records after eq. (1), tests employing a ‘zero’-flow rate, which we address as zero-flow rate tests below, yield suspicious delays of up to 10 s between flow rate at the surface and pressure in the interval exceeding (Fig. 6). These large delays presumably result from the actual time of flow in the partly empty tube. These data cannot be evaluated in a standard way. We pursue two different approaches to derive information on hydraulic properties from these tests. First, we determine the amplitude ratio of uncorrected flow rate, that is, surface-flow rate and interval pressure (Table 2) for all tests and compare their relation to pumping characteristics with that of amplitude ratios of the tests for which flow-rate correction was possible. Secondly, we test whether manual shifting of flow records guided by the prominent kinks in pressure records provides reasonable estimates for the phase shift between flow rate and pressure and thus opens the way for an inversion towards hydraulic parameters. Figure 6. View largeDownload slide Interval pressure (blue, left y-axis) and flow rate (orange, right y-axis) recorded during a periodic pumping test with a period of 180 s at a depth of 115 m (i1, see Tables 1 and 2) exhibiting a delay between start (stop) of pumping, indicated by the steep rise (fall) in flow rate, and increase in pressure, indicated by the kinks. Figure 6. View largeDownload slide Interval pressure (blue, left y-axis) and flow rate (orange, right y-axis) recorded during a periodic pumping test with a period of 180 s at a depth of 115 m (i1, see Tables 1 and 2) exhibiting a delay between start (stop) of pumping, indicated by the steep rise (fall) in flow rate, and increase in pressure, indicated by the kinks. 2.3.3 Flow regimes considered for estimation of hydraulic properties Equivalent or effective hydraulic properties, diffusivity D (m2 s−1), transmissivity T (m2 s−1) and storativity S (−), are typically estimated relying on analytical solutions of the pressure-diffusion equation for specific flow regimes, that is, the dimensional and directional characteristics of the flow pattern (e.g. Fetter 2001). Here, we considered: (i) 1-D flow, (ii) radial flow, (iii) radial flow in concentric cylindrical shells and (iv) bilinear flow. Analytical solutions of the 1-D-flow and radial-flow models were previously derived for various boundary conditions (Black & Kipp 1981; Rasmussen et al.2003; Renner & Messar 2006). Below, we investigate the effect of finite interval length on flow in homogeneous isotropic medium and present solutions for the other scenarios. Effect of interval length. The use of a double-packer probe poses a significant difference to previously performed periodic pumping tests mandating to investigate the consequences of a restricted ‘active’ length of the formation. For this purpose, we considered a radially infinite, mirror-symmetric model vertically bounded by a no-flow condition (see Appendix  B). To allow for analytical treatment, we assumed pressure to evolve linearly along the packers and to be constant above the upper (below the lower) packer. According to the model, the effects of interval length on the phase shift between flow rate and interval pressure are stronger the larger the period, that is, the penetration depth (Fig. 7). Relying on the radial-flow model for the inversion of hydraulic diffusivity from phase shift, diffusivity is the more overestimated the shorter the interval and the smaller the phase shift are. Small phase shifts correspond to large periods (small frequencies) for which a double-packer probe acts as a point source rather than a line source. Figure 7. View largeDownload slide Effect of interval length (as given in labels; radial flow corresponds to an infinitely long interval) on phase shift between flow rate and interval pressure according to the mirror-symmetric model described in Appendix  B. Curves were calculated after eq. (B17) with the summation up to n = 100. Figure 7. View largeDownload slide Effect of interval length (as given in labels; radial flow corresponds to an infinitely long interval) on phase shift between flow rate and interval pressure according to the mirror-symmetric model described in Appendix  B. Curves were calculated after eq. (B17) with the summation up to n = 100. Shell model. Analytical solutions for the injectivity analysis considering a shell model, known as multiregion composite model in reservoir engineering (Ambastha & Ramey 1992; Acosta & Ambastha 1994), are given in Appendix  C. These solutions complement the work by Ahn & Horne (2010), who reported a semi-analytical solution for interference analysis of a multicomposite radial ring system. For the analysis of the field data, we focus on a single cylindrical shell concentrically surrounding the borehole and exhibiting hydraulic properties that differ from the medium farther away from the well (comparable to the situation envisioned for the so-called skin effect, Matthews & Russell 1967). The prominent feature of this model is a perturbation of the transition from zero phase shift to a phase shift of 1/4 in comparison to radial flow in a homogeneous medium. The perturbations occur for values of ∼0.01–1 of the dimensionless parameter $${r_{\rm i}}\sqrt {\omega /{D_{{\rm{shell}}}}} $$, where Dshell is the hydraulic diffusivity of the shell, ri the radius of the injection well and ω the angular frequency. This dimensionless parameter gives the inverse of the ratio between the hydraulic penetration depth into the shell material and borehole radius. Prominent maxima, successions of maxima and minima, or a focusing of the transition occur depending on the ratios between the hydraulic parameters representing the shell and the surrounding medium (see e.g.  Appendix C2). Bilinear-flow model. Bilinear flow occurs when fluid is drained from or injected into a permeable matrix through an enclosed fracture of finite conductivity intersecting a well along its axis. A combination of two approximately linear flow regimes may result, one in the matrix with flow essentially perpendicular to the fracture and the other in the fracture itself associated with the non-negligible pressure drop or increase in it (Cinco-Ley et al.1978; Cinco-Ley & Samaniego-V 1981; Ortiz et al. 2013). In Appendix  D, we present an analytical solution of the coupled diffusion equations for bilinear flow in a homogeneous medium subjected to periodic pumping. For infinite fracture length, the phase shift is bounded by asymptotes to 1/16 and 1/8 of a cycle for large and small periods, respectively, with a smooth transition in-between (Fig. D2 in Appendix  D). The asymptotic value for long periods is consistent with the value derived by Hollaender et al. (2002). In case of fractures with finite length, only the asymptote of 1/8 of a cycle for small periods holds, too, as expected because the pressure perturbation does not reach the fracture tip. For a constant-pressure boundary at the fracture tip, phase shift monotonically decreases to 0 with increasing period. The relation between phase shift and period is not monotonic for a no-flow boundary, but minima (and possibly also maxima exceeding 1/8) occur at intermediate periods before 1/8 is asymptotically approached for very large periods when the pressure in the fracture is homogeneous and the total response is dominated by linear flow into the surrounding medium. 2.3.4 Flow-regime diagnosis from spectral analyses For conventional tests, flow regimes can be inferred from diagnostic log–log plots of pressure or pressure derivative versus time (e.g. Renard et al.2009). Analogously, slopes of amplitude spectra of non-harmonic periodic tests allow for a diagnosis of flow regimes (Hollaender et al.2002). In either type of diagram, slopes of 1/4, 1/2, 0 and 1 are indicative of bilinear flow, linear flow, radial flow and pseudo-steady-state flow, respectively. Since our excitation is not truly harmonic, spectra of flow rate and pressure contain significant amplitudes at periods shorter than the nominally excited one, too. Yet, simple division of entire flow-rate and pressure spectra proved to be impractical (Fig. 8). Thus, we restricted to specific local maxima of pressure and flow rate, as did Fokker et al. (2013). Figure 8. View largeDownload slide Spectrum of the amplitude ratio modulus for the interval at 105 m tested with a period of 540 s (i5–4). Coloured lines represent exponents of 1/4, 1/2 and 1. The restriction to local maxima of pressure and flow rate in the frequency domain (red asterisks) reduces the scatter and suggests bilinear (slope 1/4) to linear flow (1/2). Figure 8. View largeDownload slide Spectrum of the amplitude ratio modulus for the interval at 105 m tested with a period of 540 s (i5–4). Coloured lines represent exponents of 1/4, 1/2 and 1. The restriction to local maxima of pressure and flow rate in the frequency domain (red asterisks) reduces the scatter and suggests bilinear (slope 1/4) to linear flow (1/2). 2.4 Vertical-interference analysis When placing a double-packer probe in an open borehole, the pumping operations in the interval may also cause detectable pressure variations above or below the probe due to probe-parallel flow through the formation bearing additional information on its hydraulic properties. To take advantage of these pressure records, we analysed an axisymmetric model (Fig. 9) and treated it numerically using COMSOL. We determined values of phase shift and amplitude ratio of pressures measured at different locations of a double-packer probe considering different diffusivities and combinations of boundary conditions (Table 3). The effect of the investigated scenarios (1–8) is subordinate for the interference response between below lower packer and in the interval (Fig. 10). The results for pressure above upper packer separate into two groups according to the top boundary conditions labeled as ‘no flow’ and ‘constant pressure’. Figure 9. View largeDownload slide Illustration of the model used for numerical analysis of vertical interference. Blue colour indicates water-filled well sections and the two packers are represented in grey. The sketch is not to scale: distance of the lateral boundaries from the symmetry axis is 1000ri; distance between top and bottom boundary is 240 m, that is, the bottom boundary is more than 100 m away from the bottom of the borehole. Figure 9. View largeDownload slide Illustration of the model used for numerical analysis of vertical interference. Blue colour indicates water-filled well sections and the two packers are represented in grey. The sketch is not to scale: distance of the lateral boundaries from the symmetry axis is 1000ri; distance between top and bottom boundary is 240 m, that is, the bottom boundary is more than 100 m away from the bottom of the borehole. Figure 10. View largeDownload slide Phase shift and amplitude ratio of the interference response between pressure recorded (a) below the lower packer and (b) above the upper packer and in the interval calculated using the model described in Fig. 9 with different diffusivities and boundary conditions as specified in Table 3 (expressed by the numbers used as labels). The dots represent individual calculations labeled by the used diffusivity (m2 s−1). The lines intended to help the eye to follow the trends. The red square with error bars represents the observations at a depth of 85 m for a period of 180 s (i5–3). Figure 10. View largeDownload slide Phase shift and amplitude ratio of the interference response between pressure recorded (a) below the lower packer and (b) above the upper packer and in the interval calculated using the model described in Fig. 9 with different diffusivities and boundary conditions as specified in Table 3 (expressed by the numbers used as labels). The dots represent individual calculations labeled by the used diffusivity (m2 s−1). The lines intended to help the eye to follow the trends. The red square with error bars represents the observations at a depth of 85 m for a period of 180 s (i5–3). Table 3. Combinations of boundary conditions considered for the performed eight calculations employing the model sketched in Fig. 9. Boundary conditions  1  2  3  4  5  6  7  8  Top boundary  No flow  No flow  No flow  No flow  Constant pressure  Constant pressure  Constant pressure  Constant pressure  Bottom boundary  No flow  Constant pressure  Constant pressure  No flow  No flow  Constant pressure  Constant pressure  No flow  Lateral boundary  Constant pressure  Constant pressure  No flow  No flow  Constant pressure  Constant pressure  No flow  No flow  Boundary conditions  1  2  3  4  5  6  7  8  Top boundary  No flow  No flow  No flow  No flow  Constant pressure  Constant pressure  Constant pressure  Constant pressure  Bottom boundary  No flow  Constant pressure  Constant pressure  No flow  No flow  Constant pressure  Constant pressure  No flow  Lateral boundary  Constant pressure  Constant pressure  No flow  No flow  Constant pressure  Constant pressure  No flow  No flow  View Large Effective diffusivity values Dδ or Dφ are estimated from amplitude ratio or phase shift, respectively, by linearly interpolating between the modeling results. Deviation of model results from observations could be explained by several aspects. The model assumes homogeneity and isotropy for the subsurface which may not be true in reality. Furthermore, the model only considers relations between pressures but a more comprehensive model should also integrate the injectivity analysis. Finally, hydromechanical effects were observed but not considered in the current model. 3 RESULTS 3.1 Injectivity analysis Amplitude ratios of data sets resulting from performing the flow-rate correction after eq. (1) are systematically smaller than those of uncorrected data (Fig. 11), but the relative relations among the ratios for the different test intervals as well as their trends with the parameters characterizing a test, mean pressure and period, remain qualitatively similar. We therefore use the larger set of uncorrected data for the subsequent analysis and refer to this approach as simplified injectivity analysis. Figure 11. View largeDownload slide Injectivity quantified by the amplitude ratio of flow rate and injection pressure as dependent on mean interval pressure. Amplitude ratios were calculated for all tests before (labeled uncorrected and represented by solid circles) and after (labeled corrected and represented by open circles) applying the correction to flow rate for the storage of the tubing. The colour coding represents the five investigated test intervals (see Table 1). Data points with uncorrected flow rate are also labeled by the imposed oscillation period in seconds. Vertical error bars indicate the uncertainty in amplitude ratio gained from the sliding-window fast Fourier transform (FFT) analyses. Horizontal error bars represent the range of pressures at which the tests were conducted. Figure 11. View largeDownload slide Injectivity quantified by the amplitude ratio of flow rate and injection pressure as dependent on mean interval pressure. Amplitude ratios were calculated for all tests before (labeled uncorrected and represented by solid circles) and after (labeled corrected and represented by open circles) applying the correction to flow rate for the storage of the tubing. The colour coding represents the five investigated test intervals (see Table 1). Data points with uncorrected flow rate are also labeled by the imposed oscillation period in seconds. Vertical error bars indicate the uncertainty in amplitude ratio gained from the sliding-window fast Fourier transform (FFT) analyses. Horizontal error bars represent the range of pressures at which the tests were conducted. 3.1.1 Simplified injectivity analysis Amplitude ratios of uncorrected flow rates and interval pressures, here treated as tentative measures of injectivity, tend to increase with increasing mean pressure but variations at a given pressure are significant (Fig. 11). A period of 180 s was chosen for more than half of the periodic pumping tests allowing us to compare the hydraulic behaviour of the four tested intervals. The deepest interval with a depth of 115 m (i1) exhibits the lowest injectivity despite its association with the largest absolute mean pressure. The shallowest interval at 85 m depth (i3) shows a strong positive correlation between injectivity and mean pressure (also found for a period of 60 s applied in i5) indicating that the dominant hydraulic conduits are perceptible to hydromechanical effects. For a depth of 105 m (i2, i5), we neither find a pronounced dependence of amplitude ratio on mean pressure nor on interval length. Amplitude ratios decrease systematically with increasing pumping period for i5 as well as for i3. 3.1.2 Effective hydraulic parameters Effective transmissivity is the least variable of all hydraulic parameters when evaluating the results of the Fourier analyses of corrected flow rate and interval pressure relying on the analytical solution for radial flow in a homogeneous medium (Table 4). It varies by about half an order of magnitude from 3 to 9 × 10−6 m2 s−1. Apart from one outlier (i3–5), the values for hydraulic diffusivity and storativity stay within about one order of magnitude (4–40 × 10−5 m2 s−1 and 2–20 × 10−3, respectively). The outlier may to some extent be explained by the finite length of the interval. The modeling of the length effect (see Appendix  B) predicts the largest error for test i3–5, with diffusivity potentially overestimated by about one order of magnitude (0.37 m2 s−1 from conventional radial-flow model versus 0.04 m2 s−1 from mirror-symmetric double-packer model). For the other tests, yielding phase shifts of at least 0.074 cycles, the error due to neglecting the finite interval length is negligible (Fig. 7). Table 4. Hydraulic properties derived from the results of the Fourier analyses of flow rate and interval pressure relying on the analytical solution for radial flow in a homogeneous medium. Int.  Period  Nominal flow rate  Depth  Diffusivity  Transmissivity  Storativity    (s)  (L min−1)  (m)  (m2 s−1)  (m2 s−1)  (–)  i1  180  3/0  117.1  1.3 ± 0.9E−05  2.7 ± 0.8E−07  2.7 ± 1.3E−02  i2–1  180  10/7  106.7  3.9 ± 0.9E−05  3.5 ± 0.3E−06  9.4 ± 1.5E−02  i3–3  480  3/0  81.9  5.0 ± 3.1E−07  1.8 ± 0.5E−07  4.4 ± 1.8E−01  i3–5  1080  7/4  81.9  3.7 ± 1.3E−01  5.1 ± 0.2E−06  1.5 ± 0.5E−05  i5–2  60  10/7  105.0  8.4 ± 3.4E−05  4.9 ± 0.8E−06  6.4 ± 1.5E−02  i5–3  180  10/7  105.0  3.9 ± 0.3E−04  9.0 ± 0.2E−06  2.3 ± 0.2E−02  i5–4  540  10/7  105.0  5.7 ± 0.6E−05  4.8 ± 0.1E−06  8.5 ± 0.6E−02  Int.  Period  Nominal flow rate  Depth  Diffusivity  Transmissivity  Storativity    (s)  (L min−1)  (m)  (m2 s−1)  (m2 s−1)  (–)  i1  180  3/0  117.1  1.3 ± 0.9E−05  2.7 ± 0.8E−07  2.7 ± 1.3E−02  i2–1  180  10/7  106.7  3.9 ± 0.9E−05  3.5 ± 0.3E−06  9.4 ± 1.5E−02  i3–3  480  3/0  81.9  5.0 ± 3.1E−07  1.8 ± 0.5E−07  4.4 ± 1.8E−01  i3–5  1080  7/4  81.9  3.7 ± 1.3E−01  5.1 ± 0.2E−06  1.5 ± 0.5E−05  i5–2  60  10/7  105.0  8.4 ± 3.4E−05  4.9 ± 0.8E−06  6.4 ± 1.5E−02  i5–3  180  10/7  105.0  3.9 ± 0.3E−04  9.0 ± 0.2E−06  2.3 ± 0.2E−02  i5–4  540  10/7  105.0  5.7 ± 0.6E−05  4.8 ± 0.1E−06  8.5 ± 0.6E−02  View Large Phase-shift values are particularly diagnostic for hydraulic boundaries. For example, phase shifts larger than 1/8 unequivocally require bounded reservoirs (Fig. 12). The phase-shift value for the test at a depth of 85 m and a period of 180 s (i3–2) remains larger than the upper bound of 1/4 of all considered bounded models even after manual shifting. Only models with a no-flow boundary can explain the phase shift of the test i3–1 at 85 m and a period of 180 s, the largest phase shift observed for non-zero-flow rate tests (Fig. 12). However, the increase of amplitude ratios with phase shift between i3–1 and the other test at the same depth for a period of 1080 s (i3–5) cannot be explained by any of the no-flow models. Independent of interval length the phase-shift values for a depth of ∼105 m (i2, i5) are all smaller than 1/8 of a cycle, excluding the 1-D-flow model with a no-flow boundary or no boundary. The non-monotonically varying amplitude ratios for i5 also exclude a constant pressure boundary or no boundary for all the flow types. For the interval with a depth of ∼115 m (i1), the only phase shift lies between 1/16 and 1/8, which can be fit by all the models except 1-D flow with no-flow boundary or no boundary. Figure 12. View largeDownload slide Amplitude ratio versus phase shift for the injectivity analysis. Red, blue and black lines represent the solutions for radial flow, 1-D flow and bilinear flow, respectively. Solid lines are for infinite reservoirs and dashed and dashed–dotted lines represent constant pressure (pb) and no-flow (fb) boundary at distances of 20ri, respectively, for radial and bilinear flows. Changing the assumed hydraulic parameters shifts the theoretical curves vertical without distortion. Filled and open symbols represent non-zero-flow-rate sequences and zero-flow-rate sequences for which manual shift was applied, respectively (see Table 2). Errors for the data points do not exceed symbol size. Two periods are also labeled for the convenience of discussion (see Table 4). Figure 12. View largeDownload slide Amplitude ratio versus phase shift for the injectivity analysis. Red, blue and black lines represent the solutions for radial flow, 1-D flow and bilinear flow, respectively. Solid lines are for infinite reservoirs and dashed and dashed–dotted lines represent constant pressure (pb) and no-flow (fb) boundary at distances of 20ri, respectively, for radial and bilinear flows. Changing the assumed hydraulic parameters shifts the theoretical curves vertical without distortion. Filled and open symbols represent non-zero-flow-rate sequences and zero-flow-rate sequences for which manual shift was applied, respectively (see Table 2). Errors for the data points do not exceed symbol size. Two periods are also labeled for the convenience of discussion (see Table 4). At a depth of ∼105 m, three tests were performed with different periods using the short interval (60, 180 and 540 s corresponding to i5–2, i5–3, and i5–4 in Table 2). Since they all included only non-zero flow rates their results are not affected by the problems associated with the finite distance between the location of the flow meter and the top of the water column loading the interval. The non-monotonic succession of phase-shift values, with the lowest occurring for the intermediate period, can be modeled by a shell with a thickness of about 6–15 times the borehole radius of ∼0.05 m, and with a diffusivity between just below 10−3 to above 10−2 m2 s−1 (Table 5) exceeding that for the enclosing medium by a factor between about 5 and 10 (Fig. 13). Figure 13. View largeDownload slide Examples phase-shift–frequency relations for the shell model calculated using eq. (C2) in relation to observations for the periodic tests at 105 m for a period of 180 s (i5–3) and for a period of 540 s (i5–4). The labels give the ratios between the diffusivity of the shell and the surrounding medium and the lateral extension of the shell in multiples of the borehole radius. Figure 13. View largeDownload slide Examples phase-shift–frequency relations for the shell model calculated using eq. (C2) in relation to observations for the periodic tests at 105 m for a period of 180 s (i5–3) and for a period of 540 s (i5–4). The labels give the ratios between the diffusivity of the shell and the surrounding medium and the lateral extension of the shell in multiples of the borehole radius. Table 5. Parameters of the shell model (see Appendix  C) derived from the phase-shift values observed for the tests i5–3 and i5–4 (see also Fig. 13). k1/k2  D1/D2  r1/ri  Diffusivity (D1)  (–)  (–)  (–)  (m2 s−1)  6.7  8  15  8.1E−04  7.1  10  15  7.4E−04  6.3  10  10  7.1E−03  5.9  9  6.1  2.3E−02  5.6  9  6.1  2.1E−02  k1/k2  D1/D2  r1/ri  Diffusivity (D1)  (–)  (–)  (–)  (m2 s−1)  6.7  8  15  8.1E−04  7.1  10  15  7.4E−04  6.3  10  10  7.1E−03  5.9  9  6.1  2.3E−02  5.6  9  6.1  2.1E−02  View Large 3.2 Flow regime diagnosis from spectral analyses Flow-regime identification from full spectra proved difficult (Fig. 8), not the least because the spectra of flow rate and pressure contain little significant contributions for high frequencies. We thus restricted to the first two odd multiples of the test frequency that correspond to pronounced local maxima in the spectra. As the analysis is sensitive to the timing of flow-rate changes (Hollaender et al.2002), we only analyse non-zero-flow-rate sequences. The results for overtones from the various tests in a specific interval agree closely with those for nominally excited periods (Fig. 14). The slope for the interval at 105 m (i5) increases from less than 1/4 and approaches 1/2 as period increases. In contrast, the slope for interval i3 decreases from around 1 to a value just above 1/4. At face value, these results indicate a transition from bilinear to linear flow for i5 and a transition from pseudo-steady-state flow to bilinear flow for i3. The physical interpretation for a transition from bilinear to linear flow is that the pressure at the fracture tip has risen to a substantial fraction of the associated pressure increase in the well at the end of bilinear flow and the flow field transforms towards formation linear flow. The transition from pseudo-steady-state flow to bilinear flow is consistent with a succession of the domination of wellbore storage at early time by bilinear flow (finite conductivity fracture). Figure 14. View largeDownload slide Amplitude ratio in frequency domain for the deepest interval (i5) represented by asterisks and the interval at depth 81.9 m (i3) represented by circles at the main period of the pumping tests and the first overtone, that is, a third of the imposed period. The dotted lines represent exponents of 1/4, 1/2 and 1. Figure 14. View largeDownload slide Amplitude ratio in frequency domain for the deepest interval (i5) represented by asterisks and the interval at depth 81.9 m (i3) represented by circles at the main period of the pumping tests and the first overtone, that is, a third of the imposed period. The dotted lines represent exponents of 1/4, 1/2 and 1. 3.3 Vertical-interference analysis The periodic excitations in the injection intervals lead to detectable periodic responses in the pressures recorded below the lower and above the upper packer (Table 2). The vertical-interference analysis based on these pressure records has the advantage that all tests can be analysed even those involving the problematic zero-flow-rate sequences. Yet, the magnitude of interference response in a single borehole critically depends on the storage capacity of the borehole section in which it is determined. Above the upper packer we had a free, unconstrained water column corresponding to a ‘large’ storage capacity determined by the cross-section of the borehole. The borehole section below the lower packer, in contrast, constitutes an enclosed fluid volume with a ‘small’ storage capacity determined by the fluid compressibility and the deformability of packer and borehole wall. In agreement with these qualitative storage considerations, the interference response is much larger below the probe than above the probe where it is in fact detectable only in some cases (Table 2). All amplitude ratios derived from pressures below the lower packer, but that for the deepest test interval (i1), follow the expected inverse correlation with the length of the enclosed borehole section (Fig. 15): the larger the enclosed volume the larger the storage capacity and the less sensitive the section, that is, the amplitude ratio between interference response and excitation decreases with the length of the enclosed section. Relative to the general trend, the response below the probe observed for the test in the deepest interval is weak demonstrating exceptional hydraulic properties near the borehole bottom. Figure 15. View largeDownload slide Correlation between interference responses, quantified by the amplitude ratio of pressure oscillations, and the length of the enclosed section below the lower packer of the probe. We present the amplitude ratio between the pressure above the interval and in the interval in this graph, too, though it should not be controlled by the length of the enclosed section below the probe. The joint presentation of results is chosen to highlight the difference in sensitivity of vertical response between the enclosed section below the probe and the unconstrained water column above it. Figure 15. View largeDownload slide Correlation between interference responses, quantified by the amplitude ratio of pressure oscillations, and the length of the enclosed section below the lower packer of the probe. We present the amplitude ratio between the pressure above the interval and in the interval in this graph, too, though it should not be controlled by the length of the enclosed section below the probe. The joint presentation of results is chosen to highlight the difference in sensitivity of vertical response between the enclosed section below the probe and the unconstrained water column above it. Including phase shift in the comparison, observations for the interference with the borehole section above the interval are only in agreement with model results for a ‘no-flow’ condition (Fig. 10). The effective diffusivities estimated from phase shift (denoted Dφ in Table 6) tend to be larger than the ones estimated from amplitude ratio (Dδ) for i2 and i5, the long and short intervals at 105 m. The reverse seems to hold for the shallowest interval i3. Effective diffusivity values using pbel–pint data or consistently on the order of a few m2 s−1 for the long and the short interval (i2 and i5) at the same depth, except for two relatively small values (i2–1, i5–4). Those estimated by pabo–pint analysis are smaller than those estimated by pbel–pint for the short interval i5 without considering the extremely small Dφ value for i5–4. A decreasing trend of diffusivity with increasing period is indicated for i5. Table 6. Effective diffusivity Dδ estimated by amplitude ratio and Dφ by phase shift of vertical-interference analysis. Top and bottom boundaries are set as no flow. Int.  Period   Nominal flow rate  Dδ (m2 s−1)  Dφ (m2 s−1)    (s)  (L min−1)  pbel–pint  pabo–pint  pbel–pint  pabo–pint  i1  180  3/0  5E−3  5E−3  3E−3  4E−3  i2–1  180  10/7  1  –  3E−2  –  i2–2  180  3/0  4  –  2  –  i3–1  180  7/4  2E−2  –  3E−1  –  i3–2  180  3/0  3E−2  –  4E−1  –  i3–3  480  3/0  3E−2  –  3E−1  –  i3–5  1080  7/4  2E−2  –  1E−1  –  i5–1  60  3/0  4  –  8  –  i5–2  60  10/7  7  –  6  –  i5–3  180  10/7  5  9E−2  3  3E−2  i5–4  540  10/7  3  1E−2  3E−4  7E−3  Int.  Period   Nominal flow rate  Dδ (m2 s−1)  Dφ (m2 s−1)    (s)  (L min−1)  pbel–pint  pabo–pint  pbel–pint  pabo–pint  i1  180  3/0  5E−3  5E−3  3E−3  4E−3  i2–1  180  10/7  1  –  3E−2  –  i2–2  180  3/0  4  –  2  –  i3–1  180  7/4  2E−2  –  3E−1  –  i3–2  180  3/0  3E−2  –  4E−1  –  i3–3  480  3/0  3E−2  –  3E−1  –  i3–5  1080  7/4  2E−2  –  1E−1  –  i5–1  60  3/0  4  –  8  –  i5–2  60  10/7  7  –  6  –  i5–3  180  10/7  5  9E−2  3  3E−2  i5–4  540  10/7  3  1E−2  3E−4  7E−3  View Large 3.4 Conventional methods Conventional tests yield the smallest transmissivity and storativity for the shallowest interval (i4, Table 7). Disregarding the exceptionally low value found for interval i4, effective transmissivity is the least variable hydraulic property, consistent with the observations from periodic radial-flow analysis. The hydraulic parameters derived from the conventional methods and periodic radial-flow analysis are of similar order of magnitude (Table 4). Flow regimes were constrained using diagnostic log–log plots of pressure and pressure derivative versus elapsed time. Intervals i1 und i3 indicate bilinear flow and interval i2 and i5 display features of linear flow. Table 7. Hydraulic properties derived from conventional tests. Int.  Test phase  Diffusivity (m2 s−1)  Transmissivity (m2 s−1)  Storativity (–)  Remarksa  i1  Shut-in   <1.14E−04   <1.6E−06  1.4E−02  RSLA  i2  Constant flow rate   <1.39E−04   <4.3E−06  3.1E−02  SLA  i3  Shut-in   <2.29E−04   <1.9E−06  8.3E−03  RSLA  i4  Shut-in   <1.45E−03   <2.9E−08  2.0E−05  RSLA  i5  Shut-in   <3.18E−05   <2.7E−06  8.5E−02  RSLA  Int.  Test phase  Diffusivity (m2 s−1)  Transmissivity (m2 s−1)  Storativity (–)  Remarksa  i1  Shut-in   <1.14E−04   <1.6E−06  1.4E−02  RSLA  i2  Constant flow rate   <1.39E−04   <4.3E−06  3.1E−02  SLA  i3  Shut-in   <2.29E−04   <1.9E−06  8.3E−03  RSLA  i4  Shut-in   <1.45E−03   <2.9E−08  2.0E−05  RSLA  i5  Shut-in   <3.18E−05   <2.7E−06  8.5E−02  RSLA  aSLA: ‘Straight Line Analysis’ after Jacob & Lohman (1952). RSLA: ‘Recovery Straight Line Analysis’ after Agarwal (1980). View Large 4 DISCUSSION The variable differences between equilibrium pressure for a certain depth interval and a hydrostatic gradient indicate heterogeneity of the hydraulic system. The fractured aquifer is dominated by steep conduits that appear to have limited interconnectivity despite their small lateral distance (at most a few metres in the case of the borehole's bottom section). The diffusivity variations with pumping periods suggest spatial heterogeneity or complexity of the conduits in a specific interval, too, that is, simple radial flow is a poor approximation of the flow excited during the tests. The prominent role of discrete hydraulic conduits is documented by the similarity of injectivity values for the two tests at a borehole depth of about 105 m (i2 and i5) despite their differences in interval length. In the following, we discuss to what extent the employed methods yield comparable values for hydraulic parameters, comment on the constraints for the pressure dependence of the hydraulic response, and finally speculate on the flow-regime model that is most consistent with the entirety of observations in a specific interval. 4.1 Comparison of methods Periodic pumping tests can be easily implemented using field equipment for conventional testing. The excitation signals do not have to be perfectly harmonic but Fourier analysis of non-harmonic signals might actually provide information regarding periods shorter than the imposed main period (spectral analyses, e.g. Hollaender et al.2002; Fokker et al.2013). The periodic excitation can be applied even when the hydraulic pressure is not equilibrated, in fact it can be superposed to any transient, be it associated with a terminated or ongoing pumping operation. The superposition is valid as long as linearity can be assumed, that is, the diffusion equation in its basic form holds. If it does not hold, the results of conventional methods become questionable, too. As a consequence of the freedom to start periodic testing at any time, operational time is reduced and to a larger degree plannable than for conventional tests that involve long waiting times for equilibration of a priori unknown duration. During the current test campaign, a total of 11 hr was spent on conventional tests and 7 hr for periodic tests. Our results show that the analyses of periodic pumping provide effective hydraulic parameters comparable to those gained from conventional methods when assuming radial flow (Fig. 16). Only for the shallowest of the tested intervals (i3), diffusivity values derived from periodic testing deviate significantly from the one gained from conventional testing. The lower diffusivity value results from a test sequence with zero flow rate (i3–3) and thus has a limited reliability. The higher diffusivity value corresponds to test i3–5 giving a phase shift so low that it falls within the range for which the finite length of the probe cannot be neglected, but the flow field has a significant axial component (Fig. 7). The estimate based on the mirror-symmetric double-packer model is more than one order of magnitude smaller than the value assuming simple radial flow. Furthermore, the analyses of periodic tests demonstrates a dependence of the effective hydraulic parameters on period and mean pressure. Thus, differences between the methods likely reflect the simplifications inherent in the radial-flow model that are not justified for the investigated structures necessitating to investigate more complex flow models. Figure 16. View largeDownload slide (a) Radial diffusivity determined from periodic injectivity analysis and conventional methods. (b) Axial diffusivity determined by vertical-interference analysis. Labels refer to the effective diffusivity values estimated from amplitude ratio (Ddel: Dδ) and phase shift (Dphi: Dφ) using pressure records above (pabo), in (pint), and below (pbel) the interval. Figure 16. View largeDownload slide (a) Radial diffusivity determined from periodic injectivity analysis and conventional methods. (b) Axial diffusivity determined by vertical-interference analysis. Labels refer to the effective diffusivity values estimated from amplitude ratio (Ddel: Dδ) and phase shift (Dphi: Dφ) using pressure records above (pabo), in (pint), and below (pbel) the interval. Vertical-interference analysis is a valuable by-product of periodic testing when packer-separated intervals are investigated and pressure is monitored at several depth levels. We refer to its results as effective axial diffusivity to distinguish them from results of the injectivity analysis and conventional methods, addressed as effective radial diffusivity. The attributes ‘axial’ and ‘radial’ are chosen from the perspective of the pumping well, addressing flow dominantly subparallel and normal to well axis, respectively. Effective axial diffusivities are here found to be generally larger than effective radial diffusivities (Fig. 16). The axial diffusivities scatter significantly for the intervals i2 and i5 at similar mean depths but with different probe lengths. The scatter is larger for the shorter interval i5 tested for a larger spread in periods and apart from the result for the longest period of 540 s, communication to the deeper borehole section is characterized by larger diffusivity values than to the shallower sections. The vertical response to shallower depth sections is actually lost with the increase in probe length suggesting that faults outside the short interval but inside the long interval are hydraulically connected to the prominent fault isolated by the short interval. The decrease in diffusivity with increasing period or penetration depth (Table 6) also indicates that the axial connectivity is a local phenomenon. 4.2 Pressure dependence of hydraulic properties The tested intervals exhibit a general trend of increasing injectivity with increasing pressure, with the prominent exception of the deepest interval (i1) that also has the lowest injectivity. In the current geological situation, an increase in injectivity with mean pressure likely results from the pressure-induced increase in effective hydraulic aperture of the joints and faults intersecting the borehole. The more pronounced pressure dependence of the short interval i5 compared to the long interval i2 supports the notion that the pressure dependence reflects hydromechanical behaviour of discrete faults. Observations for i3 (∼85 m depth) with a relatively large period of 180 s suggest that the hydromechanical effects are not restricted to the very intersection between faults and the well. The significance of hydromechanical effects is also evidenced by the pronounced reverse pressure response observed below the lower packer when injecting in i1 (not documented here). 4.3 Constraints on flow regime Details of our observations are obviously at conflict with simple radial flow. From periodic tests, flow regimes can be constrained by (i) investigating spectra of non-harmonic signals (Fig. 14) or by (ii) comparing the variation patterns of amplitude ratio and phase shift with theoretical flow models (Fig. 12). The few constraints we derived from the analysis of spectra yield flow regimes identical to those indicated by the conventional methods. The periodic injectivity analysis tends to indicate more complex flow regimes than the conventional methods (Table 8) which we interpret as evidence for the superior resolution of the periodic approach. Table 8. Constraints on flow regimes gained from injectivity analysis, spectral analysis, or conventional methods: possible (o), excluded (x), no constraint (–) for radial (R), bilinear (BL) and 1-D (1-D) flow.     Periodic pumping method  Conventional methods      Injectivity analysis  Spectral analysis (see Fig. 14)        No flow  Constant pressure  No boundary  No boundary  No boundary  i1  R  o  o  o  –  x    1-D  x  o  x  –  x    BL  o  o  o  –  o  i2  R  o  o  o  –  x    1-D  x  o  x  –  o    BL  x  x  x  –  x  i3  R  x  x  x  x  x    1-D  x  x  x  x  x    BL  x  x  x  o  o  i5  R  o  x  x  x  x    1-D  x  x  x  o  o    BL  o  x  x  o  x      Periodic pumping method  Conventional methods      Injectivity analysis  Spectral analysis (see Fig. 14)        No flow  Constant pressure  No boundary  No boundary  No boundary  i1  R  o  o  o  –  x    1-D  x  o  x  –  x    BL  o  o  o  –  o  i2  R  o  o  o  –  x    1-D  x  o  x  –  o    BL  x  x  x  –  x  i3  R  x  x  x  x  x    1-D  x  x  x  x  x    BL  x  x  x  o  o  i5  R  o  x  x  x  x    1-D  x  x  x  o  o    BL  o  x  x  o  x  View Large Effective injectivity values determined relying on simple radial flow (Fig. 11) systematically decrease with increasing pumping period for intervals i3 and i5 (factor of 3–5 for less than one order of magnitude in period). The conventional scaling relation for diffusion processes (e.g. Weir 1999) suggests that the covered range in period corresponds to a change of half an order of magnitude in nominal penetration depth. Thus, the period dependence of injectivity indicates a change in hydraulic characteristics with distance from the borehole. When addressing the period dependence by the shell model for the deepest interval (i5), we find a decrease in effective diffusivity and permeability by about one order of magnitude beyond a zone of about several decimetres to one metre (Fig. 13). This change in hydraulic characteristics could be related to drilling-associated damage of the borehole wall or remnants of the thixotropic polymer fluid used during drilling. The flushing of the well during logging preceding the hydraulic testing may have well resulted in a cleansing of only the hydraulic conduits nearest to the borehole. For the short interval i5, the conventional method and the spectral analysis (Fig. 14) suggest a linear flow regime while injectivity analysis indicates radial or bilinear flow with a no-flow boundary. For the longer interval at the same depth (i2), the injectivity analysis suggests radial flow with any boundary condition or linear flow with constant pressure boundary. Assuming the two intervals to be governed by the same boundary condition, the only regime in common by i5 and i2 is radial flow with a no-flow boundary. For the deepest interval (i1), the injectivity is the lowest of all intervals though it was tested at the highest mean injection pressure. The low equilibrium pressure below i1 (Fig. 4) suggests that the section below i1 is disconnected from the ones above. In addition, higher axial diffusivity than radial diffusivity (Fig. 16) implies axially oriented or subvertical upflow in the sections above i1. 4.4 Synopsis The periodic testing revealed that the flow geometry in the penetrated subsurface is quite variable on the metre (lateral) to decametre (axial) scale and that the conduit system exhibits hydromechanical effects, two characteristics that appear quite reasonable for a strike-slip fault in crystalline rock. Yet, at this point, it is difficult to distinguish between a situation characterized by a lateral no-flow boundary, possibly related to relics of the polymer drilling fluid with high viscosity or drilling-associated near-welbore damage, and a truly anisotropic hydraulic system with the high axial diffusivity due to subvertically oriented fluid pathways. The latter notion is not only supported by the logging observations (Appendix  A) that document a dominance of steeply dipping fractures intersecting the well, but is also consistent with the conduit geometry deduced by Belgrano et al. (2016), who suggested that the fault zone is controlled by localized subvertically oriented ‘pipe’-like upflow zones. The vertical variability observed for the undisturbed hydraulic heads and the low hydraulic transmissivity of the deepest test interval compared to the major fault zone support this suggestion derived solely from structural investigations. The performed pumping tests correspond to nominal lateral penetration depths of up to a few tens of metres, that is, an order of magnitude smaller than the length scale of a subsurface heat exchanger required for an economically relevant petrothermal energy system. The hydraulic properties near the borehole, however, control the very access to a heat exchanger and thus their characterization, as performed here, bears important consequences for operation. The complexity of flow between formation and well indicated by our results bears significant implications for the average near-well pressure gradient and thus for effective inlet-pressure drops (e.g. Haraden 1992) that in turn bear consequences for operational parameters and bulk flow rates. Also, the realization of the relevance of hydromechanical effects should find observance in modeling and simulation, since their operation is not limited to the vicinity of boreholes (Vinci et al.2015). 5 CONCLUSIONS Periodic pumping tests were conducted using a double-packer probe placed at four different depth levels in a borehole penetrating a hydrothermally active fault. The results of the periodic testing, expressed as amplitude ratio and phase shift between recorded pressures and/or flow rate, permit to constrain effective hydraulic parameters and the flow geometry. We performed and extended the previously established injectivity analysis, but also introduced a new approach, addressed as vertical-interference analysis. We demonstrated that the analyses of the periodic pumping tests yield results comparable to conventional methods when relying on the same model. The field campaign revealed, however, several methodological advantages of the periodic testing, for example, its plannable testing time and the ease in separating the response to the actual pumping operation from noise and transients gained from data processing in frequency space. Our analytical and numerical modeling of periodic testing extends the suite of evaluation tools and specifically allowed us to decipher more details of the encountered complex flow geometry and further conduit characteristics than the conventional methods. The two approaches, the previously employed injectivity analysis and the vertical-interference analysis, are complementary by providing information on the hydraulic properties in different directions relative to the well axis. The variations in effective ‘axial’ and ‘radial’ diffusivity with pumping periods and mean pumping pressure indicate significant spatial heterogeneity of conduits that also exhibit hydromechanical effects for the investigated site. While we showed the potential of the vertical-interference analysis, clearly more work is needed in the future regarding the modeling of records available when using double- or multipacker probes. ACKNOWLEDGEMENTS The financial support by BfE, Switzerland is gratefully acknowledged as is the coordinator, Marco Herwegh, for involving us in project NFP70-P1. Field tests were performed together with Sacha Reinhardt and Markus Bosshard (SolExperts), and Daniel Eggli (Uni Bern) who is particularly thanked for generously sharing the results of the structural core and log analyses. The presented treatment of the shell model builds on an unpublished analyses of Eugen Petkau. We are grateful to the three anonymous reviewers for their valuable comments. REFERENCES Acosta L.G., Ambastha A.K., 1994. Thermal well test analysis using an analytical multi-region composite reservoir model, in SPE Annual Technical Conference and Exhibition , 25–28 September, New Orleans, Louisiana. doi:10.2118/28422-MS. Agarwal R.G., 1980. A new method to account for producing time effects when drawdown type curves are used to analyze pressure buildup and other test data, in SPE Annual Technical Conference and Exhibition , 21–24 September, Dallas, Texas. doi:10.2118/9289-MS. Ahn S., Horne R.N., 2010. Estimating permeability distributions from pressure pulse testing, in SPE Annual Technical Conference and Exhibition , 19–22 September, Florence, Italy. doi:10.2118/134391-MS. Ambastha A.K., Ramey H.J.Jr., 1992. Pressure transient analysis for a three-region composite reservoir, in SPE Rocky Mountain Regional Meeting , 18–21 May, Casper, Wyoming. doi:10.2118/24378-MS. Bakhos T., Cardiff M., Barrash W., Kitanidis P.K., 2014. Data processing for oscillatory pumping tests, J. Hydrol. , 511, 310– 319. Google Scholar CrossRef Search ADS   Baria R., Baumgärtner J., Rummel F., Pine R.J., Sato Y., 1999. HDR/HWR reservoirs: concepts, understanding and creation, Geothermics , 28, 533– 552. Google Scholar CrossRef Search ADS   Becker M.W., Guiltinan E., 2010. Cross-hole periodic hydraulic testing of inter-well connectivity, in Proceedings Thirty-Fifth Workshop on Geothermal Reservoir Engineering, California , Vol. 35, pp. 292– 297, Stanford University, Stanford Geothermal Program. Belgrano T.M., Herwegh M., Berger A., 2016. Inherited structural controls on fault geometry, architecture and hydrothermal activity: an example from Grimsel Pass, Switzerland, Swiss J. Geosci. , 109, 345– 364. Google Scholar CrossRef Search ADS   Bizzarri A., 2012. Effects of permeability and porosity evolution on simulated earthquakes, J. Struct. Geol. , 38, 243– 253. Google Scholar CrossRef Search ADS   Black J.H., Kipp K.L., 1981. Determination of hydrogeological parameters using sinusoidal pressure tests: a theoretical appraisal, Water Resour. Res. , 17, 686– 692. Google Scholar CrossRef Search ADS   Bourdet D., Ayoub J., Pirard Y., 1989. Use of pressure derivative in well test interpretation, Soc. Petrol. Eng. Format. Eval. , 4, 293– 302. Caine J.S., Evans J.P., Forster C.B., 1996. Fault zone architecture and permeability structure, Geology , 24, 1025– 1028. Google Scholar CrossRef Search ADS   Cardiff M., Bakhos T., Kitanidis P.K., Barrash W., 2013. Aquifer heterogeneity characterization with oscillatory pumping: sensitivity analysis and imaging potential, Water Resour. Res. , 49, 5395– 5410. Google Scholar CrossRef Search ADS   Carslaw H.S., Jaeger J.C., 1986. Conduction of Heat in Solids , 2nd edn, pp. 510, Oxford Univ. Press. Cinco-Ley H., Samaniego-V F., 1981. Transient pressure analysis for fractured Wells, J. Petrol. Technol. , 33, 1749– 1766. Google Scholar CrossRef Search ADS   Cinco-Ley H., Samaniego-V F., Dominguez-A N., 1978. Transient pressure behavior for a well with a finite-conductivity vertical fracture, Soc. Petrol. Eng. J. , 18, 253– 264. Google Scholar CrossRef Search ADS   Economides M.J., Nolte K.G., 2000. Reservoir Stimulation , 3rd edn, pp. 856, Wiley and Sons Ltd. Fetter C.W., 2001. Applied Hydrogeology , 4th edn, pp. 598, Prentice Hall. Fisher Q.J., Knipe R.J., 2001. The permeability of faults within siliciclastic petroleum reservoirs of the North Sea and Norwegian Continental Shelf, Mar. Petrol. Geol. , 18, 1063– 1081. Google Scholar CrossRef Search ADS   Fokker P.A., Renner J., Verga F., 2012. Applications of harmonic pulse testing to field cases, in SPE Europec/EAGE Annual Conference , 4–7 June, Copenhagen, Denmark. doi:10.2118/154048-MS. Fokker P.A., Renner J., Verga F., 2013. Numerical modeling of periodic pumping tests in wells penetrating a heterogeneous aquifer, Am. J. Environ. Sci. , 9, 1– 13. Google Scholar CrossRef Search ADS   Garcia J., Hartline C., Walters M., Wright M., Rutqvist J., Dobson P.F., Jeanne P., 2016. The Northwest Geysers EGS Demonstration Project, California: Part 1: Characterization and reservoir response to injection, Geothermics , 63, 97– 119. Google Scholar CrossRef Search ADS   Genter A., Goerke X., Graff J.-J., Cuenot N., Krall G., Schindler M., Ravier G., 2010. Current status of the EGS Soultz geothermal project (France), in Proceedings World Geothermal Congress 2010 , Bali, Indonesia. Available at: https://www.geothermal-energy.org/pdf/IGAstandard/WGC/2010/3124.pdf, last accessed 6 October 2017. Guiltinan E., Becker M.W., 2015. Measuring well hydraulic connectivity in fractured bedrock using periodic slug tests, J. Hydrol. , 521, 100– 107. Google Scholar CrossRef Search ADS   Haraden J., 1992. The status of hot dry rock as an energy source, Energy , 17, 777– 786. Google Scholar CrossRef Search ADS   Hollaender F., Hammond P.S., Gringarten A.C., 2002. Harmonic testing for continuous well and reservoir monitoring, in SPE Annual Technical Conference and Exhibition , Society of Petroleum Engineers, San Antonio, Texas. Houwers M.E., Heijnen L.J., Becker A., Rijkers R., 2015. A workflow for the estimation of fault zone permeability for geothermal production: a general model applied on the Roer Valley Graben in the Netherlands, in World Geothermal Congress 2015 , Melbourne, Australia. Available at: https://pangea.stanford.edu/ERE/db/WGC/papers/WGC/2015/12072.pdf, last accessed 6 October 2017. Jacob C.E., Lohman S.W., 1952. Nonsteady flow to a well of constant drawdown in an extensive aquifer, Trans. Am. Geophys. Union , 33, 559– 569. Google Scholar CrossRef Search ADS   López D.L., Smith L., 1996. Fluid flow in fault zones: influence of hydraulic anisotropy and heterogeneity on the fluid flow and heat transfer regime, Water Resour. Res. , 32, 3227– 3235. Google Scholar CrossRef Search ADS   Lyons W.C., 2010. Chapter 4 fluid movement in waterflooded reservoirs, in Working Guide to Reservoir Engineering , pp. 241– 277, Gulf Professional Publishing. Mathias S.A., Butler A.P., 2007. Shape factors for constant-head double-packer permeameters, Water Resour. Res. , 43, W06430, doi:10.1029/2006WR005279. Matthews C.S., Russell D.G., 1967. Pressure Buildup and Flow Tests in Wells , pp. 163, Society of Petroleum Engineers. Mitchell T.M., Faulkner D.R., 2009. The nature and origin of off-fault damage surrounding strike-slip fault zones with a wide range of displacements: a field study from the Atacama fault system, northern Chile. J. Struct. Geol. , 31, 802– 816. Google Scholar CrossRef Search ADS   Ortiz R.A.E., Jung R., Renner J., 2013. Two-dimensional numerical investigations on the termination of bilinear flow in fractures, Solid Earth , 4, 331– 345. Google Scholar CrossRef Search ADS   Rabinovich A., Barrash W., Cardiff M., Hochstetler D.L., Bakhos T., Dagan G., Kitanidis P.K., 2015. Frequency dependent hydraulic properties estimated from oscillatory pumping tests in an unconfined aquifer, J. Hydrol. , 531, 2– 16. Google Scholar CrossRef Search ADS   Rasmussen T.C., Haborak K.G., Young M.H., 2003. Estimating aquifer hydraulic properties using sinusoidal pumping at the Savannah River site, South Carolina, USA, Hydrol. J. , 11, 466– 482. Rehbinder G., 1996. The double packer permeameter with narrow packers. Analytical solution for non steady flow, Appl. Sci. Res. , 56, 255– 279. Google Scholar CrossRef Search ADS   Renard P., Glenz D., Mejias M., 2009. Understanding diagnostic plots for well-test interpretation, Hydrol. J. , 17, 589– 600. Renner J., Messar M., 2006. Periodic pumping tests, Geophys. J. Int. , 167, 479– 493. Google Scholar CrossRef Search ADS   Schindler M., Nami P., Schellschmidt R., Teza D., Tischner T., 2008. Summary of hydraulic stimulation operations in the 5 km deep crystalline HDR/EGS reservoir at Soultz-sous-Forêts, in Thirty-Third Workshop on Geothermal Reservoir Engineering , 28–30 January, Stanford University, Stanford, CA, USA. Available at: https://www.geothermal-energy.org/pdf/IGAstandard/SGW/2008/schindle.pdf, last accessed 6 October 2017. Shaoul J.R., Ross M.J., Spitzer W.J., Wheaton S.R., Mayland P.J., Singh A.P., 2007. Massive hydraulic fracturing unlocks deep tight gas reserves in India, in European Formation Damage Conference , 30 May–1 June, Scheveningen, The Netherlands. doi:10.2118/107337-MS. Sheldon H.A., Micklethwaite S., 2007. Damage and permeability around faults: implications for mineralization, Geology , 35, 903– 906. Google Scholar CrossRef Search ADS   Sibson R.H., Moore J.M.M., Rankin A.H., 1975. Seismic pumping—a hydrothermal fluid transport mechanism, J. geol. Soc. , 131, 635– 659. Google Scholar CrossRef Search ADS   Song I., Renner J., 2007. Analysis of oscillatory fluid flow through rock samples, Geophys. J. Int. , 170, 195– 204. Google Scholar CrossRef Search ADS   Vinci C., Steeb H., Renner J., 2015. The imprint of hydro-mechanics of fractures in periodic pumping tests, Geophys. J. Int. , 202, 1613– 1626. Google Scholar CrossRef Search ADS   Weir G.J., 1999. Single-phase flow regimes in a discrete fracture model, Water Resour. Res. , 35, 65– 73. Google Scholar CrossRef Search ADS   Wibberley C.A.J., Shimamoto T., 2003. Internal structure and permeability of major strike-slip fault zones: the Median Tectonic Line in Mie Prefecture, Southwest Japan, J. Struct. Geol. , 25, 59– 78. Google Scholar CrossRef Search ADS   APPENDIX A: STRUCTURAL INFORMATION FOR INTERVALS Figure A1. View largeDownload slide Images of the borehole wall for the five intervals as gained from logging with an optical camera. For i2 (i5), a section of about 1 m length is also represented by a photograph of the recovered cores. Figure A1. View largeDownload slide Images of the borehole wall for the five intervals as gained from logging with an optical camera. For i2 (i5), a section of about 1 m length is also represented by a photograph of the recovered cores. APPENDIX B: MIRROR-SYMMETRIC DOUBLE-PACKER MODEL For an isotropic and homogenous subsurface, the governing diffusion equation reads in cylindrical coordinates:   \begin{equation}\frac{{{\partial ^2}p(r,z,t)}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial p(r,z,t)}}{{\partial r}} + \frac{{{\partial ^2}p(r,z,t)}}{{\partial {z^2}}} = \frac{1}{D}\frac{{\partial p(r,z,t)}}{{\partial t}}\end{equation} (B1)where r is radial distance from the centre of the borehole and z denotes the vertical direction (upright positive). Employing separation of variables, that is,   \begin{equation*}p(r,z,t) = P(r,z)\Phi (t) = P(r,z)\exp (i\omega t),\end{equation*}the spatial variation of pressure obeys   \begin{equation}\frac{{{\partial ^2}P(r,z)}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial P(r,z)}}{{\partial r}} + \frac{{{\partial ^2}P(r,z)}}{{\partial {z^2}}} = \frac{{i\omega }}{D}P(r,z).\end{equation} (B2) We consider a double-packer probe enclosing an interval with a total height of 2h0 in a well with radius ri. Instead of modeling two packers, we assume mirror symmetry to the horizontal plane bisecting the interval (Fig. B1) and thus we can restrict to treating the upper half of the model space. The locations of the bottom and top of the packer are z = h0 and z = hp, respectively, and the aquifer extends vertically to z = b (Fig. B1). We adapt boundary conditions from Rehbinder (1996) and Mathias & Butler (2007): BC1: $${\displaystyle\frac{{\partial P(r,z)}}{{\partial z}} = 0}\quad {z = b}$$ no-flow condition for vertical bounds of model BC2: $${P(r,z) = 0}\quad {r \to \infty }$$ radially infinite reservoir BC3: $${Q = 4\pi {r_{\rm i}}\int\nolimits_0^{{h_0}} {q({r_{\rm i}},z)} {\rm{d}}z\quad {\rm{ with }}\quad q({r_{\rm i}},z) = - \displaystyle\frac{T}{{2\rho g{h_0}}}{{ {\displaystyle\frac{{\partial P(r,z)}}{{\partial r}}} |}_{{r_{\rm{i}}}}}}\quad {0 \le z} \le {h_0}\quad $$ integrating flux q over the injection interval has to yield the constant flow rate Q BC4: $${P({r_{\rm i}},z) = {p_1}}\quad {{h_{\rm p}} \le z} \le b \quad $$ average pressure above packer (neglecting effect of fluid column) BC5: $${P({r_{\rm i}},z) = {p_2}}\quad {0 \le z} \le {h_0} \quad $$ average pressure in injection interval (neglecting effect of fluid column) BC6: $${P({r_{\rm i}},z) = ({p_2} - {p_1})\displaystyle\frac{{{h_{\rm p}} - z}}{{{h_{\rm p}} - {h_0}}} + {p_1}}\quad {{h_0} \le z} \le {h_{\rm p}}$$ linear pressure change along packer section where the pressures in the well above and in the interval are denoted p1 and p2, respectively. Figure B1. View largeDownload slide Model geometry for the analytical solution of a double-packer test. The model has mirror symmetry with respect to the horizontal plane that bisects the interval. Figure B1. View largeDownload slide Model geometry for the analytical solution of a double-packer test. The model has mirror symmetry with respect to the horizontal plane that bisects the interval. After Fourier cosine transforming of z, that is,   \begin{eqnarray}{{P_n}(r) = \frac{2}{b}\int\nolimits_0^b {P(r,z){\rm {cos}}({a_n}z)\,{\rm{d}}z,}}\quad {{a_n} = \frac{{n\pi }}{b}} \end{eqnarray} (B3)and accounting for BC1, eq. (B2) becomes   \begin{equation}\frac{{{\partial ^2}{P_n}(r)}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial {P_n}(r)}}{{\partial r}} - \left( {{a}_n^{2} + \frac{{i\omega }}{D}} \right){P_n}(r) = 0\end{equation} (B4)with a general solution in the form of   \begin{eqnarray}{{P_n}(r) = {A_n}{K_0}({\eta _n}r) + {B_n}{I_0}({\eta _n}r),}\quad {{\eta _n} = \sqrt {{a}_n^{2} + \frac{{i\omega }}{D}} } \end{eqnarray} (B5)where K0 and I0 denote modified Bessel functions of the first and second kind of zero order. Since I0 → ∞ when r → ∞ (BC2), we have to require Bn = 0, that is,   \begin{equation}{P_n}(r) = {A_n}{K_0}({\eta _n}r).\end{equation} (B6) The backtransform reads   \begin{equation}P\left( {r,z} \right) = \frac{{{A_0}{K_0}({\eta _0}r)}}{2} + \sum\limits_{n = 1}^\infty {{A_n}{K_0}({\eta _n}r){\rm {cos}}\left( {{a_n}z} \right)} .\end{equation} (B7) Application of BC4–BC6 gives   \begin{equation}\frac{{{A_0}{K_0}({\eta _0}{r_{\rm i}})}}{2} + \sum\limits_{n = 1}^\infty {{A_n}{K_0}({\eta _n}{r_{\rm i}}){\rm {cos}}\left( {{a_n}z} \right)} = \left\{ {\begin{array}{c@\quad c} {{p_2}}& {0 \le z \le {h_0}}\\ {({p_2} - {p_1})(\frac{{{h_{\rm p}} - z}}{{{h_{\rm p}} - {h_0}}}) + {p_1}}& {{h_0} \le z \le {h_{\rm p}}}\\ {{p_1}}& {{h_{\rm p}} \le z \le b} \end{array}} \right.\end{equation} (B8)which upon Fourier-cosine transformation of z yields   \begin{eqnarray}&&{\frac{{{A_0}{K_0}({\eta _0}{r_{\rm i}})}}{2}\underbrace {\int\nolimits_0^b {{\rm {cos}}\left( {{a_m}z} \right)}\, {\rm{d}}z}_{ = 0} + \sum\limits_{n = 1}^\infty {{A_n}{K_0}({\eta _n}{r_{\rm i}})\int\nolimits_0^b {{\rm {cos}}\left( {{a_n}z} \right){\rm {cos}}\left( {{a_m}z} \right){\rm{d}}z} }}\nonumber\\ &=& \int\nolimits_0^{{h_0}} {{p_2}{\rm {cos}}\left( {{a_m}z} \right){\rm{d}}z} + \int\nolimits_{{h_0}}^{{h_{\rm p}}} {\left[ {({p_2} - {p_1})\frac{{{h_{\rm p}} - z}}{{{h_{\rm p}} - {h_0}}} + {p_1}} \right]{\rm {cos}}\left( {{a_m}z} \right){\rm{d}}z} + \int\nolimits_{{h_{\rm p}}}^b {{p_1}{\rm {cos}}\left( {{a_m}z} \right){\rm{d}}z} .\end{eqnarray} (B9) For an − am ≠ 0, the second integrand on the left-hand side can be modified as   \begin{equation*}{\rm {cos}}\left( {{a_n}z} \right){\rm {cos}}\left( {{a_m}z} \right) = \frac{{\cos \left[ {({a_n} - {a_m})z} \right] + \cos \left[ {({a_n} + {a_m})z} \right]}}{2}\end{equation*} and thus   \begin{equation*}\int\nolimits_0^b {\displaystyle\frac{{\cos \left[ {({a_n} - {a_m})z} \right] + \cos \left[ {({a_n} + {a_m})z} \right]}}{2}{\rm{d}}z} = \left. {\left[ {\displaystyle\frac{{\sin \left[ {(n - n)\pi \frac{z}{b}} \right]}}{{2({a_n} - {a_m})}} + \frac{{\sin \left[ {(n + n)\pi \frac{z}{b}} \right]}}{{2({a_n} + {a_m})}}} \right]} \right|_0^b = 0,\end{equation*} that is, all elements of the sum vanish but the one for which an − am = 0. In this case, the left-hand side of eq. (B9) can be algebraically manipulated as   \begin{equation}{A_m}{K_0}({\eta _m}{r_{\rm i}})\int\nolimits_0^b {{\rm cos}^2}\left( {{a_m}z} \right){\rm{d}}z = {A_m}{K_0}({\eta _m}{r_{\rm i}})\frac{b}{2}.\end{equation} (B10) After integration, the right-hand side sums up to   \begin{equation}{p_2}\frac{{\cos \left( {{a_m}{h_0}} \right)}}{{\left( {{h_{\rm p}} - {h_0}} \right){a}_m^{2}}} - {p_2}\frac{{\cos \left( {{a_m}{h_{\rm p}}} \right)}}{{\left( {{h_{\rm p}} - {h_0}} \right){a}_m^{2}}} - {p_1}\left[ {\frac{{\cos \left( {{a_m}{h_0}} \right)}}{{\left( {{h_{\rm p}} - {h_0}} \right){a}_m^{2}}} - \frac{{\cos \left( {{a_m}{h_{\rm p}}} \right)}}{{\left( {{h_{\rm p}} - {h_0}} \right){a}_m^{2}}}} \right].\end{equation} (B11) With both, left- and right-hand sides, eqs (B10) and (B11), respectively, we can solve for the unknown coefficients   \begin{equation}{A_m} = \displaystyle\frac{{{p_2}\left[ {\displaystyle\frac{{\cos \left( {{a_m}{h_0}} \right)}}{{\left( {{h_{\rm p}} - {h_0}} \right){a}_m^{2}}} - \displaystyle\frac{{\cos \left( {{a_m}{h_{\rm p}}} \right)}}{{\left( {{h_{\rm p}} - {h_0}} \right){a}_m^{2}}}} \right] - {p_1}\left[ {\displaystyle\frac{{\cos \left( {{a_m}{h_0}} \right)}}{{\left( {{h_{\rm p}} - {h_0}} \right){a}_m^{2}}} - \displaystyle\frac{{\cos \left( {{a_m}{h_{\rm p}}} \right)}}{{\left( {{h_{\rm p}} - {h_0}} \right){a}_m^{2}}}} \right]}}{{\displaystyle\frac{b}{2}{K_0}({\eta _m}{r_{\rm i}})}}.\end{equation} (B12) To derive A0, we apply Darcy's law to eq. (B7) and obtain the flow rate according to (BC3), that is,   \begin{equation}Q = 2\pi {r_{\rm i}}\frac{T}{{\rho g{h_0}}}\left( {\frac{{{A_0}{\eta _0}{K_1}({\eta _0}{r_{\rm i}})}}{2}{h_0} + \sum\limits_{n = 1}^\infty {{A_n}{\eta _n}{K_1}({\eta _n}{r_{\rm i}})\frac{{\sin \left( {{a_n}{h_0}} \right)}}{{{a_n}}}} } \right).\end{equation} (B13)and thus   \begin{equation}{A_0} = \frac{{Q\rho g}}{{\pi {r_{\rm i}}T{\eta _0}{K_1}({\eta _0}{r_{\rm i}})}} - 2\sum\limits_{n = 1}^\infty {{A_n}\frac{{{\eta _n}{K_1}({\eta _n}{r_{\rm i}})}}{{{\eta _0}{K_1}({\eta _0}{r_{\rm i}})}}\frac{{\sin \left( {{a_n}{h_0}} \right)}}{{{a_n}{h_0}}}} .\end{equation} (B14) Finally, the pressure function eq. (B7) reads   \begin{equation}P\left( {r,z} \right) = \frac{{Q\rho g}}{{2\pi {r_{\rm i}}T}}\frac{{{K_0}({\eta _0}r)}}{{{\eta _0}{K_1}({\eta _0}{r_{\rm i}})}} - \frac{{{K_0}({\eta _0}r)}}{{{\eta _0}{K_1}({\eta _0}{r_{\rm i}})}}\sum\limits_{n = 1}^\infty {{A_n}{\eta _n}{K_1}({\eta _n}{r_{\rm i}})\frac{{\sin \left( {{a_n}{h_0}} \right)}}{{{a_n}{h_0}}}} + \sum\limits_{n = 1}^\infty {{A_n}{K_0}({\eta _n}r){\rm {cos}}\left( {{a_n}z} \right)} \end{equation} (B15)where $${A_n} = \displaystyle\frac{{2( {{p_2} - {p_1}} )[ {\cos ( {{a_n}{h_0}} ) - \cos ( {{a_n}{h_{\rm p}}} )} ]}}{{b{K_0}({\eta _n}{r_{\rm i}})( {{h_{\rm p}} - {h_0}} ){a}_n^{2}}}.$$ For the injection interval, with P(ri, z < h0) = p2, we now find a relation between the ratios of flow rate and interval pressure on the one hand and between interval pressure and pressure outside the interval on the other hand:   \begin{eqnarray}{\frac{Q}{{{p_2}}} = \frac{{2\pi {r_{\rm i}}T}}{{\rho g}}\left\{ {\frac{{{\eta _0}{K_1}({\eta _0}{r_{\rm i}})}}{{{K_0}({\eta _0}{r_{\rm i}})}} + 2\frac{{1 - \frac{{{p_1}}}{{{p_2}}}}}{{b\left( {{h_{\rm p}} - {h_0}} \right)}}} {\sum\limits_{n = 1}^\infty {\frac{{\cos \left( {{a_n}{h_0}} \right) - \cos \left( {{a_n}{h_{\rm p}}} \right)}}{{{K_0}({\eta _n}{r_{\rm i}}){a}_n^{2}}}\left[ {\frac{{{\eta _n}{K_1}({\eta _n}{r_{\rm i}})}}{{{h_0}{a_n}}}\sin \left( {{a_n}{h_0}} \right) - \frac{{{\eta _0}{K_0}({\eta _n}{r_{\rm i}}){K_1}({\eta _0}{r_{\rm i}})}}{{{K_0}({\eta _0}{r_{\rm i}})}}\cos\left( {{a_n}z} \right)} \right]} } \right\}} .\end{eqnarray} (B16) It can be seen that solution consists of two parts: standard radial flow for infinite line source and an additional term depending on packer positioning. When the pressure outside of the interval remains unaltered by the pumping operation, that is, p1 = p0 = const., we can—without loss of generality—chose the reference pressure to p0 = 0. In this case, eq. (B16) simplifies to   \begin{eqnarray}\frac{Q}{{{p_2}}} = \frac{{2\pi {r_{\rm i}}T}}{{\rho g}}\left\{ {\frac{{{\eta _0}{K_1}({\eta _0}{r_{\rm i}})}}{{{K_0}({\eta _0}{r_{\rm i}})}} + \frac{2}{{b\left( {{h_{\rm p}} - {h_0}} \right)}}} {\sum\limits_{n = 1}^\infty {\frac{{\cos \left( {{a_n}{h_0}} \right) - \cos \left( {{a_n}{h_{\rm p}}} \right)}}{{{K_0}({\eta _n}{r_{\rm i}}){a}_n^{2}}}\left[ {\frac{{{\eta _n}{K_1}({\eta _n}{r_{\rm i}})}}{{{h_0}{a_n}}}\sin \left( {{a_n}{h_0}} \right) - \frac{{{\eta _0}{K_0}({\eta _n}{r_{\rm i}}){K_1}({\eta _0}{r_{\rm i}})}}{{{K_0}({\eta _0}{r_{\rm i}})}}{\rm {cos}}\left( {{a_n}z} \right)} \right]} } \right\}, \end{eqnarray} (B17)an analytical expression for the ratio of flow rate and interval pressure. APPENDIX C: CONCENTRIC CYLINDRICAL SHELLS SURROUNDING THE BOREHOLE C1 Analytical solution We consider radial flow from a borehole into and beyond a single concentric cylindrical shell located between ri, the radius of the injection well and r1 (Fig. C1). The diffusion equation can be written as:   \begin{equation}\frac{{{\partial ^2}p(r,t)}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial p(r,t)}}{{\partial r}} = \frac{1}{D}\frac{{\partial p(r,t)}}{{\partial t}}.\end{equation} (C1) Figure C1. View largeDownload slide Geometry of radial flow from a well into N concentric shells and the surrounding medium. Figure C1. View largeDownload slide Geometry of radial flow from a well into N concentric shells and the surrounding medium. Here, p is fluid pressure, t is time, r is radial distance and D is hydraulic diffusivity. The hydraulic properties of the shell (subscript 1) and the surrounding medium (subscript 2) are described by permeability k1, 2 = μT1, 2/(ρgh) and diffusivity D1, 2 = k1, 2/(s1, 2μ) = T1, 2/S1, 2 where μ denotes the fluid viscosity, T1, 2, S1, 2 and s1, 2 transmissivity, storativity and specific storage capacity of the two media, respectively, g gravitational acceleration and h the length of the injection interval. Solving the diffusion equation for periodic pressure variations in the well by separation of variables (compare Renner Messar 2006) and account for boundary conditions $${{p_2}(r \ge {r_1},t) = 0}\quad {r \to \infty }\quad $$ radially infinite reservoir $${{p_1}(r,t) = {p_2}(r,t)}\quad {r = {r_1}} \quad $$ continuity of pressure at outer shell radius $${{T_1}\displaystyle\frac{{\partial {p_1}(r,t)}}{{\partial r}} = {T_2}\frac{{\partial {p_2}(r,t)}}{{\partial r}}}\quad {r = {r_{\rm{1}}}}\quad $$ continuity of flow rate at outer shell radius yields amplitude ratio and phase shift between flow rate and interval pressure of   \begin{equation}{\delta _{{\rm{qp}}}} = 2\pi {r_{\rm i}}\frac{{{T_1}}}{{\rho g}}\left| {\frac{{{\eta _1}\left( { - {m_1}{K_1}\left( {{\eta _1}{r_{\rm i}}} \right) + {I_1}\left( {{\eta _1}{r_{\rm i}}} \right)} \right)}}{{{m_1}{K_0}\left( {{\eta _1}{r_{\rm i}}} \right) + {I_0}\left( {{\eta _1}{r_{\rm i}}} \right)}}} \right|,\end{equation} (C2)and   \begin{equation}{\varphi _{{\rm{qp}}}} = {\rm{arg}}\left( {\frac{{{\eta _1}\left( { - {m_1}{K_1}\left( {{\eta _1}{r_{\rm i}}} \right) + {I_1}\left( {{\eta _1}{r_{\rm i}}} \right)} \right)}}{{{m_1}{K_0}\left( {{\eta _1}{r_{\rm i}}} \right) + {I_0}\left( {{\eta _1}{r_{\rm i}}} \right)}}} \right),\end{equation} (C3)respectively. The constant m1 is related to the ratios of permeability (or equivalently transmissivity) and diffusivity according to   \begin{equation}{m_1} = \displaystyle\frac{{{I_1}\left( {{\eta _1}{r_1}} \right){K_0}\left( {{\eta _2}{r_1}} \right) + {I_0}\left( {{\eta _1}{r_1}} \right){K_1}\left( {{\eta _2}{r_1}} \right)\sqrt {\displaystyle\frac{{{D_1}}}{{{D_2}}}} \frac{{{k_2}}}{{{k_1}}}}}{{{K_1}\left( {{\eta _1}{r_1}} \right){K_0}\left( {{\eta _2}{r_1}} \right) - {K_0}\left( {{\eta _1}{r_1}} \right){K_1}\left( {{\eta _2}{r_1}} \right)\sqrt {\displaystyle\frac{{{D_1}}}{{{D_2}}}} \frac{{{k_2}}}{{{k_1}}}}},\end{equation} (C4)where I0 and I1 represent the modified Bessel functions of the first kind of zero and first order and K0 and K1 the modified Bessel functions of the second kind of zero and first order. The arguments of these Bessel functions contain the complex parameters $${\eta _{1,2}} = \sqrt {i\omega /{D_{1,2}}} $$ that depend on angular frequency ω and diffusivity of the two media. The analytical solution for N concentric cylindrical shells (Fig. C1) is given as   \begin{equation}{\delta _{{\rm{qp}}}} = 2\pi {r_{\rm i}}\frac{{{T_1}}}{{\rho g}}\left| {\frac{{{\eta _1}\left( { - {c_1}{K_1}\left( {{\eta _1}{r_{\rm i}}} \right) + {I_1}\left( {{\eta _1}{r_{\rm i}}} \right)} \right)}}{{{c_1}{K_0}\left( {{\eta _1}{r_{\rm i}}} \right) + {I_0}\left( {{\eta _1}{r_{\rm i}}} \right)}}} \right|,\end{equation} (C5)and   \begin{equation}{\varphi _{{\rm{qp}}}} = {\rm{arg}}\left( {\frac{{{\eta _1}\left( { - {c_1}{K_1}\left( {{\eta _1}{r_{\rm i}}} \right) + {I_1}\left( {{\eta _1}{r_{\rm i}}} \right)} \right)}}{{{c_1}{K_0}\left( {{\eta _1}{r_{\rm i}}} \right) + {I_0}\left( {{\eta _1}{r_{\rm i}}} \right)}}} \right),\end{equation} (C6)and j = 1, …N−1. For N = 1, c1 = m1 but for N > 1, to obtain c1 one has to successively determine c2, c3… and cN from   \begin{equation*}{c_N} = - \displaystyle\frac{{\frac{{{k_{N + 1}}}}{{{k_N}}}\sqrt {\displaystyle\frac{{{D_N}}}{{{D_{N + 1}}}}} {K_1}({\eta _{N + 1}}{r_N}){I_0}({\eta _N}{r_N}) + {I_1}({\eta _N}{r_N}){K_0}({\eta _{N + 1}}{r_N})}}{{\displaystyle\frac{{{k_{N + 1}}}}{{{k_N}}}\sqrt {\frac{{{D_N}}}{{{D_{N + 1}}}}} {K_1}({\eta _{N + 1}}{r_N}){K_0}({\eta _N}{r_N}) - {K_1}({\eta _N}{r_N}){K_0}({\eta _{N + 1}}{r_N})}}\end{equation*}   \begin{equation*}{c_j} = - \displaystyle\frac{{{I_0}({\eta _j}{r_j})\displaystyle\frac{{{k_{j + 1}}}}{{{k_j}}}\sqrt {\frac{{{D_j}}}{{{D_{j + 1}}}}} \left( {{c_{j + 1}}{K_1}({\eta _{j + 1}}{r_j}) - {I_1}({\eta _{j + 1}}{r_j})} \right) + {I_1}({\eta _j}{r_j})\left( {{c_{j + 1}}{K_0}({\eta _{j + 1}}{r_j}) + {I_0}({\eta _{j + 1}}{r_j})} \right)}}{{{K_0}({\eta _j}{r_j})\displaystyle\frac{{{k_{j + 1}}}}{{{k_j}}}\sqrt {\displaystyle\frac{{{D_j}}}{{{D_{j + 1}}}}} \left( {{c_{j + 1}}{K_1}({\eta _{j + 1}}{r_j}) - {I_1}({\eta _{j + 1}}{r_j})} \right) - {K_1}({\eta _j}{r_j})\left( {{c_{j + 1}}{K_0}({\eta _{j + 1}}{r_j}) + {I_0}({\eta _{j + 1}}{r_j})} \right)}}\end{equation*} C2 Parameter study For the parameter study, we restrict to the simple case of a single shell. We investigated the sensitivity of phase shift to the various model parameters by systematic parameter variation. Phase shift is reported as a function of $${r_{\rm i}}\sqrt {\omega /{D_{{\rm{shell}}}}} $$ where Dshell = D1 denotes the hydraulic diffusivity of the shell (Fig. C1). This expression can be understood either as a dimensionless frequency or an inverse dimensionless penetration depth given in multiples of the borehole radius. A sigmoidal shape of the relation between phase shift and dimensionless frequency is characteristic for a homogeneous medium. Small phase shifts occur at low frequencies and vice versa with a rather steep switch between the two limits around a dimensionless frequency of 1. This shape is significantly perturbed when ratios between the parameters of shell and medium between one-tenth and a factor of 10 are allowed. The high-frequency limit of 1/8 for phase shift is not strictly valid for the compound model but values as high as 1/5 are found. The phase shift for a shell with a low storage property compared to the medium exhibits a pronounced maximum that forms a shoulder on the steep increase characteristic for the homogeneous medium (Fig. C2a). The location of the maximum moves towards lower frequencies with increasing radius of the shell. A variation of the ratio of storage parameters for shell and medium for a fixed shell radius results in qualitatively similar phase-shift curves as the variation in shell radius (Fig. C2b). Yet, the maximum in phase transforms into a minimum preceded by a new maximum at lower frequency when the storage capacity in the shell switches from lower to higher than the storage capacity in the medium. Figure C2. View largeDownload slide Results of the shell model given with phase shift as a function of dimensionless frequency (a) for various radius ratios and indicated fixed permeability (transmissivity) and diffusivity ratios, (b) for various diffusivity ratios and indicated fixed permeability (transmissivity) and radius ratios, (c) for various permeability (transmissivity) ratios when storativity is covaried such that the diffusivity of shell and medium are identical and (d) for various identical diffusivity and permeability (transmissivity) ratios that correspond to an identical storativity for shell and medium at indicated fixed radius ratio. For identical permeability (transmissivity) of shell and medium, the diffusivity ratio is related to an inverse ratio of specific storage capacity (storativity), that is, D1/D2 = s2/s1 = S2/S1. Figure C2. View largeDownload slide Results of the shell model given with phase shift as a function of dimensionless frequency (a) for various radius ratios and indicated fixed permeability (transmissivity) and diffusivity ratios, (b) for various diffusivity ratios and indicated fixed permeability (transmissivity) and radius ratios, (c) for various permeability (transmissivity) ratios when storativity is covaried such that the diffusivity of shell and medium are identical and (d) for various identical diffusivity and permeability (transmissivity) ratios that correspond to an identical storativity for shell and medium at indicated fixed radius ratio. For identical permeability (transmissivity) of shell and medium, the diffusivity ratio is related to an inverse ratio of specific storage capacity (storativity), that is, D1/D2 = s2/s1 = S2/S1. When permeability (transmissivity) and specific storage capacity (storativity) are modified in the same way, that is, diffusivity is kept constant, the transition from the low- to the high-frequency regime is significantly perturbed. Low permeability and storativity in the shell relative to the surrounding medium leads to a transition that resembles a Heavyside function (Fig. C2c). High permeability and storativity in the shell gives rise to a pronounced maximum in phase shift at a dimensionless frequency of about 0.1 that reaches values as large as 1/5, that is, well above the limit for a homogeneous medium of 1/8. Covarying peremability and diffusivity and keeping the storage parameter uniform for shell and medium has qualitatively similar effects (Fig. C2d). APPENDIX D: BILINEAR FLOW MODEL The model consists of a slit with distinct hydraulic properties embedded in a homogeneous matrix and assumes that the direct volume flow from the source into the matrix along the x-direction is negligible (Fig. D1), but the fracture boundary is considered a harmonic source for flow into the matrix. The governing equations in the matrix and the fracture correspond to diffusion equations:   \begin{equation}\frac{{\partial {p_{\rm m}}(x,y,t)}}{{\partial t}} = {D_{\rm m}}{\nabla ^2}{p_{\rm m}}(x,y,t)\end{equation} (D1)  \begin{equation}\frac{{\partial {p_{\rm f}}(x,t)}}{{\partial t}} = {D_{\rm f}}\frac{{{\partial ^2}{p_{\rm f}}(x,t)}}{{\partial {x^2}}} + \frac{{2{D_{\rm m}}}}{\delta }{\left. {\frac{{\partial {p_{\rm m}}(x,y,t)}}{{\partial y}}} \right|_{\left| y \right| = \delta /2}}\end{equation} (D2)where pm(x, y, t) is the pressure in the matrix, pf(x, t) the pressure in the fracture, δ the width of the fracture and Dm and Df the diffusivity of matrix and fracture, respectively. Figure D1. View largeDownload slide Model for bilinear flow: the fracture boundary is a source for flow into the matrix in y-direction; no direct volume flow from the source into the matrix along the x-direction. Figure D1. View largeDownload slide Model for bilinear flow: the fracture boundary is a source for flow into the matrix in y-direction; no direct volume flow from the source into the matrix along the x-direction. As the direct volume flow from the source into the matrix along the x-direction is negligible, ∂2pm(x, y, t)/∂x2 ≈ 0, so that eq. (D1) reduces to   \begin{equation}\displaystyle\frac{{\partial {p_{\rm m}}(x,y,t)}}{{\partial t}} = {D_{\rm m}}\displaystyle\frac{{{\partial ^2}{p_{\rm m}}(x,y,t)}}{{\partial {y^2}}}.\end{equation} (D3) The pressure in the matrix, that is, the solution of eq. (D1), is given by the standard solution for diffusion from a harmonic source into a homogeneous semi-infinite medium (Carslaw & Jaeger 1986):   \begin{equation}{p_{\rm m}}\left( {x,y,t} \right) = \displaystyle\frac{{\left| y \right| - \delta /2}}{{2\sqrt {\pi {D_{\rm m}}} }}\mathop \int \nolimits_0^t {p_{\rm f}}\left( {x,\lambda } \right)\displaystyle\frac{{\exp \left[ { - \displaystyle\frac{{{{(\left| y \right| - \delta /2)}^2}}}{{4{D_{\rm m}}\left( {t - \lambda } \right)}}} \right]}}{{{{\left( {t - \lambda } \right)}^{3/2}}}}{\rm{d}}\lambda \end{equation} (D4)whose insertion into eq. (D2) yields   \begin{equation}{D_{\rm f}}\displaystyle\frac{{{\partial ^2}{p_{\rm f}}\left( {x,t} \right)}}{{\partial {x^2}}} + \displaystyle\frac{{2{D_{\rm m}}}}{\delta }\displaystyle\frac{\partial }{{\partial y}}{\left[ {\displaystyle\frac{{\left| y \right| - \delta /2}}{{2\sqrt {\pi {D_{\rm m}}} }}\mathop \int \nolimits_0^t {p_{\rm f}}\left( {x,\lambda } \right)\displaystyle\frac{{\exp \left[ { - \displaystyle\frac{{{{(\left| y \right| - \delta /2)}^2}}}{{4{D_{\rm m}}\left( {t - \lambda } \right)}}} \right]}}{{{{\left( {t - \lambda } \right)}^{3/2}}}}{\rm{d}}\lambda } \right]_{\left| y \right| - \delta /2 = 0}} = \displaystyle\frac{{\partial {p_{\rm f}}\left( {x,t} \right)}}{{\partial t}}.\end{equation} (D5) Differentiation of the second term by part restricting to the section y ≥ 0 because of the symmetry of the problem gives   \begin{equation}{D_{\rm f}}\displaystyle\frac{{{\partial ^2}{p_{\rm f}}(x,t)}}{{\partial {x^2}}} + \displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\sqrt \pi \delta }}{\left[ {\int\nolimits_0^t {{p_{\rm f}}(x,\lambda )\left( {1 + \displaystyle\frac{{{{(y - {\delta / 2})}^2}}}{{2{D_{\rm m}}(t - \lambda )}}} \right)\displaystyle\frac{{\exp \left[ { - \displaystyle\frac{{{{(\left| y \right| - {\delta / 2})}^2}}}{{4{D_{\rm m}}(t - \lambda )}}} \right]}}{{{{\left( {t - \lambda } \right)}^{3/2}}}}} {\rm{d}}\lambda } \right]_{y = {\delta / 2}}} = \displaystyle\frac{{\partial {p_{\rm f}}(x,t)}}{{\partial t}}.\end{equation} (D6) Applying y = δ/2 we get   \begin{equation}{D_{\rm f}}\displaystyle\frac{{{\partial ^2}{p_{\rm f}}(x,t)}}{{\partial {x^2}}} + \displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\sqrt \pi \delta }}\int\nolimits_0^t {\displaystyle\frac{{{p_{\rm f}}(x,\lambda )}}{{{{\left( {t - \lambda } \right)}^{\frac{3}{2}}}}}} {\rm{d}}\lambda = \displaystyle\frac{{\partial {p_{\rm f}}(x,t)}}{{\partial t}}\end{equation} (D7)and thus for harmonic pressure variations in the fracture   \begin{equation}{D_{\rm f}}\displaystyle\frac{{{\partial ^2}{p_{\rm f}}(x)}}{{\partial {x^2}}}A\exp (i\omega t) + \displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\sqrt \pi \delta }}{p_{\rm f}}(x)\int\nolimits_0^t {\displaystyle\frac{{A\exp (i\omega \lambda )}}{{{{\left( {t - \lambda } \right)}^{\frac{3}{2}}}}}} {\rm{d}}\lambda = i\omega A\exp (i\omega t){p_{\rm f}}(x).\end{equation} (D8) The integral solution is   \begin{equation}\int\nolimits_0^t {\displaystyle\frac{{A\exp (i\omega \lambda )}}{{{{\left( {t - \lambda } \right)}^{\frac{3}{2}}}}}{\rm{d}}\lambda } = \gamma \left( { - 0.5,i\omega t} \right)\exp (i\omega t) = \left[ {\Gamma \left( { - 0.5} \right) - \Gamma \left( { - 0.5,i\omega t} \right)} \right]\sqrt {i\omega } \exp (i\omega t)\end{equation} (D9)where γ( − 0.5, iωt), Γ( − 0.5) and Γ( − 0.5, iωt) represent the lower incomplete gamma function, the gamma function and the upper incomplete gamma function, respectively. For times longer than 1.25 times the period, the deviation of γ( − 0.5, iωt) from Γ( − 0.5) is less than 1 per cent and beyond 2.2 times the period the deviation falls below 0.5 per cent. Thus, we can approximate the integral as   \begin{equation}\mathop \int \nolimits_0^t \displaystyle\frac{{\exp \left( {i\omega \lambda } \right)}}{{{{\left( {t - \lambda } \right)}^{3/2}}}}{\rm{d}}\lambda \cong {\rm{\Gamma }}\left( { - 0.5} \right)\sqrt {i\omega } \exp \left( {i\omega t} \right) = - 2\sqrt {i\pi \omega } \exp \left( {i\omega t} \right)\end{equation} (D10)to arrive at   \begin{equation}{D_{\rm f}}\displaystyle\frac{{{\partial ^2}{p_{\rm f}}\left( x \right)}}{{\partial {x^2}}} - \displaystyle\frac{2}{\delta }\sqrt {i{D_{\rm m}}\omega } {p_{\rm f}}\left( x \right) = i\omega {p_{\rm f}}\left( x \right).\end{equation} (D11)when canceling the harmonic time dependence. Eq. (D11) is a homogeneous second-order differential equation with the solution   \begin{equation}{p_{\rm f}}\left( x \right) = c\exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } {\rm{ + }}i\omega } \right)} \right]}^{1/2}}x} \right\} + {\rm{d}}H\left( L \right)\end{equation} (D12)subject to the boundary conditions pf(x)→ 0 as x → ∞ (infinite fracture), or pf(x) = p0 as x = L (finite fracture with constant-pressure boundary), or $${{\partial {p_{\rm f}}( x )} / {\partial x}} = 0$$ as x = L (finite fracture with no-flow boundary) where c is a constant, H(L) is a parameter defined by   \begin{equation} H\left( L \right) = \left\{ {\begin{array} {l@\quad l} 0 & L \to \infty\\ 1 & L\ {\rm{finite}} \end{array}} \right.\end{equation} (D13)and   \begin{equation}d = \left\{ {\begin{array}{l} {{p_0} - c\exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i\omega } \right)} \right]}^{1/2}}L} \right\}{\rm{ for\,\, constant \hbox{-} pressure\,\, boundary}}}\\\\ {{{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i\omega } \right)} \right]}^{1/2}}cx\exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i\omega } \right)} \right]}^{1/2}}L} \right\}\,{\rm{ +\, const.}}\,\,{\rm{for\,\, no \hbox{-} flow\,\, boundary }}} \end{array}} \right.\end{equation} (D14) To simplify the calculation, we assume the constant appearing for the no-flow boundary to equal zero without loss of generality because it represents a constant shift in pressure that can be accommodated by scale shifting. Using Darcy's law, we get the flow rate as   \begin{eqnarray} {Q_{\rm f}}\left( {x,t} \right) = 2\pi {r_{\rm i}}hq\left( {x,t} \right) &=& - 2\pi {r_{\rm i}}\displaystyle\frac{T}{{\rho g}}\left( {\displaystyle\frac{{\partial {p_{\rm f}}\left( {x,t} \right)}}{{\partial x}}} \right)\nonumber\\ &=& 2\pi {r_{\rm i}}\displaystyle\frac{{TcA}}{{\rho g}}{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i\omega } \right)} \right]^{1/2}}\nonumber\\ && \cdot \left( {\exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } {\rm{ + }}i\omega } \right)} \right]}^{1/2}}x} \right\} - {H_{{\rm{bc}}}}\exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i\omega } \right)} \right]}^{1/2}}L} \right\}} \right){\rm{exp}}\left( {i\omega t} \right) \end{eqnarray} (D15)where T is transmissivity, ρ is fluid density and h is the height of the injection zone. We introduce a parameter reflecting the boundary condition at the fracture tip as   \begin{equation} {H_{{\rm{bc}}}} = \left\{ {\begin{array}{l} {0\ \ \ {\rm{constant\,\, pressure\,\, boundary }}\ }\\\\ {1\ \ \ {\rm{no\,\, flow\,\, boundary }}\ } \end{array}} \right.. \end{equation} (D16) Finally, the complex injectivity results to   \begin{eqnarray} \displaystyle\frac{{{Q_{\rm f}}\left( {x,t} \right)}}{{{p_{\rm f}}\left( {x,t} \right)}} &=& 2\pi {r_{\rm i}}\displaystyle\frac{{Tc}}{{\rho g}}{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i{\rm{\omega }}} \right)} \right]^{1/2}}\nonumber\\ &&{\rm{ }} \cdot \displaystyle\frac{{\exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i\omega } \right)} \right]}^{1/2}}x} \right\} - {H_{{\rm{bc}}}}\exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i{\rm{\omega }}} \right)} \right]}^{1/2}}L} \right\}}}{{c\exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i\omega } \right)} \right]}^{1/2}}x} \right\} + {\rm{d}}H\left( L \right)}}. \end{eqnarray} (D17) Specifically, for an infinite fracture, we find   \begin{equation}{\left. {\displaystyle\frac{{{Q_{\rm f}}(x = 0,t)}}{{{p_{\rm f}}(x = 0,t)}}} \right|_{H\left( L \right) = 0}} = 2\pi {r_{\rm i}}\displaystyle\frac{T}{{\rho g}}{\left( {\displaystyle\frac{{2\sqrt {i\pi {D_{\rm m}}\omega } + i\omega \sqrt {{D_{\rm m}}} \ }}{{\sqrt \pi \delta {D_{\rm f}}}}} \right)^{1/2}}\end{equation} (D18)with amplitude ratio and phase shift of   \begin{equation}{\delta _{{\rm{qp}}}} = 2\pi {r_{\rm i}}\displaystyle\frac{T}{{\rho g}}{\left( {\displaystyle\frac{{2\sqrt {\pi {D_{\rm m}}\omega } + \omega \sqrt {{D_{\rm m}}} \ }}{{\delta {D_{\rm f}}}}} \right)^{1/2}}\end{equation} (D19)and   \begin{equation}{\varphi _{{\rm{qp}}}} = \arg \left( {\sqrt {2\sqrt {i\pi \omega } + i\omega } } \right),\end{equation} (D20)respectively, at the injection point x = 0 (Fig. D2a). Figure D2. View largeDownload slide Phase shift for (a) constant-pressure boundary and (b) no-flow boundary at the tip of a fracture with finite length L and width δ = 0.096 m. (I) Dm = 1 m2 s−1, Df = 10 m2 s−1 and L = 20δ, (II) Dm = 1 m2 s−1, Df = 100 m2 s−1 and L = 20δ, (III) Dm = 1 m2 s−1, Df = 100 m2 s−1 and L = 100δ and (IV) Dm = 1 m2 s−1, Df = 10 m2 s−1 and L = 4δ. The relation for an infinite fracture is shown for comparison in (a). Figure D2. View largeDownload slide Phase shift for (a) constant-pressure boundary and (b) no-flow boundary at the tip of a fracture with finite length L and width δ = 0.096 m. (I) Dm = 1 m2 s−1, Df = 10 m2 s−1 and L = 20δ, (II) Dm = 1 m2 s−1, Df = 100 m2 s−1 and L = 20δ, (III) Dm = 1 m2 s−1, Df = 100 m2 s−1 and L = 100δ and (IV) Dm = 1 m2 s−1, Df = 10 m2 s−1 and L = 4δ. The relation for an infinite fracture is shown for comparison in (a). For a fracture with finite length and a constant-pressure boundary at its tip and assuming p0 = 0, we find   \begin{equation}\displaystyle\frac{{{Q_{\rm f}}(x = 0,t)}}{{{p_{\rm f}}(x = 0,t)}} = \displaystyle\frac{{{{\left. {\displaystyle\frac{{{Q_{\rm f}}(x = 0,t)}}{{{p_{\rm f}}(x = 0,t)}}} \right|}_{H\left( L \right) = 0}}}}{{1 - \exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i\omega } \right)} \right]}^{1/2}}L} \right\}}}\end{equation} (D21)at x = 0 (Fig. D2a). Analogously, we find   \begin{equation}\displaystyle\frac{{{Q_{\rm f}}(x = 0,t)}}{{{p_{\rm f}}(x = 0,t)}} = {\left. {\displaystyle\frac{{{Q_{\rm f}}(x = 0,t)}}{{{p_{\rm f}}(x = 0,t)}}} \right|_{H\left( L \right) = 0}}\left( {1 - \exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i{\rm{\omega }}} \right)} \right]}^{1/2}}L} \right\}} \right)\end{equation} (D22)at x = 0 (Fig. D2b). © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

Exploratory use of periodic pumping tests for hydraulic characterization of faults

Loading next page...
 
/lp/ou_press/exploratory-use-of-periodic-pumping-tests-for-hydraulic-R3T8ytRRXl
Publisher
The Royal Astronomical Society
Copyright
© The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx390
Publisher site
See Article on Publisher Site

Abstract

SUMMARY Periodic pumping tests were conducted using a double-packer probe placed at four different depth levels in borehole GDP-1 at Grimselpass, Central Swiss Alps, penetrating a hydrothermally active fault. The tests had the general objective to explore the potential of periodic testing for hydraulic characterization of faults, representing inherently complex heterogeneous hydraulic features that pose problems for conventional approaches. Site selection reflects the specific question regarding the value of this test type for quality control of hydraulic stimulations of potential geothermal reservoirs. The performed evaluation of amplitude ratio and phase shift between pressure and flow rate in the pumping interval employed analytical solutions for various flow regimes. In addition to the previously presented 1-D and radial-flow models, we extended the one for radial flow in a system of concentric shells with varying hydraulic properties and newly developed one for bilinear flow. In addition to these injectivity analyses, we pursued a vertical-interference analysis resting on observed amplitude ratio and phase shift between the periodic pressure signals above or below packers and in the interval by numerical modeling of the non-radial-flow situation. When relying on the same model the order of magnitude of transmissivity values derived from the analyses of periodic tests agrees with that gained from conventional hydraulic tests. The field campaign confirmed several advantages of the periodic testing, for example, reduced constraints on testing time relative to conventional tests since a periodic signal can easily be separated from changing background pressure by detrending and Fourier transformation. The discrepancies between aspects of the results from the periodic tests and the predictions of the considered simplified models indicate a hydraulically complex subsurface at the drill site that exhibits also hydromechanical features in accord with structural information gained from logging. The exploratory modeling of vertical injectivity shows its potential for analysing hydraulic anisotropy. Yet, more comprehensive modeling will be required to take full advantage of all the pressure records typically acquired when using a double-packer probe for periodic tests. Fracture and flow, Hydrothermal systems, Europe, Fourier analysis, Numerical modelling, Mechanics, theory, and modelling 1 INTRODUCTION Analyses of pressure and/or flow transients associated with pumping operations in wells, that is, injection or production, constitute the primary means for hydraulic characterization of the subsurface (e.g. Fetter 2001). Effective hydraulic properties are determined relying on models considering specific flow geometries and boundary conditions. Specific aspects of observations may guide the selection of one out of a number of alternative models (e.g. Matthews & Russell 1967; Bourdet et al.1989). Periodic pumping tests involve consecutive periods of injection and/or production resulting in alternating flow and pressure and promise several practical and analytical benefits. They may be seen as either an alternative or a complement to conventional testing for constraining effective hydraulic properties (e.g. Rasmussen et al.2003; Renner & Messar 2006), diagnosing the flow regime (e.g. Hollaender et al.2002) and providing valuable information on subsurface heterogeneity (Renner & Messar 2006; Ahn & Horne 2010; Fokker et al.2012; Cardiff et al.2013). The periodic signals can be easily separated from noise (Renner & Messar 2006; Bakhos et al.2014) or background variations associated with long-term operation (e.g. Guiltinan & Becker 2015). In addition, numerical modeling of periodic tests in the frequency domain is much faster than that of typical transients in the time domain (Cardiff et al.2013; Fokker et al.2013). As is the case for conventional testing, effective (or equivalent) hydraulic parameters may be derived from periodic pressure and flow-rate observations by relying on a simple diffusion equation for a homogeneous and isotropic medium in the first step of analyses. Yet, their variation with pumping period can be attributed to—among other possible reasons—spatial heterogeneity (e.g. Renner & Messar 2006; Rabinovich et al.2015). Renner & Messar (2006) performed comprehensive field tests and found that diffusivity values estimated from the analysis of flow-rate and pressure responses at the pumping well increased with increasing period. In contrast, the diffusivity gained from analysis of pressure responses at the pumping and a monitoring well decreased with increasing period, as also found in field tests by Becker & Guiltinan (2010), Rabinovich et al. (2015) and Guiltinan & Becker (2015). The estimated diffusivity values of the two analyses agree at a critical period corresponding to a penetration depth on the order of the well distances. A similar phenomenon was observed in oscillatory hydraulic tests on rock samples (Song & Renner 2007). Fokker et al. (2013) successfully modeled the field observations using a heterogeneous subsurface suggesting that the disparity of hydraulic properties below the critical period reflects effective hydraulic properties of different parts of the subsurface probed depending on period. The work described here aimed at assessing the use of periodic testing for in situ hydraulic characterization of fault structures on the scale of metres to decametres which is of considerable interest in fundamental science and industrial applications. Previous fundamental research covered the relation between fault zone architecture as constrained by structure investigations (e.g. Wibberley & Shimamoto 2003), hydraulic properties (Caine et al.1996; López & Smith 1996) in cases derived from laboratory experiments (e.g. Mitchell & Faulkner 2009), mineralization (Sheldon & Micklethwaite 2007) and earthquake-source mechanics (Sibson et al.1975; Bizzarri 2012). From a perspective of industrial applications, the hydraulic characteristics of faults play an important role for assessing the performance of a reservoir in general (e.g. Fisher & Knipe 2001) but also for stimulations. The term hydraulic stimulation refers to an operation in wells initially developed by the hydrocarbon industry (Economides & Nolte 2000; Shaoul et al.2007) but now also used for developing deep geothermal reservoirs (Baria et al.1999; Schindler et al.2008; Genter et al.2010; Houwers et al.2015; Garcia et al.2016). We present and analyse results from periodic pumping tests in a well penetrating a strike-slip fault considered an analogue for a heat exchanger as needed for economic recovery of petrothermal energy. Our work addresses the need to develop testing protocols and evaluation methods suitable for heterogeneous reservoirs, such as faulted rocks. Equivalent hydraulic parameters are estimated by comparing the relation between flow rate and pressure with a range of previously presented and newly derived analytical solutions. The pressure record in the borehole section separated by packers is compared to that above or below the probe and the potential of the found amplitude ratios and phase shifts for the derivation of hydraulic parameters is explored relying on a numerical model. 2 EXPERIMENTAL DETAILS AND DATA ANALYSIS We use the term injectivity analysis when analysis involves flow rate and pressure in the actively tested section of the pumping well, addressed as interval. Advanced double-packer tests provide the pressure responses above and below the interval, too. We refer to the comparison of these various pressures measured in a single well as vertical-interference test following the convention of addressing the analyses of pressure responses in monitoring wells as interference tests (Matthews & Russell 1967). 2.1 Pumping tests 2.1.1 Location The close to E–W striking Grimsel Breccia Fault zone (GBFZ, Fig. 1) crops out at Grimselpass, Central Swiss Alps, with an altitude of 2160 m above sea level. To the north and the south, two shear zones subparallel to GBFZ, namely Grimsel Shear Zone (GSZ 2) and the Southern Boundary Shear Zone, crop out at distances of about 80 and 500 m, respectively (Fig. 1). The fault zone is hydrothermally active as evidenced by a 250 m wide zone of hydrothermal discharge in the Transitgas AG Tunnel that runs in NS direction at 200 m depth. The origin of the warm fluid is estimated at 4 ± 1 km depth with a reservoir temperature of at least 100 °C (Belgrano et al.2016). Borehole GDP-1 (coordinates of 669΄464/157΄025) with a radius of 0.048 m was drilled about 100 m west of the freshwater reservoir Totensee with a deviation of 24° from the vertical and an azimuth of 164° to intersect the GBFZ approximately normal to its strike. The Transitgas AG tunnel passes by the borehole at about 180 m distance to the west. The deviated borehole was drilled to a final depth of 125.3 m using a thixotropic polymer fluid. The true vertical depth is about 114.5 m. The borehole is open except for the casing stabilizing the first 2.1 m at the top of the well (Fig. 2). Figure 1. View largeDownload slide Sketch of the borehole site. The trajectory of borehole GDP-1 is indicated by the orange line. The depth of the Transitgas AG Tunnel (trajectory in light blue) is 200 m; the zone of inflow is indicated in dark blue. Estimated surface traces of shear zones and faults are represented by grey and dashed black lines, respectively. GBFZ: Grimsel Breccia Fault zone. GSZ 2: Grimsel Shear Zone. SBSZ: Southern Boundary Shear Zone. Figure 1. View largeDownload slide Sketch of the borehole site. The trajectory of borehole GDP-1 is indicated by the orange line. The depth of the Transitgas AG Tunnel (trajectory in light blue) is 200 m; the zone of inflow is indicated in dark blue. Estimated surface traces of shear zones and faults are represented by grey and dashed black lines, respectively. GBFZ: Grimsel Breccia Fault zone. GSZ 2: Grimsel Shear Zone. SBSZ: Southern Boundary Shear Zone. Figure 2. View largeDownload slide Sketch of the measurement setup. The flow rate recorded by the flow meter at the surface is denoted as Qsurface. Pressure recorded above the upper packer, in the interval and below the lower packer are represented as pabo, $${p_{{\mathop{\rm int}} }}$$ and pbel, respectively. Figure 2. View largeDownload slide Sketch of the measurement setup. The flow rate recorded by the flow meter at the surface is denoted as Qsurface. Pressure recorded above the upper packer, in the interval and below the lower packer are represented as pabo, $${p_{{\mathop{\rm int}} }}$$ and pbel, respectively. 2.1.2 Equipment and procedures The test system (Fig. 2) consisted of a reservoir filled with water from Totensee, a rotary pump and several magnetic-inductive flow meters covering a range from 0.35 to 400 L min−1 (accuracy better than 0.5 per cent of full scale) at the surface. A variable-length double-packer probe with downhole pressure gauges (full scale 2 MPa, accuracy 0.1 per cent of full scale) was lowered into the borehole. The packers with 1 m length were inflated from their initial diameter of 72 mm using nitrogen with pressures of up to 2.5 MPa. Data were digitally recorded with a resolution of 16 bit and a time step of 1 or 2 s. Four depth levels were selected for the testing relying on the preliminary results of borehole logging and inspection of drilled cores (Table 1 and Appendix  A). Tests with an interval length of 9.4 m at all four depth levels (i1–i4) were complemented by one test with a short interval of 1.7 m (i5) as part of i2. The latter isolated the prominent fault at a depth of ∼105 m with an orientation of 150/85 estimated from an optical log and core analyses (Appendix  A). The orientation of this fracture is actually also representative of the rather steeply dipping fractures striking subparallel to GBFZ that dominate in interval i1 and all of i2 (Table 1). Structural characterization of the main fault (i3) was not possible owing to extensive borehole breakout and core loss. The shallowest interval at ∼45 m (i4) exhibits a slightly higher fracture density than i2 and i3 mainly resulting from two perpendicular sets of E–W striking fractures, of which one is steeply dipping. Conventional tests (e.g. constant head, constant flow, slug, pulse and recovery tests) were performed at all intervals, while periodic pumping tests were conducted at all depth sections but i4. Table 1. Specifics of the test intervals. Int.  Packer position  Depth at centre of probe  Fracture density  Comment    Lower (m)  Upper (m)  (m)  (1 m−1)    i1  121.8  112.4  117.1  7.6  Several small breccia  i2  111.4  102.0  106.7  5.5  Prominent, isolated fault  i3  86.6  77.2  81.9  –*  Main fault with breccia  i4  49.6  40.2  44.9  8.8  Several fault strands  i5  105.8  104.1  105.0  6.4  Prominent, isolated fault  Int.  Packer position  Depth at centre of probe  Fracture density  Comment    Lower (m)  Upper (m)  (m)  (1 m−1)    i1  121.8  112.4  117.1  7.6  Several small breccia  i2  111.4  102.0  106.7  5.5  Prominent, isolated fault  i3  86.6  77.2  81.9  –*  Main fault with breccia  i4  49.6  40.2  44.9  8.8  Several fault strands  i5  105.8  104.1  105.0  6.4  Prominent, isolated fault  *Severe breakouts and core loss, neither logging (see Appendix  A) nor core analysis revealed structural information. View Large Initial selection of oscillation period and flow rate was based on forward calculations assuming radial flow in a homogeneous medium with diffusivity values considered typical for a fractured formation. The programme was then adjusted on the basis of the first results such that ultimately tests were performed at periods τ between 60 and 1080 s (Table 2). A single period consisted of two injection phases, each with a duration of half of the period, with nominal flow rates that differed by about 3 L min−1. By varying the lower rate in a period between 0 and 7 L min−1, we also modified the average pressure at which periodic testing took place (Fig. 3). According to the simple scaling relation $${z_p} \sim \sqrt {D\tau } $$, the chosen periods correspond to penetration depths zp of several decimetres to tens of metres for assumed values of hydraulic diffusivityD between 1E–5 and 1E–2 m2 s−1. Thus, neither Totensee nor the tunnel are expected to constitute relevant boundary conditions for the exerted flow. Figure 3. View largeDownload slide Pressures (left y-axis) and flow rate (green, right y-axis) recorded during the periodic testing at a depth of about 105 m (i5) for oscillation periods of (a) 60 s (i5–1, i5–2, see Table 2) and (b) 180 and 540 s (i5–3, i5–4). Figure 3. View largeDownload slide Pressures (left y-axis) and flow rate (green, right y-axis) recorded during the periodic testing at a depth of about 105 m (i5) for oscillation periods of (a) 60 s (i5–1, i5–2, see Table 2) and (b) 180 and 540 s (i5–3, i5–4). Table 2. Results of the Fourier analyses of recorded pressures and flow rate. Amplitude ratios and phase shifts are only reported when the amplitude spectrum confirmed that the imposed frequency is the dominant frequency in the signal. The results for tests that included sequences of zero-flow rate and therefore were manually shifted along the time axis are indicated by italic numbers. Int.  Period  Nominal flow rate  Injectivity analysis  Vertical-interference analysis        Qcor–pint  Qsurface–pint  pbel–pint  pabo–pint        Phase shift  Amplitude ratio  Amplitude ratio  Phase shift  Amplitude ratio  Phase shift  Amplitude ratio    (s)  (L min−1)  (cycles)  (m3 s−1 Pa−1)  (m3 s−1 Pa−1)  (cycles)  (–)  (cycles)  (–)  i1  180  3/0  3.65 ± 0.03E−01  4.54 ± 0.26E−10  2.05 ± 0.02E−09  4.43 ± 0.31E−01  6.90 ± 0.59E−03  6.14 ± 0.35E−01  1.40 ± 0.23E−03        1.08 ± 0.05E−01  5.21 ± 0.12E−10            i2–1  180  10/7  9.89 ± 0.22E−02  4.06 ± 0.09E−09  5.23 ± 0.08E−09  7.71 ± 0.25E−01  2.01 ± 0.30E−01  –  –  i2–2  180  3/0  1.70 ± 0.03E−01  3.76 ± 0.05E−09  5.26 ± 0.01E−09  8.74 ± 0.05E−01  3.44 ± 0.09E−01  –  –  i3–1  180  7/4  2.03 ± 0.05E−01  2.86 ± 0.30E−09  4.23 ± 0.09E−09  8.62 ± 0.03E−01  1.22 ± 0.04E−02  –  –  i3–2  180  3/0  3.24 ± 0.05E−01  1.06 ± 0.03E−09  2.67 ± 0.02E−09  8.91 ± 0.03E−01  1.51 ± 0.03E−02  –  –        2.65 ± 0.05E−01  9.88 ± 0.31E−10            i3–3  480  3/0  1.19 ± 0.02E−01  9.64 ± 0.10E−10  1.47 ± 0.02E−09  8.97 ± 0.04E−01  2.74 ± 0.07E−02  –  –  i3–5  1080  7/4  2.39 ± 0.08E−02  6.19 ± 0.20E−10  7.19 ± 0.23E−10  8.80 ± 0.10E−01  3.85 ± 0.17E−02  –  –  i5–1  60  3/0  3.80 ± 0.08E−01  5.19 ± 0.13E−09  9.43 ± 0.11E−09  9.46 ± 0.06E−01  1.53 ± 0.05E−01  –  –        1.40 ± 0.01E−01  4.90 ± 0.03E−09            i5–2  60  10/7  1.02 ± 0.03E−01  6.60 ± 0.13E−09  1.04 ± 0.01E−08  9.16 ± 0.17E−01  2.24 ± 0.20E−01  –  –  i5–3  180  10/7  7.40 ± 0.17E−02  4.65 ± 0.05E−09  5.63 ± 0.04E−09  8.94 ± 0.28E−01  2.85 ± 0.46E−01  7.04 ± 0.77E−01  7.60 ± 3.40E−03  i5–4  540  10/7  8.40 ± 0.11E−02  3.21 ± 0.03E−09  3.50 ± 0.04E−09  2.62 ± 1.67E−02  2.70 ± 0.15E−01  7.15 ± 0.51E−01  3.50 ± 0.59E−03  Int.  Period  Nominal flow rate  Injectivity analysis  Vertical-interference analysis        Qcor–pint  Qsurface–pint  pbel–pint  pabo–pint        Phase shift  Amplitude ratio  Amplitude ratio  Phase shift  Amplitude ratio  Phase shift  Amplitude ratio    (s)  (L min−1)  (cycles)  (m3 s−1 Pa−1)  (m3 s−1 Pa−1)  (cycles)  (–)  (cycles)  (–)  i1  180  3/0  3.65 ± 0.03E−01  4.54 ± 0.26E−10  2.05 ± 0.02E−09  4.43 ± 0.31E−01  6.90 ± 0.59E−03  6.14 ± 0.35E−01  1.40 ± 0.23E−03        1.08 ± 0.05E−01  5.21 ± 0.12E−10            i2–1  180  10/7  9.89 ± 0.22E−02  4.06 ± 0.09E−09  5.23 ± 0.08E−09  7.71 ± 0.25E−01  2.01 ± 0.30E−01  –  –  i2–2  180  3/0  1.70 ± 0.03E−01  3.76 ± 0.05E−09  5.26 ± 0.01E−09  8.74 ± 0.05E−01  3.44 ± 0.09E−01  –  –  i3–1  180  7/4  2.03 ± 0.05E−01  2.86 ± 0.30E−09  4.23 ± 0.09E−09  8.62 ± 0.03E−01  1.22 ± 0.04E−02  –  –  i3–2  180  3/0  3.24 ± 0.05E−01  1.06 ± 0.03E−09  2.67 ± 0.02E−09  8.91 ± 0.03E−01  1.51 ± 0.03E−02  –  –        2.65 ± 0.05E−01  9.88 ± 0.31E−10            i3–3  480  3/0  1.19 ± 0.02E−01  9.64 ± 0.10E−10  1.47 ± 0.02E−09  8.97 ± 0.04E−01  2.74 ± 0.07E−02  –  –  i3–5  1080  7/4  2.39 ± 0.08E−02  6.19 ± 0.20E−10  7.19 ± 0.23E−10  8.80 ± 0.10E−01  3.85 ± 0.17E−02  –  –  i5–1  60  3/0  3.80 ± 0.08E−01  5.19 ± 0.13E−09  9.43 ± 0.11E−09  9.46 ± 0.06E−01  1.53 ± 0.05E−01  –  –        1.40 ± 0.01E−01  4.90 ± 0.03E−09            i5–2  60  10/7  1.02 ± 0.03E−01  6.60 ± 0.13E−09  1.04 ± 0.01E−08  9.16 ± 0.17E−01  2.24 ± 0.20E−01  –  –  i5–3  180  10/7  7.40 ± 0.17E−02  4.65 ± 0.05E−09  5.63 ± 0.04E−09  8.94 ± 0.28E−01  2.85 ± 0.46E−01  7.04 ± 0.77E−01  7.60 ± 3.40E−03  i5–4  540  10/7  8.40 ± 0.11E−02  3.21 ± 0.03E−09  3.50 ± 0.04E−09  2.62 ± 1.67E−02  2.70 ± 0.15E−01  7.15 ± 0.51E−01  3.50 ± 0.59E−03  View Large 2.1.3 Undisturbed hydraulic heads We constrained the natural hydraulic heads from the pressure levels above, in, and below the interval asymptotically approached with time after packer inflation (Fig. 4). The recorded pressures suggest that the unperturbed water level of about 32 m below the surface, corresponding to about 700 L of water in the well, is actually the consequence of a dynamic situation with a prominent outflow at the borehole bottom, that is, at a true vertical depth of more than 110 m, and an inflow between about 73 and 110 m. Relative to the unperturbed pressure the pumping operations employed overpressures of at most 300 kPa (Fig. 4). Pumping pressures reached the level corresponding to a water column filling the entire borehole only for the interval i5. Figure 4. View largeDownload slide Equilibrium pressure of intervals (pint) and below probe (pbel) compared to pumping pressure in the interval and two hydrostatic water pressures for water levels at the surface (0 m) or at 32 m depth in the well. Figure 4. View largeDownload slide Equilibrium pressure of intervals (pint) and below probe (pbel) compared to pumping pressure in the interval and two hydrostatic water pressures for water levels at the surface (0 m) or at 32 m depth in the well. 2.2 Fourier analysis We performed two basic steps of data processing before Fourier transformation. Recorded time-series were cut to lengths corresponding to integer multiples of the pumping period and pressure records were automatically detrended to avoid artefacts due to differences between pressure levels at start and end of a pumping sequence (Fig. 5). Transformation yield amplitude and phase spectra from which we determined amplitude ratios and phase shifts between corrected flow rate (see Section 2.3.1) and interval pressure but also between interval pressure and pressure below and above the interval for the imposed period (see Table 2). The ratio between the amplitudes of flow rate and interval pressure corresponds to the injectivity index determined from conventional tests (e.g. Lyons 2010). Figure 5. View largeDownload slide Comparison of (a) recorded and (b) detrended pressure during the periodic testing at a depth of about 105 m (i5–3). For detrending, the average of the upper and lower envelopes is subtracted from the signal. Figure 5. View largeDownload slide Comparison of (a) recorded and (b) detrended pressure during the periodic testing at a depth of about 105 m (i5–3). For detrending, the average of the upper and lower envelopes is subtracted from the signal. To constrain the uncertainty of phase shift and amplitude ratio, we performed a sliding-window analysis, that is, a window with a length of three nominal periods was successively moved over the entire data set using a step size of 10–20 per cent of a period. For each window position, phase shift and amplitude ratio was determined and the standard deviation of all values is here reported as uncertainty. 2.3 Injectivity analysis 2.3.1 Correction for the storage capacity of the tubing Injectivity analysis requires the true flow rate into the rock that differs from the flow rate measured at the surface owing to the storage of fluid in the hydraulic system, composed of the tubing and the borehole section enclosed by the packers. For a closed system entirely filled with fluid, its capacity to store fluid stems from the finite fluid compressibility, the deformability of the tubing and of the pressurized rock section in the interval, and changes in packer seating due to changes in interval pressure. For an open system, one has to additionally account for a geometrical storage capacity that specifies how much fluid has to be added to/removed from the water column for a unit change in its height and thus in pressure it exerts. The used tubing with an inner diameter of 23.7 mm exhibits a storage capacity dominated by its geometrical component of $${S_{{\rm{tube}}}} = {{{A_{{\rm{tube}}}}} / {\rho g}} = 4.8 \times {10^{ - 8}}\,{{\rm{m}}^3}\,{\rm{P}}{{\rm{a}}^{ - 1}}$$ where Atube, ρ and g denote cross-section of the tubes (m2), fluid density (kg m−3) and gravitational acceleration (m s−2), respectively. The contribution from the compressibility of water is at least three orders of magnitude smaller, as is the contribution from the deformation of the steel tubes with a thickness of almost 5 mm. We constrained the combined contribution of rock deformation and changes in packer seating to interval-storage capacity by a rapid pulse test for the long interval of 9.4 m length (corresponding to an interval volume of about 60 L). The valve of the probe was briefly opened to pressurize the interval by the water column in the tubing. The change in volume was determined from the change in water-column height in the tubing and the change in pressure was recorded by the interval sensor. The ratio between these changes yields an interval-storage capacity of $${S_{{\mathop{\rm int}} }}( {9.4\ {\rm{m}}} ) = 6.0 \times {10^{ - 10}}\,{{\rm{m}}^3}\,{\rm{P}}{{\rm{a}}^{ - 1}}$$, that is, almost two orders of magnitude smaller than the geometrical storage capacity of the tubing (or 1.3 per cent). The shorter interval length of 1.7 m (corresponding to an interval volume of about 10 L) was not explicitly tested, but the reduced length should lead to a further reduction of the storage capacity roughly proportional to interval length, that is, Sint(1.7 m) ≃ 1.1 × 10 − 10m3Pa − 1. Thus, we neglected all contributions to storage capacity but the geometrical one of the tubing when determining true flow rates into or out of the rock formation from the flow rate recorded at the surface, Qsurface(L min−1), according to   \begin{equation}{Q_{\rm cor}} = {Q_{{\rm{surface}}}} - {S_{{\rm{tube}}}}{\dot{p}_{{\rm{int}}}}.\end{equation} (1) The time derivative of the interval pressure, $${\dot{p}_{{\rm{int}}}}$$, was calculated using a Fourier transformation of the recorded pressure. 2.3.2 Flow delay All tests with a water level in the tube below surface level are potentially affected by the conventional approach of remotely measuring flow rate at the surface and the associated time delay between flow in the gauge and actual addition of fluid to the water column loading the interval. Yet, when correcting flow records after eq. (1), tests employing a ‘zero’-flow rate, which we address as zero-flow rate tests below, yield suspicious delays of up to 10 s between flow rate at the surface and pressure in the interval exceeding (Fig. 6). These large delays presumably result from the actual time of flow in the partly empty tube. These data cannot be evaluated in a standard way. We pursue two different approaches to derive information on hydraulic properties from these tests. First, we determine the amplitude ratio of uncorrected flow rate, that is, surface-flow rate and interval pressure (Table 2) for all tests and compare their relation to pumping characteristics with that of amplitude ratios of the tests for which flow-rate correction was possible. Secondly, we test whether manual shifting of flow records guided by the prominent kinks in pressure records provides reasonable estimates for the phase shift between flow rate and pressure and thus opens the way for an inversion towards hydraulic parameters. Figure 6. View largeDownload slide Interval pressure (blue, left y-axis) and flow rate (orange, right y-axis) recorded during a periodic pumping test with a period of 180 s at a depth of 115 m (i1, see Tables 1 and 2) exhibiting a delay between start (stop) of pumping, indicated by the steep rise (fall) in flow rate, and increase in pressure, indicated by the kinks. Figure 6. View largeDownload slide Interval pressure (blue, left y-axis) and flow rate (orange, right y-axis) recorded during a periodic pumping test with a period of 180 s at a depth of 115 m (i1, see Tables 1 and 2) exhibiting a delay between start (stop) of pumping, indicated by the steep rise (fall) in flow rate, and increase in pressure, indicated by the kinks. 2.3.3 Flow regimes considered for estimation of hydraulic properties Equivalent or effective hydraulic properties, diffusivity D (m2 s−1), transmissivity T (m2 s−1) and storativity S (−), are typically estimated relying on analytical solutions of the pressure-diffusion equation for specific flow regimes, that is, the dimensional and directional characteristics of the flow pattern (e.g. Fetter 2001). Here, we considered: (i) 1-D flow, (ii) radial flow, (iii) radial flow in concentric cylindrical shells and (iv) bilinear flow. Analytical solutions of the 1-D-flow and radial-flow models were previously derived for various boundary conditions (Black & Kipp 1981; Rasmussen et al.2003; Renner & Messar 2006). Below, we investigate the effect of finite interval length on flow in homogeneous isotropic medium and present solutions for the other scenarios. Effect of interval length. The use of a double-packer probe poses a significant difference to previously performed periodic pumping tests mandating to investigate the consequences of a restricted ‘active’ length of the formation. For this purpose, we considered a radially infinite, mirror-symmetric model vertically bounded by a no-flow condition (see Appendix  B). To allow for analytical treatment, we assumed pressure to evolve linearly along the packers and to be constant above the upper (below the lower) packer. According to the model, the effects of interval length on the phase shift between flow rate and interval pressure are stronger the larger the period, that is, the penetration depth (Fig. 7). Relying on the radial-flow model for the inversion of hydraulic diffusivity from phase shift, diffusivity is the more overestimated the shorter the interval and the smaller the phase shift are. Small phase shifts correspond to large periods (small frequencies) for which a double-packer probe acts as a point source rather than a line source. Figure 7. View largeDownload slide Effect of interval length (as given in labels; radial flow corresponds to an infinitely long interval) on phase shift between flow rate and interval pressure according to the mirror-symmetric model described in Appendix  B. Curves were calculated after eq. (B17) with the summation up to n = 100. Figure 7. View largeDownload slide Effect of interval length (as given in labels; radial flow corresponds to an infinitely long interval) on phase shift between flow rate and interval pressure according to the mirror-symmetric model described in Appendix  B. Curves were calculated after eq. (B17) with the summation up to n = 100. Shell model. Analytical solutions for the injectivity analysis considering a shell model, known as multiregion composite model in reservoir engineering (Ambastha & Ramey 1992; Acosta & Ambastha 1994), are given in Appendix  C. These solutions complement the work by Ahn & Horne (2010), who reported a semi-analytical solution for interference analysis of a multicomposite radial ring system. For the analysis of the field data, we focus on a single cylindrical shell concentrically surrounding the borehole and exhibiting hydraulic properties that differ from the medium farther away from the well (comparable to the situation envisioned for the so-called skin effect, Matthews & Russell 1967). The prominent feature of this model is a perturbation of the transition from zero phase shift to a phase shift of 1/4 in comparison to radial flow in a homogeneous medium. The perturbations occur for values of ∼0.01–1 of the dimensionless parameter $${r_{\rm i}}\sqrt {\omega /{D_{{\rm{shell}}}}} $$, where Dshell is the hydraulic diffusivity of the shell, ri the radius of the injection well and ω the angular frequency. This dimensionless parameter gives the inverse of the ratio between the hydraulic penetration depth into the shell material and borehole radius. Prominent maxima, successions of maxima and minima, or a focusing of the transition occur depending on the ratios between the hydraulic parameters representing the shell and the surrounding medium (see e.g.  Appendix C2). Bilinear-flow model. Bilinear flow occurs when fluid is drained from or injected into a permeable matrix through an enclosed fracture of finite conductivity intersecting a well along its axis. A combination of two approximately linear flow regimes may result, one in the matrix with flow essentially perpendicular to the fracture and the other in the fracture itself associated with the non-negligible pressure drop or increase in it (Cinco-Ley et al.1978; Cinco-Ley & Samaniego-V 1981; Ortiz et al. 2013). In Appendix  D, we present an analytical solution of the coupled diffusion equations for bilinear flow in a homogeneous medium subjected to periodic pumping. For infinite fracture length, the phase shift is bounded by asymptotes to 1/16 and 1/8 of a cycle for large and small periods, respectively, with a smooth transition in-between (Fig. D2 in Appendix  D). The asymptotic value for long periods is consistent with the value derived by Hollaender et al. (2002). In case of fractures with finite length, only the asymptote of 1/8 of a cycle for small periods holds, too, as expected because the pressure perturbation does not reach the fracture tip. For a constant-pressure boundary at the fracture tip, phase shift monotonically decreases to 0 with increasing period. The relation between phase shift and period is not monotonic for a no-flow boundary, but minima (and possibly also maxima exceeding 1/8) occur at intermediate periods before 1/8 is asymptotically approached for very large periods when the pressure in the fracture is homogeneous and the total response is dominated by linear flow into the surrounding medium. 2.3.4 Flow-regime diagnosis from spectral analyses For conventional tests, flow regimes can be inferred from diagnostic log–log plots of pressure or pressure derivative versus time (e.g. Renard et al.2009). Analogously, slopes of amplitude spectra of non-harmonic periodic tests allow for a diagnosis of flow regimes (Hollaender et al.2002). In either type of diagram, slopes of 1/4, 1/2, 0 and 1 are indicative of bilinear flow, linear flow, radial flow and pseudo-steady-state flow, respectively. Since our excitation is not truly harmonic, spectra of flow rate and pressure contain significant amplitudes at periods shorter than the nominally excited one, too. Yet, simple division of entire flow-rate and pressure spectra proved to be impractical (Fig. 8). Thus, we restricted to specific local maxima of pressure and flow rate, as did Fokker et al. (2013). Figure 8. View largeDownload slide Spectrum of the amplitude ratio modulus for the interval at 105 m tested with a period of 540 s (i5–4). Coloured lines represent exponents of 1/4, 1/2 and 1. The restriction to local maxima of pressure and flow rate in the frequency domain (red asterisks) reduces the scatter and suggests bilinear (slope 1/4) to linear flow (1/2). Figure 8. View largeDownload slide Spectrum of the amplitude ratio modulus for the interval at 105 m tested with a period of 540 s (i5–4). Coloured lines represent exponents of 1/4, 1/2 and 1. The restriction to local maxima of pressure and flow rate in the frequency domain (red asterisks) reduces the scatter and suggests bilinear (slope 1/4) to linear flow (1/2). 2.4 Vertical-interference analysis When placing a double-packer probe in an open borehole, the pumping operations in the interval may also cause detectable pressure variations above or below the probe due to probe-parallel flow through the formation bearing additional information on its hydraulic properties. To take advantage of these pressure records, we analysed an axisymmetric model (Fig. 9) and treated it numerically using COMSOL. We determined values of phase shift and amplitude ratio of pressures measured at different locations of a double-packer probe considering different diffusivities and combinations of boundary conditions (Table 3). The effect of the investigated scenarios (1–8) is subordinate for the interference response between below lower packer and in the interval (Fig. 10). The results for pressure above upper packer separate into two groups according to the top boundary conditions labeled as ‘no flow’ and ‘constant pressure’. Figure 9. View largeDownload slide Illustration of the model used for numerical analysis of vertical interference. Blue colour indicates water-filled well sections and the two packers are represented in grey. The sketch is not to scale: distance of the lateral boundaries from the symmetry axis is 1000ri; distance between top and bottom boundary is 240 m, that is, the bottom boundary is more than 100 m away from the bottom of the borehole. Figure 9. View largeDownload slide Illustration of the model used for numerical analysis of vertical interference. Blue colour indicates water-filled well sections and the two packers are represented in grey. The sketch is not to scale: distance of the lateral boundaries from the symmetry axis is 1000ri; distance between top and bottom boundary is 240 m, that is, the bottom boundary is more than 100 m away from the bottom of the borehole. Figure 10. View largeDownload slide Phase shift and amplitude ratio of the interference response between pressure recorded (a) below the lower packer and (b) above the upper packer and in the interval calculated using the model described in Fig. 9 with different diffusivities and boundary conditions as specified in Table 3 (expressed by the numbers used as labels). The dots represent individual calculations labeled by the used diffusivity (m2 s−1). The lines intended to help the eye to follow the trends. The red square with error bars represents the observations at a depth of 85 m for a period of 180 s (i5–3). Figure 10. View largeDownload slide Phase shift and amplitude ratio of the interference response between pressure recorded (a) below the lower packer and (b) above the upper packer and in the interval calculated using the model described in Fig. 9 with different diffusivities and boundary conditions as specified in Table 3 (expressed by the numbers used as labels). The dots represent individual calculations labeled by the used diffusivity (m2 s−1). The lines intended to help the eye to follow the trends. The red square with error bars represents the observations at a depth of 85 m for a period of 180 s (i5–3). Table 3. Combinations of boundary conditions considered for the performed eight calculations employing the model sketched in Fig. 9. Boundary conditions  1  2  3  4  5  6  7  8  Top boundary  No flow  No flow  No flow  No flow  Constant pressure  Constant pressure  Constant pressure  Constant pressure  Bottom boundary  No flow  Constant pressure  Constant pressure  No flow  No flow  Constant pressure  Constant pressure  No flow  Lateral boundary  Constant pressure  Constant pressure  No flow  No flow  Constant pressure  Constant pressure  No flow  No flow  Boundary conditions  1  2  3  4  5  6  7  8  Top boundary  No flow  No flow  No flow  No flow  Constant pressure  Constant pressure  Constant pressure  Constant pressure  Bottom boundary  No flow  Constant pressure  Constant pressure  No flow  No flow  Constant pressure  Constant pressure  No flow  Lateral boundary  Constant pressure  Constant pressure  No flow  No flow  Constant pressure  Constant pressure  No flow  No flow  View Large Effective diffusivity values Dδ or Dφ are estimated from amplitude ratio or phase shift, respectively, by linearly interpolating between the modeling results. Deviation of model results from observations could be explained by several aspects. The model assumes homogeneity and isotropy for the subsurface which may not be true in reality. Furthermore, the model only considers relations between pressures but a more comprehensive model should also integrate the injectivity analysis. Finally, hydromechanical effects were observed but not considered in the current model. 3 RESULTS 3.1 Injectivity analysis Amplitude ratios of data sets resulting from performing the flow-rate correction after eq. (1) are systematically smaller than those of uncorrected data (Fig. 11), but the relative relations among the ratios for the different test intervals as well as their trends with the parameters characterizing a test, mean pressure and period, remain qualitatively similar. We therefore use the larger set of uncorrected data for the subsequent analysis and refer to this approach as simplified injectivity analysis. Figure 11. View largeDownload slide Injectivity quantified by the amplitude ratio of flow rate and injection pressure as dependent on mean interval pressure. Amplitude ratios were calculated for all tests before (labeled uncorrected and represented by solid circles) and after (labeled corrected and represented by open circles) applying the correction to flow rate for the storage of the tubing. The colour coding represents the five investigated test intervals (see Table 1). Data points with uncorrected flow rate are also labeled by the imposed oscillation period in seconds. Vertical error bars indicate the uncertainty in amplitude ratio gained from the sliding-window fast Fourier transform (FFT) analyses. Horizontal error bars represent the range of pressures at which the tests were conducted. Figure 11. View largeDownload slide Injectivity quantified by the amplitude ratio of flow rate and injection pressure as dependent on mean interval pressure. Amplitude ratios were calculated for all tests before (labeled uncorrected and represented by solid circles) and after (labeled corrected and represented by open circles) applying the correction to flow rate for the storage of the tubing. The colour coding represents the five investigated test intervals (see Table 1). Data points with uncorrected flow rate are also labeled by the imposed oscillation period in seconds. Vertical error bars indicate the uncertainty in amplitude ratio gained from the sliding-window fast Fourier transform (FFT) analyses. Horizontal error bars represent the range of pressures at which the tests were conducted. 3.1.1 Simplified injectivity analysis Amplitude ratios of uncorrected flow rates and interval pressures, here treated as tentative measures of injectivity, tend to increase with increasing mean pressure but variations at a given pressure are significant (Fig. 11). A period of 180 s was chosen for more than half of the periodic pumping tests allowing us to compare the hydraulic behaviour of the four tested intervals. The deepest interval with a depth of 115 m (i1) exhibits the lowest injectivity despite its association with the largest absolute mean pressure. The shallowest interval at 85 m depth (i3) shows a strong positive correlation between injectivity and mean pressure (also found for a period of 60 s applied in i5) indicating that the dominant hydraulic conduits are perceptible to hydromechanical effects. For a depth of 105 m (i2, i5), we neither find a pronounced dependence of amplitude ratio on mean pressure nor on interval length. Amplitude ratios decrease systematically with increasing pumping period for i5 as well as for i3. 3.1.2 Effective hydraulic parameters Effective transmissivity is the least variable of all hydraulic parameters when evaluating the results of the Fourier analyses of corrected flow rate and interval pressure relying on the analytical solution for radial flow in a homogeneous medium (Table 4). It varies by about half an order of magnitude from 3 to 9 × 10−6 m2 s−1. Apart from one outlier (i3–5), the values for hydraulic diffusivity and storativity stay within about one order of magnitude (4–40 × 10−5 m2 s−1 and 2–20 × 10−3, respectively). The outlier may to some extent be explained by the finite length of the interval. The modeling of the length effect (see Appendix  B) predicts the largest error for test i3–5, with diffusivity potentially overestimated by about one order of magnitude (0.37 m2 s−1 from conventional radial-flow model versus 0.04 m2 s−1 from mirror-symmetric double-packer model). For the other tests, yielding phase shifts of at least 0.074 cycles, the error due to neglecting the finite interval length is negligible (Fig. 7). Table 4. Hydraulic properties derived from the results of the Fourier analyses of flow rate and interval pressure relying on the analytical solution for radial flow in a homogeneous medium. Int.  Period  Nominal flow rate  Depth  Diffusivity  Transmissivity  Storativity    (s)  (L min−1)  (m)  (m2 s−1)  (m2 s−1)  (–)  i1  180  3/0  117.1  1.3 ± 0.9E−05  2.7 ± 0.8E−07  2.7 ± 1.3E−02  i2–1  180  10/7  106.7  3.9 ± 0.9E−05  3.5 ± 0.3E−06  9.4 ± 1.5E−02  i3–3  480  3/0  81.9  5.0 ± 3.1E−07  1.8 ± 0.5E−07  4.4 ± 1.8E−01  i3–5  1080  7/4  81.9  3.7 ± 1.3E−01  5.1 ± 0.2E−06  1.5 ± 0.5E−05  i5–2  60  10/7  105.0  8.4 ± 3.4E−05  4.9 ± 0.8E−06  6.4 ± 1.5E−02  i5–3  180  10/7  105.0  3.9 ± 0.3E−04  9.0 ± 0.2E−06  2.3 ± 0.2E−02  i5–4  540  10/7  105.0  5.7 ± 0.6E−05  4.8 ± 0.1E−06  8.5 ± 0.6E−02  Int.  Period  Nominal flow rate  Depth  Diffusivity  Transmissivity  Storativity    (s)  (L min−1)  (m)  (m2 s−1)  (m2 s−1)  (–)  i1  180  3/0  117.1  1.3 ± 0.9E−05  2.7 ± 0.8E−07  2.7 ± 1.3E−02  i2–1  180  10/7  106.7  3.9 ± 0.9E−05  3.5 ± 0.3E−06  9.4 ± 1.5E−02  i3–3  480  3/0  81.9  5.0 ± 3.1E−07  1.8 ± 0.5E−07  4.4 ± 1.8E−01  i3–5  1080  7/4  81.9  3.7 ± 1.3E−01  5.1 ± 0.2E−06  1.5 ± 0.5E−05  i5–2  60  10/7  105.0  8.4 ± 3.4E−05  4.9 ± 0.8E−06  6.4 ± 1.5E−02  i5–3  180  10/7  105.0  3.9 ± 0.3E−04  9.0 ± 0.2E−06  2.3 ± 0.2E−02  i5–4  540  10/7  105.0  5.7 ± 0.6E−05  4.8 ± 0.1E−06  8.5 ± 0.6E−02  View Large Phase-shift values are particularly diagnostic for hydraulic boundaries. For example, phase shifts larger than 1/8 unequivocally require bounded reservoirs (Fig. 12). The phase-shift value for the test at a depth of 85 m and a period of 180 s (i3–2) remains larger than the upper bound of 1/4 of all considered bounded models even after manual shifting. Only models with a no-flow boundary can explain the phase shift of the test i3–1 at 85 m and a period of 180 s, the largest phase shift observed for non-zero-flow rate tests (Fig. 12). However, the increase of amplitude ratios with phase shift between i3–1 and the other test at the same depth for a period of 1080 s (i3–5) cannot be explained by any of the no-flow models. Independent of interval length the phase-shift values for a depth of ∼105 m (i2, i5) are all smaller than 1/8 of a cycle, excluding the 1-D-flow model with a no-flow boundary or no boundary. The non-monotonically varying amplitude ratios for i5 also exclude a constant pressure boundary or no boundary for all the flow types. For the interval with a depth of ∼115 m (i1), the only phase shift lies between 1/16 and 1/8, which can be fit by all the models except 1-D flow with no-flow boundary or no boundary. Figure 12. View largeDownload slide Amplitude ratio versus phase shift for the injectivity analysis. Red, blue and black lines represent the solutions for radial flow, 1-D flow and bilinear flow, respectively. Solid lines are for infinite reservoirs and dashed and dashed–dotted lines represent constant pressure (pb) and no-flow (fb) boundary at distances of 20ri, respectively, for radial and bilinear flows. Changing the assumed hydraulic parameters shifts the theoretical curves vertical without distortion. Filled and open symbols represent non-zero-flow-rate sequences and zero-flow-rate sequences for which manual shift was applied, respectively (see Table 2). Errors for the data points do not exceed symbol size. Two periods are also labeled for the convenience of discussion (see Table 4). Figure 12. View largeDownload slide Amplitude ratio versus phase shift for the injectivity analysis. Red, blue and black lines represent the solutions for radial flow, 1-D flow and bilinear flow, respectively. Solid lines are for infinite reservoirs and dashed and dashed–dotted lines represent constant pressure (pb) and no-flow (fb) boundary at distances of 20ri, respectively, for radial and bilinear flows. Changing the assumed hydraulic parameters shifts the theoretical curves vertical without distortion. Filled and open symbols represent non-zero-flow-rate sequences and zero-flow-rate sequences for which manual shift was applied, respectively (see Table 2). Errors for the data points do not exceed symbol size. Two periods are also labeled for the convenience of discussion (see Table 4). At a depth of ∼105 m, three tests were performed with different periods using the short interval (60, 180 and 540 s corresponding to i5–2, i5–3, and i5–4 in Table 2). Since they all included only non-zero flow rates their results are not affected by the problems associated with the finite distance between the location of the flow meter and the top of the water column loading the interval. The non-monotonic succession of phase-shift values, with the lowest occurring for the intermediate period, can be modeled by a shell with a thickness of about 6–15 times the borehole radius of ∼0.05 m, and with a diffusivity between just below 10−3 to above 10−2 m2 s−1 (Table 5) exceeding that for the enclosing medium by a factor between about 5 and 10 (Fig. 13). Figure 13. View largeDownload slide Examples phase-shift–frequency relations for the shell model calculated using eq. (C2) in relation to observations for the periodic tests at 105 m for a period of 180 s (i5–3) and for a period of 540 s (i5–4). The labels give the ratios between the diffusivity of the shell and the surrounding medium and the lateral extension of the shell in multiples of the borehole radius. Figure 13. View largeDownload slide Examples phase-shift–frequency relations for the shell model calculated using eq. (C2) in relation to observations for the periodic tests at 105 m for a period of 180 s (i5–3) and for a period of 540 s (i5–4). The labels give the ratios between the diffusivity of the shell and the surrounding medium and the lateral extension of the shell in multiples of the borehole radius. Table 5. Parameters of the shell model (see Appendix  C) derived from the phase-shift values observed for the tests i5–3 and i5–4 (see also Fig. 13). k1/k2  D1/D2  r1/ri  Diffusivity (D1)  (–)  (–)  (–)  (m2 s−1)  6.7  8  15  8.1E−04  7.1  10  15  7.4E−04  6.3  10  10  7.1E−03  5.9  9  6.1  2.3E−02  5.6  9  6.1  2.1E−02  k1/k2  D1/D2  r1/ri  Diffusivity (D1)  (–)  (–)  (–)  (m2 s−1)  6.7  8  15  8.1E−04  7.1  10  15  7.4E−04  6.3  10  10  7.1E−03  5.9  9  6.1  2.3E−02  5.6  9  6.1  2.1E−02  View Large 3.2 Flow regime diagnosis from spectral analyses Flow-regime identification from full spectra proved difficult (Fig. 8), not the least because the spectra of flow rate and pressure contain little significant contributions for high frequencies. We thus restricted to the first two odd multiples of the test frequency that correspond to pronounced local maxima in the spectra. As the analysis is sensitive to the timing of flow-rate changes (Hollaender et al.2002), we only analyse non-zero-flow-rate sequences. The results for overtones from the various tests in a specific interval agree closely with those for nominally excited periods (Fig. 14). The slope for the interval at 105 m (i5) increases from less than 1/4 and approaches 1/2 as period increases. In contrast, the slope for interval i3 decreases from around 1 to a value just above 1/4. At face value, these results indicate a transition from bilinear to linear flow for i5 and a transition from pseudo-steady-state flow to bilinear flow for i3. The physical interpretation for a transition from bilinear to linear flow is that the pressure at the fracture tip has risen to a substantial fraction of the associated pressure increase in the well at the end of bilinear flow and the flow field transforms towards formation linear flow. The transition from pseudo-steady-state flow to bilinear flow is consistent with a succession of the domination of wellbore storage at early time by bilinear flow (finite conductivity fracture). Figure 14. View largeDownload slide Amplitude ratio in frequency domain for the deepest interval (i5) represented by asterisks and the interval at depth 81.9 m (i3) represented by circles at the main period of the pumping tests and the first overtone, that is, a third of the imposed period. The dotted lines represent exponents of 1/4, 1/2 and 1. Figure 14. View largeDownload slide Amplitude ratio in frequency domain for the deepest interval (i5) represented by asterisks and the interval at depth 81.9 m (i3) represented by circles at the main period of the pumping tests and the first overtone, that is, a third of the imposed period. The dotted lines represent exponents of 1/4, 1/2 and 1. 3.3 Vertical-interference analysis The periodic excitations in the injection intervals lead to detectable periodic responses in the pressures recorded below the lower and above the upper packer (Table 2). The vertical-interference analysis based on these pressure records has the advantage that all tests can be analysed even those involving the problematic zero-flow-rate sequences. Yet, the magnitude of interference response in a single borehole critically depends on the storage capacity of the borehole section in which it is determined. Above the upper packer we had a free, unconstrained water column corresponding to a ‘large’ storage capacity determined by the cross-section of the borehole. The borehole section below the lower packer, in contrast, constitutes an enclosed fluid volume with a ‘small’ storage capacity determined by the fluid compressibility and the deformability of packer and borehole wall. In agreement with these qualitative storage considerations, the interference response is much larger below the probe than above the probe where it is in fact detectable only in some cases (Table 2). All amplitude ratios derived from pressures below the lower packer, but that for the deepest test interval (i1), follow the expected inverse correlation with the length of the enclosed borehole section (Fig. 15): the larger the enclosed volume the larger the storage capacity and the less sensitive the section, that is, the amplitude ratio between interference response and excitation decreases with the length of the enclosed section. Relative to the general trend, the response below the probe observed for the test in the deepest interval is weak demonstrating exceptional hydraulic properties near the borehole bottom. Figure 15. View largeDownload slide Correlation between interference responses, quantified by the amplitude ratio of pressure oscillations, and the length of the enclosed section below the lower packer of the probe. We present the amplitude ratio between the pressure above the interval and in the interval in this graph, too, though it should not be controlled by the length of the enclosed section below the probe. The joint presentation of results is chosen to highlight the difference in sensitivity of vertical response between the enclosed section below the probe and the unconstrained water column above it. Figure 15. View largeDownload slide Correlation between interference responses, quantified by the amplitude ratio of pressure oscillations, and the length of the enclosed section below the lower packer of the probe. We present the amplitude ratio between the pressure above the interval and in the interval in this graph, too, though it should not be controlled by the length of the enclosed section below the probe. The joint presentation of results is chosen to highlight the difference in sensitivity of vertical response between the enclosed section below the probe and the unconstrained water column above it. Including phase shift in the comparison, observations for the interference with the borehole section above the interval are only in agreement with model results for a ‘no-flow’ condition (Fig. 10). The effective diffusivities estimated from phase shift (denoted Dφ in Table 6) tend to be larger than the ones estimated from amplitude ratio (Dδ) for i2 and i5, the long and short intervals at 105 m. The reverse seems to hold for the shallowest interval i3. Effective diffusivity values using pbel–pint data or consistently on the order of a few m2 s−1 for the long and the short interval (i2 and i5) at the same depth, except for two relatively small values (i2–1, i5–4). Those estimated by pabo–pint analysis are smaller than those estimated by pbel–pint for the short interval i5 without considering the extremely small Dφ value for i5–4. A decreasing trend of diffusivity with increasing period is indicated for i5. Table 6. Effective diffusivity Dδ estimated by amplitude ratio and Dφ by phase shift of vertical-interference analysis. Top and bottom boundaries are set as no flow. Int.  Period   Nominal flow rate  Dδ (m2 s−1)  Dφ (m2 s−1)    (s)  (L min−1)  pbel–pint  pabo–pint  pbel–pint  pabo–pint  i1  180  3/0  5E−3  5E−3  3E−3  4E−3  i2–1  180  10/7  1  –  3E−2  –  i2–2  180  3/0  4  –  2  –  i3–1  180  7/4  2E−2  –  3E−1  –  i3–2  180  3/0  3E−2  –  4E−1  –  i3–3  480  3/0  3E−2  –  3E−1  –  i3–5  1080  7/4  2E−2  –  1E−1  –  i5–1  60  3/0  4  –  8  –  i5–2  60  10/7  7  –  6  –  i5–3  180  10/7  5  9E−2  3  3E−2  i5–4  540  10/7  3  1E−2  3E−4  7E−3  Int.  Period   Nominal flow rate  Dδ (m2 s−1)  Dφ (m2 s−1)    (s)  (L min−1)  pbel–pint  pabo–pint  pbel–pint  pabo–pint  i1  180  3/0  5E−3  5E−3  3E−3  4E−3  i2–1  180  10/7  1  –  3E−2  –  i2–2  180  3/0  4  –  2  –  i3–1  180  7/4  2E−2  –  3E−1  –  i3–2  180  3/0  3E−2  –  4E−1  –  i3–3  480  3/0  3E−2  –  3E−1  –  i3–5  1080  7/4  2E−2  –  1E−1  –  i5–1  60  3/0  4  –  8  –  i5–2  60  10/7  7  –  6  –  i5–3  180  10/7  5  9E−2  3  3E−2  i5–4  540  10/7  3  1E−2  3E−4  7E−3  View Large 3.4 Conventional methods Conventional tests yield the smallest transmissivity and storativity for the shallowest interval (i4, Table 7). Disregarding the exceptionally low value found for interval i4, effective transmissivity is the least variable hydraulic property, consistent with the observations from periodic radial-flow analysis. The hydraulic parameters derived from the conventional methods and periodic radial-flow analysis are of similar order of magnitude (Table 4). Flow regimes were constrained using diagnostic log–log plots of pressure and pressure derivative versus elapsed time. Intervals i1 und i3 indicate bilinear flow and interval i2 and i5 display features of linear flow. Table 7. Hydraulic properties derived from conventional tests. Int.  Test phase  Diffusivity (m2 s−1)  Transmissivity (m2 s−1)  Storativity (–)  Remarksa  i1  Shut-in   <1.14E−04   <1.6E−06  1.4E−02  RSLA  i2  Constant flow rate   <1.39E−04   <4.3E−06  3.1E−02  SLA  i3  Shut-in   <2.29E−04   <1.9E−06  8.3E−03  RSLA  i4  Shut-in   <1.45E−03   <2.9E−08  2.0E−05  RSLA  i5  Shut-in   <3.18E−05   <2.7E−06  8.5E−02  RSLA  Int.  Test phase  Diffusivity (m2 s−1)  Transmissivity (m2 s−1)  Storativity (–)  Remarksa  i1  Shut-in   <1.14E−04   <1.6E−06  1.4E−02  RSLA  i2  Constant flow rate   <1.39E−04   <4.3E−06  3.1E−02  SLA  i3  Shut-in   <2.29E−04   <1.9E−06  8.3E−03  RSLA  i4  Shut-in   <1.45E−03   <2.9E−08  2.0E−05  RSLA  i5  Shut-in   <3.18E−05   <2.7E−06  8.5E−02  RSLA  aSLA: ‘Straight Line Analysis’ after Jacob & Lohman (1952). RSLA: ‘Recovery Straight Line Analysis’ after Agarwal (1980). View Large 4 DISCUSSION The variable differences between equilibrium pressure for a certain depth interval and a hydrostatic gradient indicate heterogeneity of the hydraulic system. The fractured aquifer is dominated by steep conduits that appear to have limited interconnectivity despite their small lateral distance (at most a few metres in the case of the borehole's bottom section). The diffusivity variations with pumping periods suggest spatial heterogeneity or complexity of the conduits in a specific interval, too, that is, simple radial flow is a poor approximation of the flow excited during the tests. The prominent role of discrete hydraulic conduits is documented by the similarity of injectivity values for the two tests at a borehole depth of about 105 m (i2 and i5) despite their differences in interval length. In the following, we discuss to what extent the employed methods yield comparable values for hydraulic parameters, comment on the constraints for the pressure dependence of the hydraulic response, and finally speculate on the flow-regime model that is most consistent with the entirety of observations in a specific interval. 4.1 Comparison of methods Periodic pumping tests can be easily implemented using field equipment for conventional testing. The excitation signals do not have to be perfectly harmonic but Fourier analysis of non-harmonic signals might actually provide information regarding periods shorter than the imposed main period (spectral analyses, e.g. Hollaender et al.2002; Fokker et al.2013). The periodic excitation can be applied even when the hydraulic pressure is not equilibrated, in fact it can be superposed to any transient, be it associated with a terminated or ongoing pumping operation. The superposition is valid as long as linearity can be assumed, that is, the diffusion equation in its basic form holds. If it does not hold, the results of conventional methods become questionable, too. As a consequence of the freedom to start periodic testing at any time, operational time is reduced and to a larger degree plannable than for conventional tests that involve long waiting times for equilibration of a priori unknown duration. During the current test campaign, a total of 11 hr was spent on conventional tests and 7 hr for periodic tests. Our results show that the analyses of periodic pumping provide effective hydraulic parameters comparable to those gained from conventional methods when assuming radial flow (Fig. 16). Only for the shallowest of the tested intervals (i3), diffusivity values derived from periodic testing deviate significantly from the one gained from conventional testing. The lower diffusivity value results from a test sequence with zero flow rate (i3–3) and thus has a limited reliability. The higher diffusivity value corresponds to test i3–5 giving a phase shift so low that it falls within the range for which the finite length of the probe cannot be neglected, but the flow field has a significant axial component (Fig. 7). The estimate based on the mirror-symmetric double-packer model is more than one order of magnitude smaller than the value assuming simple radial flow. Furthermore, the analyses of periodic tests demonstrates a dependence of the effective hydraulic parameters on period and mean pressure. Thus, differences between the methods likely reflect the simplifications inherent in the radial-flow model that are not justified for the investigated structures necessitating to investigate more complex flow models. Figure 16. View largeDownload slide (a) Radial diffusivity determined from periodic injectivity analysis and conventional methods. (b) Axial diffusivity determined by vertical-interference analysis. Labels refer to the effective diffusivity values estimated from amplitude ratio (Ddel: Dδ) and phase shift (Dphi: Dφ) using pressure records above (pabo), in (pint), and below (pbel) the interval. Figure 16. View largeDownload slide (a) Radial diffusivity determined from periodic injectivity analysis and conventional methods. (b) Axial diffusivity determined by vertical-interference analysis. Labels refer to the effective diffusivity values estimated from amplitude ratio (Ddel: Dδ) and phase shift (Dphi: Dφ) using pressure records above (pabo), in (pint), and below (pbel) the interval. Vertical-interference analysis is a valuable by-product of periodic testing when packer-separated intervals are investigated and pressure is monitored at several depth levels. We refer to its results as effective axial diffusivity to distinguish them from results of the injectivity analysis and conventional methods, addressed as effective radial diffusivity. The attributes ‘axial’ and ‘radial’ are chosen from the perspective of the pumping well, addressing flow dominantly subparallel and normal to well axis, respectively. Effective axial diffusivities are here found to be generally larger than effective radial diffusivities (Fig. 16). The axial diffusivities scatter significantly for the intervals i2 and i5 at similar mean depths but with different probe lengths. The scatter is larger for the shorter interval i5 tested for a larger spread in periods and apart from the result for the longest period of 540 s, communication to the deeper borehole section is characterized by larger diffusivity values than to the shallower sections. The vertical response to shallower depth sections is actually lost with the increase in probe length suggesting that faults outside the short interval but inside the long interval are hydraulically connected to the prominent fault isolated by the short interval. The decrease in diffusivity with increasing period or penetration depth (Table 6) also indicates that the axial connectivity is a local phenomenon. 4.2 Pressure dependence of hydraulic properties The tested intervals exhibit a general trend of increasing injectivity with increasing pressure, with the prominent exception of the deepest interval (i1) that also has the lowest injectivity. In the current geological situation, an increase in injectivity with mean pressure likely results from the pressure-induced increase in effective hydraulic aperture of the joints and faults intersecting the borehole. The more pronounced pressure dependence of the short interval i5 compared to the long interval i2 supports the notion that the pressure dependence reflects hydromechanical behaviour of discrete faults. Observations for i3 (∼85 m depth) with a relatively large period of 180 s suggest that the hydromechanical effects are not restricted to the very intersection between faults and the well. The significance of hydromechanical effects is also evidenced by the pronounced reverse pressure response observed below the lower packer when injecting in i1 (not documented here). 4.3 Constraints on flow regime Details of our observations are obviously at conflict with simple radial flow. From periodic tests, flow regimes can be constrained by (i) investigating spectra of non-harmonic signals (Fig. 14) or by (ii) comparing the variation patterns of amplitude ratio and phase shift with theoretical flow models (Fig. 12). The few constraints we derived from the analysis of spectra yield flow regimes identical to those indicated by the conventional methods. The periodic injectivity analysis tends to indicate more complex flow regimes than the conventional methods (Table 8) which we interpret as evidence for the superior resolution of the periodic approach. Table 8. Constraints on flow regimes gained from injectivity analysis, spectral analysis, or conventional methods: possible (o), excluded (x), no constraint (–) for radial (R), bilinear (BL) and 1-D (1-D) flow.     Periodic pumping method  Conventional methods      Injectivity analysis  Spectral analysis (see Fig. 14)        No flow  Constant pressure  No boundary  No boundary  No boundary  i1  R  o  o  o  –  x    1-D  x  o  x  –  x    BL  o  o  o  –  o  i2  R  o  o  o  –  x    1-D  x  o  x  –  o    BL  x  x  x  –  x  i3  R  x  x  x  x  x    1-D  x  x  x  x  x    BL  x  x  x  o  o  i5  R  o  x  x  x  x    1-D  x  x  x  o  o    BL  o  x  x  o  x      Periodic pumping method  Conventional methods      Injectivity analysis  Spectral analysis (see Fig. 14)        No flow  Constant pressure  No boundary  No boundary  No boundary  i1  R  o  o  o  –  x    1-D  x  o  x  –  x    BL  o  o  o  –  o  i2  R  o  o  o  –  x    1-D  x  o  x  –  o    BL  x  x  x  –  x  i3  R  x  x  x  x  x    1-D  x  x  x  x  x    BL  x  x  x  o  o  i5  R  o  x  x  x  x    1-D  x  x  x  o  o    BL  o  x  x  o  x  View Large Effective injectivity values determined relying on simple radial flow (Fig. 11) systematically decrease with increasing pumping period for intervals i3 and i5 (factor of 3–5 for less than one order of magnitude in period). The conventional scaling relation for diffusion processes (e.g. Weir 1999) suggests that the covered range in period corresponds to a change of half an order of magnitude in nominal penetration depth. Thus, the period dependence of injectivity indicates a change in hydraulic characteristics with distance from the borehole. When addressing the period dependence by the shell model for the deepest interval (i5), we find a decrease in effective diffusivity and permeability by about one order of magnitude beyond a zone of about several decimetres to one metre (Fig. 13). This change in hydraulic characteristics could be related to drilling-associated damage of the borehole wall or remnants of the thixotropic polymer fluid used during drilling. The flushing of the well during logging preceding the hydraulic testing may have well resulted in a cleansing of only the hydraulic conduits nearest to the borehole. For the short interval i5, the conventional method and the spectral analysis (Fig. 14) suggest a linear flow regime while injectivity analysis indicates radial or bilinear flow with a no-flow boundary. For the longer interval at the same depth (i2), the injectivity analysis suggests radial flow with any boundary condition or linear flow with constant pressure boundary. Assuming the two intervals to be governed by the same boundary condition, the only regime in common by i5 and i2 is radial flow with a no-flow boundary. For the deepest interval (i1), the injectivity is the lowest of all intervals though it was tested at the highest mean injection pressure. The low equilibrium pressure below i1 (Fig. 4) suggests that the section below i1 is disconnected from the ones above. In addition, higher axial diffusivity than radial diffusivity (Fig. 16) implies axially oriented or subvertical upflow in the sections above i1. 4.4 Synopsis The periodic testing revealed that the flow geometry in the penetrated subsurface is quite variable on the metre (lateral) to decametre (axial) scale and that the conduit system exhibits hydromechanical effects, two characteristics that appear quite reasonable for a strike-slip fault in crystalline rock. Yet, at this point, it is difficult to distinguish between a situation characterized by a lateral no-flow boundary, possibly related to relics of the polymer drilling fluid with high viscosity or drilling-associated near-welbore damage, and a truly anisotropic hydraulic system with the high axial diffusivity due to subvertically oriented fluid pathways. The latter notion is not only supported by the logging observations (Appendix  A) that document a dominance of steeply dipping fractures intersecting the well, but is also consistent with the conduit geometry deduced by Belgrano et al. (2016), who suggested that the fault zone is controlled by localized subvertically oriented ‘pipe’-like upflow zones. The vertical variability observed for the undisturbed hydraulic heads and the low hydraulic transmissivity of the deepest test interval compared to the major fault zone support this suggestion derived solely from structural investigations. The performed pumping tests correspond to nominal lateral penetration depths of up to a few tens of metres, that is, an order of magnitude smaller than the length scale of a subsurface heat exchanger required for an economically relevant petrothermal energy system. The hydraulic properties near the borehole, however, control the very access to a heat exchanger and thus their characterization, as performed here, bears important consequences for operation. The complexity of flow between formation and well indicated by our results bears significant implications for the average near-well pressure gradient and thus for effective inlet-pressure drops (e.g. Haraden 1992) that in turn bear consequences for operational parameters and bulk flow rates. Also, the realization of the relevance of hydromechanical effects should find observance in modeling and simulation, since their operation is not limited to the vicinity of boreholes (Vinci et al.2015). 5 CONCLUSIONS Periodic pumping tests were conducted using a double-packer probe placed at four different depth levels in a borehole penetrating a hydrothermally active fault. The results of the periodic testing, expressed as amplitude ratio and phase shift between recorded pressures and/or flow rate, permit to constrain effective hydraulic parameters and the flow geometry. We performed and extended the previously established injectivity analysis, but also introduced a new approach, addressed as vertical-interference analysis. We demonstrated that the analyses of the periodic pumping tests yield results comparable to conventional methods when relying on the same model. The field campaign revealed, however, several methodological advantages of the periodic testing, for example, its plannable testing time and the ease in separating the response to the actual pumping operation from noise and transients gained from data processing in frequency space. Our analytical and numerical modeling of periodic testing extends the suite of evaluation tools and specifically allowed us to decipher more details of the encountered complex flow geometry and further conduit characteristics than the conventional methods. The two approaches, the previously employed injectivity analysis and the vertical-interference analysis, are complementary by providing information on the hydraulic properties in different directions relative to the well axis. The variations in effective ‘axial’ and ‘radial’ diffusivity with pumping periods and mean pumping pressure indicate significant spatial heterogeneity of conduits that also exhibit hydromechanical effects for the investigated site. While we showed the potential of the vertical-interference analysis, clearly more work is needed in the future regarding the modeling of records available when using double- or multipacker probes. ACKNOWLEDGEMENTS The financial support by BfE, Switzerland is gratefully acknowledged as is the coordinator, Marco Herwegh, for involving us in project NFP70-P1. Field tests were performed together with Sacha Reinhardt and Markus Bosshard (SolExperts), and Daniel Eggli (Uni Bern) who is particularly thanked for generously sharing the results of the structural core and log analyses. The presented treatment of the shell model builds on an unpublished analyses of Eugen Petkau. We are grateful to the three anonymous reviewers for their valuable comments. REFERENCES Acosta L.G., Ambastha A.K., 1994. Thermal well test analysis using an analytical multi-region composite reservoir model, in SPE Annual Technical Conference and Exhibition , 25–28 September, New Orleans, Louisiana. doi:10.2118/28422-MS. Agarwal R.G., 1980. A new method to account for producing time effects when drawdown type curves are used to analyze pressure buildup and other test data, in SPE Annual Technical Conference and Exhibition , 21–24 September, Dallas, Texas. doi:10.2118/9289-MS. Ahn S., Horne R.N., 2010. Estimating permeability distributions from pressure pulse testing, in SPE Annual Technical Conference and Exhibition , 19–22 September, Florence, Italy. doi:10.2118/134391-MS. Ambastha A.K., Ramey H.J.Jr., 1992. Pressure transient analysis for a three-region composite reservoir, in SPE Rocky Mountain Regional Meeting , 18–21 May, Casper, Wyoming. doi:10.2118/24378-MS. Bakhos T., Cardiff M., Barrash W., Kitanidis P.K., 2014. Data processing for oscillatory pumping tests, J. Hydrol. , 511, 310– 319. Google Scholar CrossRef Search ADS   Baria R., Baumgärtner J., Rummel F., Pine R.J., Sato Y., 1999. HDR/HWR reservoirs: concepts, understanding and creation, Geothermics , 28, 533– 552. Google Scholar CrossRef Search ADS   Becker M.W., Guiltinan E., 2010. Cross-hole periodic hydraulic testing of inter-well connectivity, in Proceedings Thirty-Fifth Workshop on Geothermal Reservoir Engineering, California , Vol. 35, pp. 292– 297, Stanford University, Stanford Geothermal Program. Belgrano T.M., Herwegh M., Berger A., 2016. Inherited structural controls on fault geometry, architecture and hydrothermal activity: an example from Grimsel Pass, Switzerland, Swiss J. Geosci. , 109, 345– 364. Google Scholar CrossRef Search ADS   Bizzarri A., 2012. Effects of permeability and porosity evolution on simulated earthquakes, J. Struct. Geol. , 38, 243– 253. Google Scholar CrossRef Search ADS   Black J.H., Kipp K.L., 1981. Determination of hydrogeological parameters using sinusoidal pressure tests: a theoretical appraisal, Water Resour. Res. , 17, 686– 692. Google Scholar CrossRef Search ADS   Bourdet D., Ayoub J., Pirard Y., 1989. Use of pressure derivative in well test interpretation, Soc. Petrol. Eng. Format. Eval. , 4, 293– 302. Caine J.S., Evans J.P., Forster C.B., 1996. Fault zone architecture and permeability structure, Geology , 24, 1025– 1028. Google Scholar CrossRef Search ADS   Cardiff M., Bakhos T., Kitanidis P.K., Barrash W., 2013. Aquifer heterogeneity characterization with oscillatory pumping: sensitivity analysis and imaging potential, Water Resour. Res. , 49, 5395– 5410. Google Scholar CrossRef Search ADS   Carslaw H.S., Jaeger J.C., 1986. Conduction of Heat in Solids , 2nd edn, pp. 510, Oxford Univ. Press. Cinco-Ley H., Samaniego-V F., 1981. Transient pressure analysis for fractured Wells, J. Petrol. Technol. , 33, 1749– 1766. Google Scholar CrossRef Search ADS   Cinco-Ley H., Samaniego-V F., Dominguez-A N., 1978. Transient pressure behavior for a well with a finite-conductivity vertical fracture, Soc. Petrol. Eng. J. , 18, 253– 264. Google Scholar CrossRef Search ADS   Economides M.J., Nolte K.G., 2000. Reservoir Stimulation , 3rd edn, pp. 856, Wiley and Sons Ltd. Fetter C.W., 2001. Applied Hydrogeology , 4th edn, pp. 598, Prentice Hall. Fisher Q.J., Knipe R.J., 2001. The permeability of faults within siliciclastic petroleum reservoirs of the North Sea and Norwegian Continental Shelf, Mar. Petrol. Geol. , 18, 1063– 1081. Google Scholar CrossRef Search ADS   Fokker P.A., Renner J., Verga F., 2012. Applications of harmonic pulse testing to field cases, in SPE Europec/EAGE Annual Conference , 4–7 June, Copenhagen, Denmark. doi:10.2118/154048-MS. Fokker P.A., Renner J., Verga F., 2013. Numerical modeling of periodic pumping tests in wells penetrating a heterogeneous aquifer, Am. J. Environ. Sci. , 9, 1– 13. Google Scholar CrossRef Search ADS   Garcia J., Hartline C., Walters M., Wright M., Rutqvist J., Dobson P.F., Jeanne P., 2016. The Northwest Geysers EGS Demonstration Project, California: Part 1: Characterization and reservoir response to injection, Geothermics , 63, 97– 119. Google Scholar CrossRef Search ADS   Genter A., Goerke X., Graff J.-J., Cuenot N., Krall G., Schindler M., Ravier G., 2010. Current status of the EGS Soultz geothermal project (France), in Proceedings World Geothermal Congress 2010 , Bali, Indonesia. Available at: https://www.geothermal-energy.org/pdf/IGAstandard/WGC/2010/3124.pdf, last accessed 6 October 2017. Guiltinan E., Becker M.W., 2015. Measuring well hydraulic connectivity in fractured bedrock using periodic slug tests, J. Hydrol. , 521, 100– 107. Google Scholar CrossRef Search ADS   Haraden J., 1992. The status of hot dry rock as an energy source, Energy , 17, 777– 786. Google Scholar CrossRef Search ADS   Hollaender F., Hammond P.S., Gringarten A.C., 2002. Harmonic testing for continuous well and reservoir monitoring, in SPE Annual Technical Conference and Exhibition , Society of Petroleum Engineers, San Antonio, Texas. Houwers M.E., Heijnen L.J., Becker A., Rijkers R., 2015. A workflow for the estimation of fault zone permeability for geothermal production: a general model applied on the Roer Valley Graben in the Netherlands, in World Geothermal Congress 2015 , Melbourne, Australia. Available at: https://pangea.stanford.edu/ERE/db/WGC/papers/WGC/2015/12072.pdf, last accessed 6 October 2017. Jacob C.E., Lohman S.W., 1952. Nonsteady flow to a well of constant drawdown in an extensive aquifer, Trans. Am. Geophys. Union , 33, 559– 569. Google Scholar CrossRef Search ADS   López D.L., Smith L., 1996. Fluid flow in fault zones: influence of hydraulic anisotropy and heterogeneity on the fluid flow and heat transfer regime, Water Resour. Res. , 32, 3227– 3235. Google Scholar CrossRef Search ADS   Lyons W.C., 2010. Chapter 4 fluid movement in waterflooded reservoirs, in Working Guide to Reservoir Engineering , pp. 241– 277, Gulf Professional Publishing. Mathias S.A., Butler A.P., 2007. Shape factors for constant-head double-packer permeameters, Water Resour. Res. , 43, W06430, doi:10.1029/2006WR005279. Matthews C.S., Russell D.G., 1967. Pressure Buildup and Flow Tests in Wells , pp. 163, Society of Petroleum Engineers. Mitchell T.M., Faulkner D.R., 2009. The nature and origin of off-fault damage surrounding strike-slip fault zones with a wide range of displacements: a field study from the Atacama fault system, northern Chile. J. Struct. Geol. , 31, 802– 816. Google Scholar CrossRef Search ADS   Ortiz R.A.E., Jung R., Renner J., 2013. Two-dimensional numerical investigations on the termination of bilinear flow in fractures, Solid Earth , 4, 331– 345. Google Scholar CrossRef Search ADS   Rabinovich A., Barrash W., Cardiff M., Hochstetler D.L., Bakhos T., Dagan G., Kitanidis P.K., 2015. Frequency dependent hydraulic properties estimated from oscillatory pumping tests in an unconfined aquifer, J. Hydrol. , 531, 2– 16. Google Scholar CrossRef Search ADS   Rasmussen T.C., Haborak K.G., Young M.H., 2003. Estimating aquifer hydraulic properties using sinusoidal pumping at the Savannah River site, South Carolina, USA, Hydrol. J. , 11, 466– 482. Rehbinder G., 1996. The double packer permeameter with narrow packers. Analytical solution for non steady flow, Appl. Sci. Res. , 56, 255– 279. Google Scholar CrossRef Search ADS   Renard P., Glenz D., Mejias M., 2009. Understanding diagnostic plots for well-test interpretation, Hydrol. J. , 17, 589– 600. Renner J., Messar M., 2006. Periodic pumping tests, Geophys. J. Int. , 167, 479– 493. Google Scholar CrossRef Search ADS   Schindler M., Nami P., Schellschmidt R., Teza D., Tischner T., 2008. Summary of hydraulic stimulation operations in the 5 km deep crystalline HDR/EGS reservoir at Soultz-sous-Forêts, in Thirty-Third Workshop on Geothermal Reservoir Engineering , 28–30 January, Stanford University, Stanford, CA, USA. Available at: https://www.geothermal-energy.org/pdf/IGAstandard/SGW/2008/schindle.pdf, last accessed 6 October 2017. Shaoul J.R., Ross M.J., Spitzer W.J., Wheaton S.R., Mayland P.J., Singh A.P., 2007. Massive hydraulic fracturing unlocks deep tight gas reserves in India, in European Formation Damage Conference , 30 May–1 June, Scheveningen, The Netherlands. doi:10.2118/107337-MS. Sheldon H.A., Micklethwaite S., 2007. Damage and permeability around faults: implications for mineralization, Geology , 35, 903– 906. Google Scholar CrossRef Search ADS   Sibson R.H., Moore J.M.M., Rankin A.H., 1975. Seismic pumping—a hydrothermal fluid transport mechanism, J. geol. Soc. , 131, 635– 659. Google Scholar CrossRef Search ADS   Song I., Renner J., 2007. Analysis of oscillatory fluid flow through rock samples, Geophys. J. Int. , 170, 195– 204. Google Scholar CrossRef Search ADS   Vinci C., Steeb H., Renner J., 2015. The imprint of hydro-mechanics of fractures in periodic pumping tests, Geophys. J. Int. , 202, 1613– 1626. Google Scholar CrossRef Search ADS   Weir G.J., 1999. Single-phase flow regimes in a discrete fracture model, Water Resour. Res. , 35, 65– 73. Google Scholar CrossRef Search ADS   Wibberley C.A.J., Shimamoto T., 2003. Internal structure and permeability of major strike-slip fault zones: the Median Tectonic Line in Mie Prefecture, Southwest Japan, J. Struct. Geol. , 25, 59– 78. Google Scholar CrossRef Search ADS   APPENDIX A: STRUCTURAL INFORMATION FOR INTERVALS Figure A1. View largeDownload slide Images of the borehole wall for the five intervals as gained from logging with an optical camera. For i2 (i5), a section of about 1 m length is also represented by a photograph of the recovered cores. Figure A1. View largeDownload slide Images of the borehole wall for the five intervals as gained from logging with an optical camera. For i2 (i5), a section of about 1 m length is also represented by a photograph of the recovered cores. APPENDIX B: MIRROR-SYMMETRIC DOUBLE-PACKER MODEL For an isotropic and homogenous subsurface, the governing diffusion equation reads in cylindrical coordinates:   \begin{equation}\frac{{{\partial ^2}p(r,z,t)}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial p(r,z,t)}}{{\partial r}} + \frac{{{\partial ^2}p(r,z,t)}}{{\partial {z^2}}} = \frac{1}{D}\frac{{\partial p(r,z,t)}}{{\partial t}}\end{equation} (B1)where r is radial distance from the centre of the borehole and z denotes the vertical direction (upright positive). Employing separation of variables, that is,   \begin{equation*}p(r,z,t) = P(r,z)\Phi (t) = P(r,z)\exp (i\omega t),\end{equation*}the spatial variation of pressure obeys   \begin{equation}\frac{{{\partial ^2}P(r,z)}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial P(r,z)}}{{\partial r}} + \frac{{{\partial ^2}P(r,z)}}{{\partial {z^2}}} = \frac{{i\omega }}{D}P(r,z).\end{equation} (B2) We consider a double-packer probe enclosing an interval with a total height of 2h0 in a well with radius ri. Instead of modeling two packers, we assume mirror symmetry to the horizontal plane bisecting the interval (Fig. B1) and thus we can restrict to treating the upper half of the model space. The locations of the bottom and top of the packer are z = h0 and z = hp, respectively, and the aquifer extends vertically to z = b (Fig. B1). We adapt boundary conditions from Rehbinder (1996) and Mathias & Butler (2007): BC1: $${\displaystyle\frac{{\partial P(r,z)}}{{\partial z}} = 0}\quad {z = b}$$ no-flow condition for vertical bounds of model BC2: $${P(r,z) = 0}\quad {r \to \infty }$$ radially infinite reservoir BC3: $${Q = 4\pi {r_{\rm i}}\int\nolimits_0^{{h_0}} {q({r_{\rm i}},z)} {\rm{d}}z\quad {\rm{ with }}\quad q({r_{\rm i}},z) = - \displaystyle\frac{T}{{2\rho g{h_0}}}{{ {\displaystyle\frac{{\partial P(r,z)}}{{\partial r}}} |}_{{r_{\rm{i}}}}}}\quad {0 \le z} \le {h_0}\quad $$ integrating flux q over the injection interval has to yield the constant flow rate Q BC4: $${P({r_{\rm i}},z) = {p_1}}\quad {{h_{\rm p}} \le z} \le b \quad $$ average pressure above packer (neglecting effect of fluid column) BC5: $${P({r_{\rm i}},z) = {p_2}}\quad {0 \le z} \le {h_0} \quad $$ average pressure in injection interval (neglecting effect of fluid column) BC6: $${P({r_{\rm i}},z) = ({p_2} - {p_1})\displaystyle\frac{{{h_{\rm p}} - z}}{{{h_{\rm p}} - {h_0}}} + {p_1}}\quad {{h_0} \le z} \le {h_{\rm p}}$$ linear pressure change along packer section where the pressures in the well above and in the interval are denoted p1 and p2, respectively. Figure B1. View largeDownload slide Model geometry for the analytical solution of a double-packer test. The model has mirror symmetry with respect to the horizontal plane that bisects the interval. Figure B1. View largeDownload slide Model geometry for the analytical solution of a double-packer test. The model has mirror symmetry with respect to the horizontal plane that bisects the interval. After Fourier cosine transforming of z, that is,   \begin{eqnarray}{{P_n}(r) = \frac{2}{b}\int\nolimits_0^b {P(r,z){\rm {cos}}({a_n}z)\,{\rm{d}}z,}}\quad {{a_n} = \frac{{n\pi }}{b}} \end{eqnarray} (B3)and accounting for BC1, eq. (B2) becomes   \begin{equation}\frac{{{\partial ^2}{P_n}(r)}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial {P_n}(r)}}{{\partial r}} - \left( {{a}_n^{2} + \frac{{i\omega }}{D}} \right){P_n}(r) = 0\end{equation} (B4)with a general solution in the form of   \begin{eqnarray}{{P_n}(r) = {A_n}{K_0}({\eta _n}r) + {B_n}{I_0}({\eta _n}r),}\quad {{\eta _n} = \sqrt {{a}_n^{2} + \frac{{i\omega }}{D}} } \end{eqnarray} (B5)where K0 and I0 denote modified Bessel functions of the first and second kind of zero order. Since I0 → ∞ when r → ∞ (BC2), we have to require Bn = 0, that is,   \begin{equation}{P_n}(r) = {A_n}{K_0}({\eta _n}r).\end{equation} (B6) The backtransform reads   \begin{equation}P\left( {r,z} \right) = \frac{{{A_0}{K_0}({\eta _0}r)}}{2} + \sum\limits_{n = 1}^\infty {{A_n}{K_0}({\eta _n}r){\rm {cos}}\left( {{a_n}z} \right)} .\end{equation} (B7) Application of BC4–BC6 gives   \begin{equation}\frac{{{A_0}{K_0}({\eta _0}{r_{\rm i}})}}{2} + \sum\limits_{n = 1}^\infty {{A_n}{K_0}({\eta _n}{r_{\rm i}}){\rm {cos}}\left( {{a_n}z} \right)} = \left\{ {\begin{array}{c@\quad c} {{p_2}}& {0 \le z \le {h_0}}\\ {({p_2} - {p_1})(\frac{{{h_{\rm p}} - z}}{{{h_{\rm p}} - {h_0}}}) + {p_1}}& {{h_0} \le z \le {h_{\rm p}}}\\ {{p_1}}& {{h_{\rm p}} \le z \le b} \end{array}} \right.\end{equation} (B8)which upon Fourier-cosine transformation of z yields   \begin{eqnarray}&&{\frac{{{A_0}{K_0}({\eta _0}{r_{\rm i}})}}{2}\underbrace {\int\nolimits_0^b {{\rm {cos}}\left( {{a_m}z} \right)}\, {\rm{d}}z}_{ = 0} + \sum\limits_{n = 1}^\infty {{A_n}{K_0}({\eta _n}{r_{\rm i}})\int\nolimits_0^b {{\rm {cos}}\left( {{a_n}z} \right){\rm {cos}}\left( {{a_m}z} \right){\rm{d}}z} }}\nonumber\\ &=& \int\nolimits_0^{{h_0}} {{p_2}{\rm {cos}}\left( {{a_m}z} \right){\rm{d}}z} + \int\nolimits_{{h_0}}^{{h_{\rm p}}} {\left[ {({p_2} - {p_1})\frac{{{h_{\rm p}} - z}}{{{h_{\rm p}} - {h_0}}} + {p_1}} \right]{\rm {cos}}\left( {{a_m}z} \right){\rm{d}}z} + \int\nolimits_{{h_{\rm p}}}^b {{p_1}{\rm {cos}}\left( {{a_m}z} \right){\rm{d}}z} .\end{eqnarray} (B9) For an − am ≠ 0, the second integrand on the left-hand side can be modified as   \begin{equation*}{\rm {cos}}\left( {{a_n}z} \right){\rm {cos}}\left( {{a_m}z} \right) = \frac{{\cos \left[ {({a_n} - {a_m})z} \right] + \cos \left[ {({a_n} + {a_m})z} \right]}}{2}\end{equation*} and thus   \begin{equation*}\int\nolimits_0^b {\displaystyle\frac{{\cos \left[ {({a_n} - {a_m})z} \right] + \cos \left[ {({a_n} + {a_m})z} \right]}}{2}{\rm{d}}z} = \left. {\left[ {\displaystyle\frac{{\sin \left[ {(n - n)\pi \frac{z}{b}} \right]}}{{2({a_n} - {a_m})}} + \frac{{\sin \left[ {(n + n)\pi \frac{z}{b}} \right]}}{{2({a_n} + {a_m})}}} \right]} \right|_0^b = 0,\end{equation*} that is, all elements of the sum vanish but the one for which an − am = 0. In this case, the left-hand side of eq. (B9) can be algebraically manipulated as   \begin{equation}{A_m}{K_0}({\eta _m}{r_{\rm i}})\int\nolimits_0^b {{\rm cos}^2}\left( {{a_m}z} \right){\rm{d}}z = {A_m}{K_0}({\eta _m}{r_{\rm i}})\frac{b}{2}.\end{equation} (B10) After integration, the right-hand side sums up to   \begin{equation}{p_2}\frac{{\cos \left( {{a_m}{h_0}} \right)}}{{\left( {{h_{\rm p}} - {h_0}} \right){a}_m^{2}}} - {p_2}\frac{{\cos \left( {{a_m}{h_{\rm p}}} \right)}}{{\left( {{h_{\rm p}} - {h_0}} \right){a}_m^{2}}} - {p_1}\left[ {\frac{{\cos \left( {{a_m}{h_0}} \right)}}{{\left( {{h_{\rm p}} - {h_0}} \right){a}_m^{2}}} - \frac{{\cos \left( {{a_m}{h_{\rm p}}} \right)}}{{\left( {{h_{\rm p}} - {h_0}} \right){a}_m^{2}}}} \right].\end{equation} (B11) With both, left- and right-hand sides, eqs (B10) and (B11), respectively, we can solve for the unknown coefficients   \begin{equation}{A_m} = \displaystyle\frac{{{p_2}\left[ {\displaystyle\frac{{\cos \left( {{a_m}{h_0}} \right)}}{{\left( {{h_{\rm p}} - {h_0}} \right){a}_m^{2}}} - \displaystyle\frac{{\cos \left( {{a_m}{h_{\rm p}}} \right)}}{{\left( {{h_{\rm p}} - {h_0}} \right){a}_m^{2}}}} \right] - {p_1}\left[ {\displaystyle\frac{{\cos \left( {{a_m}{h_0}} \right)}}{{\left( {{h_{\rm p}} - {h_0}} \right){a}_m^{2}}} - \displaystyle\frac{{\cos \left( {{a_m}{h_{\rm p}}} \right)}}{{\left( {{h_{\rm p}} - {h_0}} \right){a}_m^{2}}}} \right]}}{{\displaystyle\frac{b}{2}{K_0}({\eta _m}{r_{\rm i}})}}.\end{equation} (B12) To derive A0, we apply Darcy's law to eq. (B7) and obtain the flow rate according to (BC3), that is,   \begin{equation}Q = 2\pi {r_{\rm i}}\frac{T}{{\rho g{h_0}}}\left( {\frac{{{A_0}{\eta _0}{K_1}({\eta _0}{r_{\rm i}})}}{2}{h_0} + \sum\limits_{n = 1}^\infty {{A_n}{\eta _n}{K_1}({\eta _n}{r_{\rm i}})\frac{{\sin \left( {{a_n}{h_0}} \right)}}{{{a_n}}}} } \right).\end{equation} (B13)and thus   \begin{equation}{A_0} = \frac{{Q\rho g}}{{\pi {r_{\rm i}}T{\eta _0}{K_1}({\eta _0}{r_{\rm i}})}} - 2\sum\limits_{n = 1}^\infty {{A_n}\frac{{{\eta _n}{K_1}({\eta _n}{r_{\rm i}})}}{{{\eta _0}{K_1}({\eta _0}{r_{\rm i}})}}\frac{{\sin \left( {{a_n}{h_0}} \right)}}{{{a_n}{h_0}}}} .\end{equation} (B14) Finally, the pressure function eq. (B7) reads   \begin{equation}P\left( {r,z} \right) = \frac{{Q\rho g}}{{2\pi {r_{\rm i}}T}}\frac{{{K_0}({\eta _0}r)}}{{{\eta _0}{K_1}({\eta _0}{r_{\rm i}})}} - \frac{{{K_0}({\eta _0}r)}}{{{\eta _0}{K_1}({\eta _0}{r_{\rm i}})}}\sum\limits_{n = 1}^\infty {{A_n}{\eta _n}{K_1}({\eta _n}{r_{\rm i}})\frac{{\sin \left( {{a_n}{h_0}} \right)}}{{{a_n}{h_0}}}} + \sum\limits_{n = 1}^\infty {{A_n}{K_0}({\eta _n}r){\rm {cos}}\left( {{a_n}z} \right)} \end{equation} (B15)where $${A_n} = \displaystyle\frac{{2( {{p_2} - {p_1}} )[ {\cos ( {{a_n}{h_0}} ) - \cos ( {{a_n}{h_{\rm p}}} )} ]}}{{b{K_0}({\eta _n}{r_{\rm i}})( {{h_{\rm p}} - {h_0}} ){a}_n^{2}}}.$$ For the injection interval, with P(ri, z < h0) = p2, we now find a relation between the ratios of flow rate and interval pressure on the one hand and between interval pressure and pressure outside the interval on the other hand:   \begin{eqnarray}{\frac{Q}{{{p_2}}} = \frac{{2\pi {r_{\rm i}}T}}{{\rho g}}\left\{ {\frac{{{\eta _0}{K_1}({\eta _0}{r_{\rm i}})}}{{{K_0}({\eta _0}{r_{\rm i}})}} + 2\frac{{1 - \frac{{{p_1}}}{{{p_2}}}}}{{b\left( {{h_{\rm p}} - {h_0}} \right)}}} {\sum\limits_{n = 1}^\infty {\frac{{\cos \left( {{a_n}{h_0}} \right) - \cos \left( {{a_n}{h_{\rm p}}} \right)}}{{{K_0}({\eta _n}{r_{\rm i}}){a}_n^{2}}}\left[ {\frac{{{\eta _n}{K_1}({\eta _n}{r_{\rm i}})}}{{{h_0}{a_n}}}\sin \left( {{a_n}{h_0}} \right) - \frac{{{\eta _0}{K_0}({\eta _n}{r_{\rm i}}){K_1}({\eta _0}{r_{\rm i}})}}{{{K_0}({\eta _0}{r_{\rm i}})}}\cos\left( {{a_n}z} \right)} \right]} } \right\}} .\end{eqnarray} (B16) It can be seen that solution consists of two parts: standard radial flow for infinite line source and an additional term depending on packer positioning. When the pressure outside of the interval remains unaltered by the pumping operation, that is, p1 = p0 = const., we can—without loss of generality—chose the reference pressure to p0 = 0. In this case, eq. (B16) simplifies to   \begin{eqnarray}\frac{Q}{{{p_2}}} = \frac{{2\pi {r_{\rm i}}T}}{{\rho g}}\left\{ {\frac{{{\eta _0}{K_1}({\eta _0}{r_{\rm i}})}}{{{K_0}({\eta _0}{r_{\rm i}})}} + \frac{2}{{b\left( {{h_{\rm p}} - {h_0}} \right)}}} {\sum\limits_{n = 1}^\infty {\frac{{\cos \left( {{a_n}{h_0}} \right) - \cos \left( {{a_n}{h_{\rm p}}} \right)}}{{{K_0}({\eta _n}{r_{\rm i}}){a}_n^{2}}}\left[ {\frac{{{\eta _n}{K_1}({\eta _n}{r_{\rm i}})}}{{{h_0}{a_n}}}\sin \left( {{a_n}{h_0}} \right) - \frac{{{\eta _0}{K_0}({\eta _n}{r_{\rm i}}){K_1}({\eta _0}{r_{\rm i}})}}{{{K_0}({\eta _0}{r_{\rm i}})}}{\rm {cos}}\left( {{a_n}z} \right)} \right]} } \right\}, \end{eqnarray} (B17)an analytical expression for the ratio of flow rate and interval pressure. APPENDIX C: CONCENTRIC CYLINDRICAL SHELLS SURROUNDING THE BOREHOLE C1 Analytical solution We consider radial flow from a borehole into and beyond a single concentric cylindrical shell located between ri, the radius of the injection well and r1 (Fig. C1). The diffusion equation can be written as:   \begin{equation}\frac{{{\partial ^2}p(r,t)}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial p(r,t)}}{{\partial r}} = \frac{1}{D}\frac{{\partial p(r,t)}}{{\partial t}}.\end{equation} (C1) Figure C1. View largeDownload slide Geometry of radial flow from a well into N concentric shells and the surrounding medium. Figure C1. View largeDownload slide Geometry of radial flow from a well into N concentric shells and the surrounding medium. Here, p is fluid pressure, t is time, r is radial distance and D is hydraulic diffusivity. The hydraulic properties of the shell (subscript 1) and the surrounding medium (subscript 2) are described by permeability k1, 2 = μT1, 2/(ρgh) and diffusivity D1, 2 = k1, 2/(s1, 2μ) = T1, 2/S1, 2 where μ denotes the fluid viscosity, T1, 2, S1, 2 and s1, 2 transmissivity, storativity and specific storage capacity of the two media, respectively, g gravitational acceleration and h the length of the injection interval. Solving the diffusion equation for periodic pressure variations in the well by separation of variables (compare Renner Messar 2006) and account for boundary conditions $${{p_2}(r \ge {r_1},t) = 0}\quad {r \to \infty }\quad $$ radially infinite reservoir $${{p_1}(r,t) = {p_2}(r,t)}\quad {r = {r_1}} \quad $$ continuity of pressure at outer shell radius $${{T_1}\displaystyle\frac{{\partial {p_1}(r,t)}}{{\partial r}} = {T_2}\frac{{\partial {p_2}(r,t)}}{{\partial r}}}\quad {r = {r_{\rm{1}}}}\quad $$ continuity of flow rate at outer shell radius yields amplitude ratio and phase shift between flow rate and interval pressure of   \begin{equation}{\delta _{{\rm{qp}}}} = 2\pi {r_{\rm i}}\frac{{{T_1}}}{{\rho g}}\left| {\frac{{{\eta _1}\left( { - {m_1}{K_1}\left( {{\eta _1}{r_{\rm i}}} \right) + {I_1}\left( {{\eta _1}{r_{\rm i}}} \right)} \right)}}{{{m_1}{K_0}\left( {{\eta _1}{r_{\rm i}}} \right) + {I_0}\left( {{\eta _1}{r_{\rm i}}} \right)}}} \right|,\end{equation} (C2)and   \begin{equation}{\varphi _{{\rm{qp}}}} = {\rm{arg}}\left( {\frac{{{\eta _1}\left( { - {m_1}{K_1}\left( {{\eta _1}{r_{\rm i}}} \right) + {I_1}\left( {{\eta _1}{r_{\rm i}}} \right)} \right)}}{{{m_1}{K_0}\left( {{\eta _1}{r_{\rm i}}} \right) + {I_0}\left( {{\eta _1}{r_{\rm i}}} \right)}}} \right),\end{equation} (C3)respectively. The constant m1 is related to the ratios of permeability (or equivalently transmissivity) and diffusivity according to   \begin{equation}{m_1} = \displaystyle\frac{{{I_1}\left( {{\eta _1}{r_1}} \right){K_0}\left( {{\eta _2}{r_1}} \right) + {I_0}\left( {{\eta _1}{r_1}} \right){K_1}\left( {{\eta _2}{r_1}} \right)\sqrt {\displaystyle\frac{{{D_1}}}{{{D_2}}}} \frac{{{k_2}}}{{{k_1}}}}}{{{K_1}\left( {{\eta _1}{r_1}} \right){K_0}\left( {{\eta _2}{r_1}} \right) - {K_0}\left( {{\eta _1}{r_1}} \right){K_1}\left( {{\eta _2}{r_1}} \right)\sqrt {\displaystyle\frac{{{D_1}}}{{{D_2}}}} \frac{{{k_2}}}{{{k_1}}}}},\end{equation} (C4)where I0 and I1 represent the modified Bessel functions of the first kind of zero and first order and K0 and K1 the modified Bessel functions of the second kind of zero and first order. The arguments of these Bessel functions contain the complex parameters $${\eta _{1,2}} = \sqrt {i\omega /{D_{1,2}}} $$ that depend on angular frequency ω and diffusivity of the two media. The analytical solution for N concentric cylindrical shells (Fig. C1) is given as   \begin{equation}{\delta _{{\rm{qp}}}} = 2\pi {r_{\rm i}}\frac{{{T_1}}}{{\rho g}}\left| {\frac{{{\eta _1}\left( { - {c_1}{K_1}\left( {{\eta _1}{r_{\rm i}}} \right) + {I_1}\left( {{\eta _1}{r_{\rm i}}} \right)} \right)}}{{{c_1}{K_0}\left( {{\eta _1}{r_{\rm i}}} \right) + {I_0}\left( {{\eta _1}{r_{\rm i}}} \right)}}} \right|,\end{equation} (C5)and   \begin{equation}{\varphi _{{\rm{qp}}}} = {\rm{arg}}\left( {\frac{{{\eta _1}\left( { - {c_1}{K_1}\left( {{\eta _1}{r_{\rm i}}} \right) + {I_1}\left( {{\eta _1}{r_{\rm i}}} \right)} \right)}}{{{c_1}{K_0}\left( {{\eta _1}{r_{\rm i}}} \right) + {I_0}\left( {{\eta _1}{r_{\rm i}}} \right)}}} \right),\end{equation} (C6)and j = 1, …N−1. For N = 1, c1 = m1 but for N > 1, to obtain c1 one has to successively determine c2, c3… and cN from   \begin{equation*}{c_N} = - \displaystyle\frac{{\frac{{{k_{N + 1}}}}{{{k_N}}}\sqrt {\displaystyle\frac{{{D_N}}}{{{D_{N + 1}}}}} {K_1}({\eta _{N + 1}}{r_N}){I_0}({\eta _N}{r_N}) + {I_1}({\eta _N}{r_N}){K_0}({\eta _{N + 1}}{r_N})}}{{\displaystyle\frac{{{k_{N + 1}}}}{{{k_N}}}\sqrt {\frac{{{D_N}}}{{{D_{N + 1}}}}} {K_1}({\eta _{N + 1}}{r_N}){K_0}({\eta _N}{r_N}) - {K_1}({\eta _N}{r_N}){K_0}({\eta _{N + 1}}{r_N})}}\end{equation*}   \begin{equation*}{c_j} = - \displaystyle\frac{{{I_0}({\eta _j}{r_j})\displaystyle\frac{{{k_{j + 1}}}}{{{k_j}}}\sqrt {\frac{{{D_j}}}{{{D_{j + 1}}}}} \left( {{c_{j + 1}}{K_1}({\eta _{j + 1}}{r_j}) - {I_1}({\eta _{j + 1}}{r_j})} \right) + {I_1}({\eta _j}{r_j})\left( {{c_{j + 1}}{K_0}({\eta _{j + 1}}{r_j}) + {I_0}({\eta _{j + 1}}{r_j})} \right)}}{{{K_0}({\eta _j}{r_j})\displaystyle\frac{{{k_{j + 1}}}}{{{k_j}}}\sqrt {\displaystyle\frac{{{D_j}}}{{{D_{j + 1}}}}} \left( {{c_{j + 1}}{K_1}({\eta _{j + 1}}{r_j}) - {I_1}({\eta _{j + 1}}{r_j})} \right) - {K_1}({\eta _j}{r_j})\left( {{c_{j + 1}}{K_0}({\eta _{j + 1}}{r_j}) + {I_0}({\eta _{j + 1}}{r_j})} \right)}}\end{equation*} C2 Parameter study For the parameter study, we restrict to the simple case of a single shell. We investigated the sensitivity of phase shift to the various model parameters by systematic parameter variation. Phase shift is reported as a function of $${r_{\rm i}}\sqrt {\omega /{D_{{\rm{shell}}}}} $$ where Dshell = D1 denotes the hydraulic diffusivity of the shell (Fig. C1). This expression can be understood either as a dimensionless frequency or an inverse dimensionless penetration depth given in multiples of the borehole radius. A sigmoidal shape of the relation between phase shift and dimensionless frequency is characteristic for a homogeneous medium. Small phase shifts occur at low frequencies and vice versa with a rather steep switch between the two limits around a dimensionless frequency of 1. This shape is significantly perturbed when ratios between the parameters of shell and medium between one-tenth and a factor of 10 are allowed. The high-frequency limit of 1/8 for phase shift is not strictly valid for the compound model but values as high as 1/5 are found. The phase shift for a shell with a low storage property compared to the medium exhibits a pronounced maximum that forms a shoulder on the steep increase characteristic for the homogeneous medium (Fig. C2a). The location of the maximum moves towards lower frequencies with increasing radius of the shell. A variation of the ratio of storage parameters for shell and medium for a fixed shell radius results in qualitatively similar phase-shift curves as the variation in shell radius (Fig. C2b). Yet, the maximum in phase transforms into a minimum preceded by a new maximum at lower frequency when the storage capacity in the shell switches from lower to higher than the storage capacity in the medium. Figure C2. View largeDownload slide Results of the shell model given with phase shift as a function of dimensionless frequency (a) for various radius ratios and indicated fixed permeability (transmissivity) and diffusivity ratios, (b) for various diffusivity ratios and indicated fixed permeability (transmissivity) and radius ratios, (c) for various permeability (transmissivity) ratios when storativity is covaried such that the diffusivity of shell and medium are identical and (d) for various identical diffusivity and permeability (transmissivity) ratios that correspond to an identical storativity for shell and medium at indicated fixed radius ratio. For identical permeability (transmissivity) of shell and medium, the diffusivity ratio is related to an inverse ratio of specific storage capacity (storativity), that is, D1/D2 = s2/s1 = S2/S1. Figure C2. View largeDownload slide Results of the shell model given with phase shift as a function of dimensionless frequency (a) for various radius ratios and indicated fixed permeability (transmissivity) and diffusivity ratios, (b) for various diffusivity ratios and indicated fixed permeability (transmissivity) and radius ratios, (c) for various permeability (transmissivity) ratios when storativity is covaried such that the diffusivity of shell and medium are identical and (d) for various identical diffusivity and permeability (transmissivity) ratios that correspond to an identical storativity for shell and medium at indicated fixed radius ratio. For identical permeability (transmissivity) of shell and medium, the diffusivity ratio is related to an inverse ratio of specific storage capacity (storativity), that is, D1/D2 = s2/s1 = S2/S1. When permeability (transmissivity) and specific storage capacity (storativity) are modified in the same way, that is, diffusivity is kept constant, the transition from the low- to the high-frequency regime is significantly perturbed. Low permeability and storativity in the shell relative to the surrounding medium leads to a transition that resembles a Heavyside function (Fig. C2c). High permeability and storativity in the shell gives rise to a pronounced maximum in phase shift at a dimensionless frequency of about 0.1 that reaches values as large as 1/5, that is, well above the limit for a homogeneous medium of 1/8. Covarying peremability and diffusivity and keeping the storage parameter uniform for shell and medium has qualitatively similar effects (Fig. C2d). APPENDIX D: BILINEAR FLOW MODEL The model consists of a slit with distinct hydraulic properties embedded in a homogeneous matrix and assumes that the direct volume flow from the source into the matrix along the x-direction is negligible (Fig. D1), but the fracture boundary is considered a harmonic source for flow into the matrix. The governing equations in the matrix and the fracture correspond to diffusion equations:   \begin{equation}\frac{{\partial {p_{\rm m}}(x,y,t)}}{{\partial t}} = {D_{\rm m}}{\nabla ^2}{p_{\rm m}}(x,y,t)\end{equation} (D1)  \begin{equation}\frac{{\partial {p_{\rm f}}(x,t)}}{{\partial t}} = {D_{\rm f}}\frac{{{\partial ^2}{p_{\rm f}}(x,t)}}{{\partial {x^2}}} + \frac{{2{D_{\rm m}}}}{\delta }{\left. {\frac{{\partial {p_{\rm m}}(x,y,t)}}{{\partial y}}} \right|_{\left| y \right| = \delta /2}}\end{equation} (D2)where pm(x, y, t) is the pressure in the matrix, pf(x, t) the pressure in the fracture, δ the width of the fracture and Dm and Df the diffusivity of matrix and fracture, respectively. Figure D1. View largeDownload slide Model for bilinear flow: the fracture boundary is a source for flow into the matrix in y-direction; no direct volume flow from the source into the matrix along the x-direction. Figure D1. View largeDownload slide Model for bilinear flow: the fracture boundary is a source for flow into the matrix in y-direction; no direct volume flow from the source into the matrix along the x-direction. As the direct volume flow from the source into the matrix along the x-direction is negligible, ∂2pm(x, y, t)/∂x2 ≈ 0, so that eq. (D1) reduces to   \begin{equation}\displaystyle\frac{{\partial {p_{\rm m}}(x,y,t)}}{{\partial t}} = {D_{\rm m}}\displaystyle\frac{{{\partial ^2}{p_{\rm m}}(x,y,t)}}{{\partial {y^2}}}.\end{equation} (D3) The pressure in the matrix, that is, the solution of eq. (D1), is given by the standard solution for diffusion from a harmonic source into a homogeneous semi-infinite medium (Carslaw & Jaeger 1986):   \begin{equation}{p_{\rm m}}\left( {x,y,t} \right) = \displaystyle\frac{{\left| y \right| - \delta /2}}{{2\sqrt {\pi {D_{\rm m}}} }}\mathop \int \nolimits_0^t {p_{\rm f}}\left( {x,\lambda } \right)\displaystyle\frac{{\exp \left[ { - \displaystyle\frac{{{{(\left| y \right| - \delta /2)}^2}}}{{4{D_{\rm m}}\left( {t - \lambda } \right)}}} \right]}}{{{{\left( {t - \lambda } \right)}^{3/2}}}}{\rm{d}}\lambda \end{equation} (D4)whose insertion into eq. (D2) yields   \begin{equation}{D_{\rm f}}\displaystyle\frac{{{\partial ^2}{p_{\rm f}}\left( {x,t} \right)}}{{\partial {x^2}}} + \displaystyle\frac{{2{D_{\rm m}}}}{\delta }\displaystyle\frac{\partial }{{\partial y}}{\left[ {\displaystyle\frac{{\left| y \right| - \delta /2}}{{2\sqrt {\pi {D_{\rm m}}} }}\mathop \int \nolimits_0^t {p_{\rm f}}\left( {x,\lambda } \right)\displaystyle\frac{{\exp \left[ { - \displaystyle\frac{{{{(\left| y \right| - \delta /2)}^2}}}{{4{D_{\rm m}}\left( {t - \lambda } \right)}}} \right]}}{{{{\left( {t - \lambda } \right)}^{3/2}}}}{\rm{d}}\lambda } \right]_{\left| y \right| - \delta /2 = 0}} = \displaystyle\frac{{\partial {p_{\rm f}}\left( {x,t} \right)}}{{\partial t}}.\end{equation} (D5) Differentiation of the second term by part restricting to the section y ≥ 0 because of the symmetry of the problem gives   \begin{equation}{D_{\rm f}}\displaystyle\frac{{{\partial ^2}{p_{\rm f}}(x,t)}}{{\partial {x^2}}} + \displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\sqrt \pi \delta }}{\left[ {\int\nolimits_0^t {{p_{\rm f}}(x,\lambda )\left( {1 + \displaystyle\frac{{{{(y - {\delta / 2})}^2}}}{{2{D_{\rm m}}(t - \lambda )}}} \right)\displaystyle\frac{{\exp \left[ { - \displaystyle\frac{{{{(\left| y \right| - {\delta / 2})}^2}}}{{4{D_{\rm m}}(t - \lambda )}}} \right]}}{{{{\left( {t - \lambda } \right)}^{3/2}}}}} {\rm{d}}\lambda } \right]_{y = {\delta / 2}}} = \displaystyle\frac{{\partial {p_{\rm f}}(x,t)}}{{\partial t}}.\end{equation} (D6) Applying y = δ/2 we get   \begin{equation}{D_{\rm f}}\displaystyle\frac{{{\partial ^2}{p_{\rm f}}(x,t)}}{{\partial {x^2}}} + \displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\sqrt \pi \delta }}\int\nolimits_0^t {\displaystyle\frac{{{p_{\rm f}}(x,\lambda )}}{{{{\left( {t - \lambda } \right)}^{\frac{3}{2}}}}}} {\rm{d}}\lambda = \displaystyle\frac{{\partial {p_{\rm f}}(x,t)}}{{\partial t}}\end{equation} (D7)and thus for harmonic pressure variations in the fracture   \begin{equation}{D_{\rm f}}\displaystyle\frac{{{\partial ^2}{p_{\rm f}}(x)}}{{\partial {x^2}}}A\exp (i\omega t) + \displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\sqrt \pi \delta }}{p_{\rm f}}(x)\int\nolimits_0^t {\displaystyle\frac{{A\exp (i\omega \lambda )}}{{{{\left( {t - \lambda } \right)}^{\frac{3}{2}}}}}} {\rm{d}}\lambda = i\omega A\exp (i\omega t){p_{\rm f}}(x).\end{equation} (D8) The integral solution is   \begin{equation}\int\nolimits_0^t {\displaystyle\frac{{A\exp (i\omega \lambda )}}{{{{\left( {t - \lambda } \right)}^{\frac{3}{2}}}}}{\rm{d}}\lambda } = \gamma \left( { - 0.5,i\omega t} \right)\exp (i\omega t) = \left[ {\Gamma \left( { - 0.5} \right) - \Gamma \left( { - 0.5,i\omega t} \right)} \right]\sqrt {i\omega } \exp (i\omega t)\end{equation} (D9)where γ( − 0.5, iωt), Γ( − 0.5) and Γ( − 0.5, iωt) represent the lower incomplete gamma function, the gamma function and the upper incomplete gamma function, respectively. For times longer than 1.25 times the period, the deviation of γ( − 0.5, iωt) from Γ( − 0.5) is less than 1 per cent and beyond 2.2 times the period the deviation falls below 0.5 per cent. Thus, we can approximate the integral as   \begin{equation}\mathop \int \nolimits_0^t \displaystyle\frac{{\exp \left( {i\omega \lambda } \right)}}{{{{\left( {t - \lambda } \right)}^{3/2}}}}{\rm{d}}\lambda \cong {\rm{\Gamma }}\left( { - 0.5} \right)\sqrt {i\omega } \exp \left( {i\omega t} \right) = - 2\sqrt {i\pi \omega } \exp \left( {i\omega t} \right)\end{equation} (D10)to arrive at   \begin{equation}{D_{\rm f}}\displaystyle\frac{{{\partial ^2}{p_{\rm f}}\left( x \right)}}{{\partial {x^2}}} - \displaystyle\frac{2}{\delta }\sqrt {i{D_{\rm m}}\omega } {p_{\rm f}}\left( x \right) = i\omega {p_{\rm f}}\left( x \right).\end{equation} (D11)when canceling the harmonic time dependence. Eq. (D11) is a homogeneous second-order differential equation with the solution   \begin{equation}{p_{\rm f}}\left( x \right) = c\exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } {\rm{ + }}i\omega } \right)} \right]}^{1/2}}x} \right\} + {\rm{d}}H\left( L \right)\end{equation} (D12)subject to the boundary conditions pf(x)→ 0 as x → ∞ (infinite fracture), or pf(x) = p0 as x = L (finite fracture with constant-pressure boundary), or $${{\partial {p_{\rm f}}( x )} / {\partial x}} = 0$$ as x = L (finite fracture with no-flow boundary) where c is a constant, H(L) is a parameter defined by   \begin{equation} H\left( L \right) = \left\{ {\begin{array} {l@\quad l} 0 & L \to \infty\\ 1 & L\ {\rm{finite}} \end{array}} \right.\end{equation} (D13)and   \begin{equation}d = \left\{ {\begin{array}{l} {{p_0} - c\exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i\omega } \right)} \right]}^{1/2}}L} \right\}{\rm{ for\,\, constant \hbox{-} pressure\,\, boundary}}}\\\\ {{{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i\omega } \right)} \right]}^{1/2}}cx\exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i\omega } \right)} \right]}^{1/2}}L} \right\}\,{\rm{ +\, const.}}\,\,{\rm{for\,\, no \hbox{-} flow\,\, boundary }}} \end{array}} \right.\end{equation} (D14) To simplify the calculation, we assume the constant appearing for the no-flow boundary to equal zero without loss of generality because it represents a constant shift in pressure that can be accommodated by scale shifting. Using Darcy's law, we get the flow rate as   \begin{eqnarray} {Q_{\rm f}}\left( {x,t} \right) = 2\pi {r_{\rm i}}hq\left( {x,t} \right) &=& - 2\pi {r_{\rm i}}\displaystyle\frac{T}{{\rho g}}\left( {\displaystyle\frac{{\partial {p_{\rm f}}\left( {x,t} \right)}}{{\partial x}}} \right)\nonumber\\ &=& 2\pi {r_{\rm i}}\displaystyle\frac{{TcA}}{{\rho g}}{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i\omega } \right)} \right]^{1/2}}\nonumber\\ && \cdot \left( {\exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } {\rm{ + }}i\omega } \right)} \right]}^{1/2}}x} \right\} - {H_{{\rm{bc}}}}\exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i\omega } \right)} \right]}^{1/2}}L} \right\}} \right){\rm{exp}}\left( {i\omega t} \right) \end{eqnarray} (D15)where T is transmissivity, ρ is fluid density and h is the height of the injection zone. We introduce a parameter reflecting the boundary condition at the fracture tip as   \begin{equation} {H_{{\rm{bc}}}} = \left\{ {\begin{array}{l} {0\ \ \ {\rm{constant\,\, pressure\,\, boundary }}\ }\\\\ {1\ \ \ {\rm{no\,\, flow\,\, boundary }}\ } \end{array}} \right.. \end{equation} (D16) Finally, the complex injectivity results to   \begin{eqnarray} \displaystyle\frac{{{Q_{\rm f}}\left( {x,t} \right)}}{{{p_{\rm f}}\left( {x,t} \right)}} &=& 2\pi {r_{\rm i}}\displaystyle\frac{{Tc}}{{\rho g}}{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i{\rm{\omega }}} \right)} \right]^{1/2}}\nonumber\\ &&{\rm{ }} \cdot \displaystyle\frac{{\exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i\omega } \right)} \right]}^{1/2}}x} \right\} - {H_{{\rm{bc}}}}\exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i{\rm{\omega }}} \right)} \right]}^{1/2}}L} \right\}}}{{c\exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i\omega } \right)} \right]}^{1/2}}x} \right\} + {\rm{d}}H\left( L \right)}}. \end{eqnarray} (D17) Specifically, for an infinite fracture, we find   \begin{equation}{\left. {\displaystyle\frac{{{Q_{\rm f}}(x = 0,t)}}{{{p_{\rm f}}(x = 0,t)}}} \right|_{H\left( L \right) = 0}} = 2\pi {r_{\rm i}}\displaystyle\frac{T}{{\rho g}}{\left( {\displaystyle\frac{{2\sqrt {i\pi {D_{\rm m}}\omega } + i\omega \sqrt {{D_{\rm m}}} \ }}{{\sqrt \pi \delta {D_{\rm f}}}}} \right)^{1/2}}\end{equation} (D18)with amplitude ratio and phase shift of   \begin{equation}{\delta _{{\rm{qp}}}} = 2\pi {r_{\rm i}}\displaystyle\frac{T}{{\rho g}}{\left( {\displaystyle\frac{{2\sqrt {\pi {D_{\rm m}}\omega } + \omega \sqrt {{D_{\rm m}}} \ }}{{\delta {D_{\rm f}}}}} \right)^{1/2}}\end{equation} (D19)and   \begin{equation}{\varphi _{{\rm{qp}}}} = \arg \left( {\sqrt {2\sqrt {i\pi \omega } + i\omega } } \right),\end{equation} (D20)respectively, at the injection point x = 0 (Fig. D2a). Figure D2. View largeDownload slide Phase shift for (a) constant-pressure boundary and (b) no-flow boundary at the tip of a fracture with finite length L and width δ = 0.096 m. (I) Dm = 1 m2 s−1, Df = 10 m2 s−1 and L = 20δ, (II) Dm = 1 m2 s−1, Df = 100 m2 s−1 and L = 20δ, (III) Dm = 1 m2 s−1, Df = 100 m2 s−1 and L = 100δ and (IV) Dm = 1 m2 s−1, Df = 10 m2 s−1 and L = 4δ. The relation for an infinite fracture is shown for comparison in (a). Figure D2. View largeDownload slide Phase shift for (a) constant-pressure boundary and (b) no-flow boundary at the tip of a fracture with finite length L and width δ = 0.096 m. (I) Dm = 1 m2 s−1, Df = 10 m2 s−1 and L = 20δ, (II) Dm = 1 m2 s−1, Df = 100 m2 s−1 and L = 20δ, (III) Dm = 1 m2 s−1, Df = 100 m2 s−1 and L = 100δ and (IV) Dm = 1 m2 s−1, Df = 10 m2 s−1 and L = 4δ. The relation for an infinite fracture is shown for comparison in (a). For a fracture with finite length and a constant-pressure boundary at its tip and assuming p0 = 0, we find   \begin{equation}\displaystyle\frac{{{Q_{\rm f}}(x = 0,t)}}{{{p_{\rm f}}(x = 0,t)}} = \displaystyle\frac{{{{\left. {\displaystyle\frac{{{Q_{\rm f}}(x = 0,t)}}{{{p_{\rm f}}(x = 0,t)}}} \right|}_{H\left( L \right) = 0}}}}{{1 - \exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i\omega } \right)} \right]}^{1/2}}L} \right\}}}\end{equation} (D21)at x = 0 (Fig. D2a). Analogously, we find   \begin{equation}\displaystyle\frac{{{Q_{\rm f}}(x = 0,t)}}{{{p_{\rm f}}(x = 0,t)}} = {\left. {\displaystyle\frac{{{Q_{\rm f}}(x = 0,t)}}{{{p_{\rm f}}(x = 0,t)}}} \right|_{H\left( L \right) = 0}}\left( {1 - \exp \left\{ { - {{\left[ {\displaystyle\frac{{\sqrt {{D_{\rm m}}} }}{{\delta {D_{\rm f}}\sqrt \pi }}\left( {2\sqrt {i\pi \omega } + i{\rm{\omega }}} \right)} \right]}^{1/2}}L} \right\}} \right)\end{equation} (D22)at x = 0 (Fig. D2b). © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Journal

Geophysical Journal InternationalOxford University Press

Published: Jan 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

Monthly Plan

  • Read unlimited articles
  • Personalized recommendations
  • No expiration
  • Print 20 pages per month
  • 20% off on PDF purchases
  • Organize your research
  • Get updates on your journals and topic searches

$49/month

Start Free Trial

14-day Free Trial

Best Deal — 39% off

Annual Plan

  • All the features of the Professional Plan, but for 39% off!
  • Billed annually
  • No expiration
  • For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588

$360/year

billed annually
Start Free Trial

14-day Free Trial