Crack Features and Shear-Wave Splitting Associated with Fracture Extension during Hydraulic Stimulation of the Geothermal Reservoir in Soultz-sous-Forêts

The recent tomography results obtained within the scope of the Enhanced Geothermal System (EGS) European Soultz project led us to revisit the meso-fracturing properties of Soultz test site. In this paper, we develop a novel approach coupling effective medium modeling and shear-wave splitting to characterize the evolution of crack properties throughout the hydraulic stimulation process. The stimulation experiment performed in 2000 consisted of 3 successive injection steps spanning over 6 days. An accurate 4-D tomographic image was first carried out based upon the travel-times measured for the induced seismicity [Calo M., Dorbath C., Cornet F.H., Cuenot N. (2011) Large-scale aseismic motion identified through 4-D P-wave tomography, Geophys. J. Int. 186 , 1295-1314]. The current study shows how to take advantage of the resulting compressional wave (Calo et al., 2011) and shear-wave velocity models. These are given as input data to an anisotropic effective medium model and converted into crack properties. In short, the effective medium model aims to estimate the impact of cracks on velocities. It refers to a crack-free matrix and 2 families of penny-shaped cracks with orientations in agreement with the main observed geological features: North-South strike and dip of 65°East and 65°West [Genter A., Traineau H. (1996) Analysis of macroscopic fractures in granite in the HDR geothermal well EPS-1, Soultz-sous-Forets, France, J. Vol. Geoth. Res. 72 , 121-141], respectively. The resulting output data are the spatial distributions of crack features (lengths and apertures) within the 3-D geological formation. We point out that a flow rate increase results in a crack shortening in the area imaged by both compressional and shear waves, especially in the upper part of the reservoir. Conversely, the crack length, estimated during continuous injection rate phases, is higher than during the increasing injection rate phases. A possible explanation for this is that cracks remain large because the system has time to relax. We also calculate the extension and opening rates during all hydraulic stimulation sets. While the opening rate is unchanged, the extension rate varies depending on the stimulation phase. It is also shown to be higher around and above the open-hole section than below. This can indicate a potential upward path that makes fluid percolation easier within the granite formation, this path being induced by the temperature gradient. We also compare the evolution of crack extension during injection with shear-wave splitting. Split shear waves were recorded at 2 stations during hydraulic stimulation and processed in terms of splitting parameters. The fast shear-wave polarization remains constant and parallel to the maximum horizontal stress orientation while the amplitude of splitting varies with time. We observe a good agreement between travel-time differences and crack extension rates during the first 4 days of the stimulation experiment. Afterwards, these two parameters depart from each other. This study emphasizes the added value of the coupling between effective medium modeling and shear-wave splitting to monitor meso-scale cracks in reservoirs submitted to hydraulic stimulation.


INTRODUCTION
The European Enhanced Geothermal System (EGS) project at Soultz-sous-Forêts (France) started in 1987. Its main objective was the extraction of energy from hot fractured granites. The test site is located in France on the western edge of the Rhine graben, some 50 km north of Strasbourg close to the German border. Large-scale faults are observed with distinct major orientations that can be characterized from outcrop measurements. The main orientation is known as the Rhenish orientation: it corresponds to normal faults striking N10°-N20°. These faults were created due to the development of the graben during the Oligocene. The analysis of outcrop observations is consistent with the imagery logs collected for different wells.
The concept of a geothermal system consists in extracting the heat from the rock mass by circulating fluids between injectors and producers. With the natural permeability of granite being too low, it was enhanced through hydraulic stimulation. The objective was to create a network of connected and permeable fractures between the wells. Two wells, named GPK1 and GPK2, were first drilled respectively in 1987 and 1995, down to depths of about 3 600 and 3 900 m. The depths reported in this paper are the True Vertical Depths (TVD) measured from the drilling platform. At these depths, the bottom hole temperature was 160°C. In 1999, the GPK2 well was deepened to about 5 000 m, increasing the bottom hole temperature to 200°C. Then, in 2003 and2004, two deviated wells reaching depths of 5 000 m each, GPK3 and GPK4, were drilled from the same platform as GPK2. The GPK2 well was subjected to a stimulation program in 2000 to enhance the rock mass permeability. Water was injected with an increasing stepwise rate over a period of 6 days and the well was then shut in (Fig. 1). The well head pressure was monitored all along the stimulation process plus a few days after shut-in.
Many geophysical studies in reservoir characterization focus on the variations in the elastic properties of rocks. They commonly involve seismic data, which are processed in terms of seismic attributes: amplitudes, time-shifts, velocities, impedances, etc. These processed data still have to be related to the physical properties of the rock mass and the fluids saturating the pore space. This is essential to understand the impact of stimulation on the rock system. Such need motivated the development of research projects referring to the Effective Medium Theory (EMT).
The EMT makes it possible to evaluate the macroscopic properties of a medium within a representative elementary volume provided that its microstructural properties are known. In particular, the EMT has focussed on fractured media (Kachanov, 1993;Sayers and Kachanov, 1991;Schoenberg and Sayers, 1995), where fracturing is a major factor affecting the velocities of the seismic waves propagating in the rock mass. In short, velocity variations can be explained in terms of variations in crack properties: crack density, aspect ratio or orientation. A very special characteristic of the Soultz site is the intense alteration of the granitic basement due to fluid flows. This has led to the identification of various fracture networks depending on their contributions to fluid flow within the reservoir (Sausse and Genter, 2005;Sausse et al., 2006Sausse et al., , 2010Dezayes et al., 2005Dezayes et al., , 2010. At the first order, the flow behavior of the fracture networks in granite is driven by the connectivity between fractures. This parameter is very sensitive to the fracture density and fracture geometry. The inverse problem consists in deriving information about the microstructure from the available seismic attributes. This can be achieved by integrating the forward effective model into an optimization process. Such idea was developed in active geological areas by Zhao and Mizuno (1999) and Mishra and Zhao (2003) who characterized the hypocenter locations of the Kobe (1995) and Buhj (2001) earthquakes in terms of crack density and saturation rate. On the other hand, Adelinet et al. (2011) focused on passive microseismicity data: they estimated the spatial distributions of the crack density and the fluid bulk modulus beneath the Reykjanes Peninsula thanks to accurate tomography data.
In this paper, we refer to EMT to analyze how hydraulic stimulations influence fracturing in the Soultz granitic reservoir. The leading idea is to take advantage of the available velocities to characterize the evolution of meso-crack features during stimulation. In other words, we develop an anistropic effective medium model that provides numerical responses for seismic velocities. These responses are then Wellhead injection rate (blue) and overpressure (red) in GPK2 during the 2000 stimulation experiment. The three main steps are tagged by A, B and C. The 14 steps used for the 4-D tomography, according to Calò et al. (2011), are annotated above the two curves with black italic numbers.
compared to the velocities obtained from tomography. The overall procedure consists in adjusting the input parameters of the EMT model, which are crack densities and crack apertures, to minimize the velocity mismatch. We then address the following questions: can seismic information be used to infer crack features (length and aperture)? Do cracks remain the same during water injection?
The first section of this paper recaps how EMT modeling is used to infer crack features from the data collected during the hydraulic stimulation of GPK2. The second section focuses on the analysis of the shear-wave splitting evidenced at 2 stations close to GPK2 in terms of anisotropy (amplitude and orientation). The original contribution of this paper is the subject of the third section where the EMT results are compared to the information derived from the shear-wave splitting analysis. Finally, it makes it possible to establish a physically-based model to explain the reservoir behavior during hydraulic stimulation.

Hydraulic Parameters and 4-D Tomography Data
The stimulation test performed in 2000 at well GPK2 aimed to enhance the connectivity between the open section of GPK2 and the surrounding fractured medium. Well GPK2 was open between 4 403 and 4 955 m over a section of about 500 m length. The seismic monitoring performed during this test includes two distinct periods: monitoring during injection over 6 days and monitoring during shut over 4 days after. The first period of time encompasses 3 consecutive steps during which the injection rate was first set to 30 L/s for 24 h (indicated as injection A in Fig. 1), then to 40 L/s for 27 h (injection B) and last to 50 L/s for 90 h (injection C). The total volume of fluid injected was about 23 000 m 3 . The injection rates and pressures collected at well GPK2 during stimulation are reported in Figure 1. At the very beginning, the overpressure increases very quickly to 11 MPa (bottom plot of Fig. 1). Then, it goes on increasing, but moderately, to about 12 MPa, before declining. The pressure behavior observed during the second injection step is very similar: a fast rise to 12 MPa followed by a gradual drop. The third step led to a different pressure response. The pressure again rises up at once to 12 MPa and then, it keeps on increasing smoothly till the well shut in, followed by a slow pressure decrease. Dorbath et al. (2009) concluded that the rock mass around well GPK2 could be still considered at this stage as quasi-closed. In other words, the stimulation test did not succeed in connecting the well to a major fracture network.
Intense microseismic activity was monitored during and after the stimulation test from a dense seismic network consisting of 3 down-hole array and 14 surface stations (Fig. 2) (Cuenot et al., 2008). This makes it possible to perform a 4-D tomography study from about 7 000 well-located events ( Fig. 2) which had duration magnitude from À0.9 to 2.5. Calò et al. (2011) reported the results for the P-wave tomography. The shear-wave velocities used in the present study are presented in Appendix A. The total set of events was split into 14 subsets based upon the time injection scheme, according to Calò et al. (2011). These subsets are indicated by the black italic numbers 1 to 14 in Figure 1. The corresponding periods of time have distinct durations (Tab. 1). This choice was preferred to better capture the influence of  The 2000 stimulation of GPK2: horizontal projection of the cloud of induced seismic events. The radius of each circle is proportional to the magnitude. The thick blue line indicates the open-hole section of the injection well. Surface seismic stations are identified by stars, the red ones are used for the shearwave splitting study.  (Zhang and Thurber, 2003) has been used to solve event locations and velocity structure simultaneously by taking into account both absolute and differential traveltime data. The initial 1-D reference velocity models used in the tomography workflow are presented in Table 2. The tomography results include velocity maps at depths between 4 000 and 5 000 m and centered on well GPK1 for every subsets. The initial horizontal nodes spacing for the tomography is 250 m in the X and Y directions near the center of the grid (Calò et al., 2011). After coupling a tomographic method (tomoDD) and an averaging method (WAM, described in Calò, 2009), the final grid has a resolution of 100 m in the three space directions (Calò et al., 2011). Whatever the subset under consideration, the maximal seismic velocity standard deviation is less than 0.03 km/s over the resolved grid. The resolved area depends on the seismic velocity models and subset studied: it is larger for compressional-wave than for shear-wave velocities (Fig. 3). We focus hereafter on the area jointly resolved in compressional and shear-wave velocities. Calò et al. (2011) observed a low compressional-wave velocity anomaly located around the GPK2 open-hole section with a magnitude of about 10% of the mean velocity. The anomaly is detected every time the injection rate is increased. Otherwise, during the periods of constant injection rates, the anomaly almost vanishes.

Log Data Available at Well GPK2
The GPK2 well was drilled in 1995 down to 3 883 m. Two preliminary hydraulic stimulation tests were carried out in 1995 and 1996. Then, GPK2 was deepened in 1999 to 4 955 m before performing the 2000 massive hydraulic stimulation. The casing of the well was installed to a depth of 4 403 m, meaning that the open-hole section was between 4 403 and 4 955 m. Sonic and Gamma-Ray logs were recorded between 1 500 and 3 720 m depth after the first drilling. Deeper, the well crosses large cavities. They have been explained as the intersection with a major kilometric fault (Sausse et al., 2008). There is no sonic log below 3 900 m due to a casing collapse close to a significant cavity. An initial set of Ultrasonic Borehole Imager (UBI) images collected at the end of the first drilling phase provides a description of the fracture pattern along the well for depths between 1 422 and 3 807 m. During the deepening phase, a second UBI set was acquired to image the borehole   Resolved area for velocities from tomography investigation at 4 600 m depth for both compressional (blue) and shear (red) waves. The area resolved for compressional waves is the same than the one presented in Calò et al. (2011). GPK1 wellhead is the geographical origin for the gridpoints. Axes X and Y correspond to geographical distances expressed in km. North direction is vertical upward. Black squares: GPK1 and GPK2 locations. between 3 479 and 3 866 m. The analysis of these UBI data sets reveals two crack families: natural fractures and Drilling Induced Tensile (DIT) fractures (Bérard and Cornet, 2003). As a result, the crack density estimated from UBI data collected in GPK2 cannot be directly used to quantify the natural fracturing within the reservoir. UBI data were also collected in deeper wells GPK3 and GPK4.

Theoretical Background
A challenging issue in many geophysical studies is to express the variations in the measured responses (e.g., seismic travel time or waveform amplitude) in terms of variations in the physical properties of rocks. In the case studied in this paper, this consists in explaining the evolution of seismic velocities with time against the evolution of fracture properties with time. The link is ensured referring to effective medium modeling. The EMT relies on the Representative Elementary Volume (REV), i.e. the smallest volume over which a property is considered as a value representative for the whole. The effective medium model used hereafter is a homogenization approach based on Eshelby's inclusion theory (Eshelby, 1957). It involves an ellipsoidal inclusion embedded within an infinite medium. The final objective is to evaluate the disruption in the strain field in the medium due to this inclusion. Following Eshelby (1957), the strain tensor e i in inclusion i can be expressed as: C i and C 0 are the elastic tensors for the inclusion and matrix, respectively. I is the four-order identity tensor, P the polarization Hill's tensor, and E inf the stress field far away from the inclusion. If the inclusion is flattened, as in cracks for instance, parts of the four-order tensor (I + P : (C i -C 0 ) approach infinity. In this case, the relevant tensor to describe strain is (Barthélémy, 2005): where n is the crack aspect ratio. In this study, we consider "penny-shaped" inclusions to approximate cracks and assume that these cracks are embedded into an isotropic matrix. Cracks are filled with water whose bulk modulus is 2.2 GPa. In addition, as there are many cracks within the matrix, we refer to homogenization to derive the macroscopic elasticity at the REV scale. A number of techniques have been studied in the literature as the diluted scheme, the self-consistent scheme (Budiansky and O'Connell, 1976) or the differential effective modeling (Henyey and Pomphrey, 1982;Norris, 1985). As granite has a very low porosity resulting from cracks only Guéguen and Schubnel, 2003), we chose to use the Mori-Tanaka approach (Mori and Tanaka, 1973) in which we consider cracks embedded into a matrix with the elastic properties of an inclusion-free medium. Such an approach is appropriate for rocks characterized by low porosities (< 10%) and clearly identified matrices. Thus, this approach is well-suited for estimating the elastic properties of the Soultz granite. Within this framework, cracks are described by two parameters: volumetric density q and aspect ratio n. The crack density is defined in a volumetric way following Walsh (1965) and Budiansky and O'Connell (1976): V is the total volume of the REV and c i the radius of the i-th penny-shaped crack. This volumetric density is related to the P 32 fracture intensity (ratio between area of fractures and volume of rock mass) and the mean radius of cracks, noted c, by: The elastic moduli being known at the REV scale, the compressional and shear-wave velocities can be derived from Christoffel's equation: with C ik the symetrical Christoffel's matrix, d ik the Kronecker symbol, l the matrix density and V model the velocity vector. The subscripts i and k correspond to the three main directions of space (ranging from 1 to 3). Solving this equation provides three real and positive eigenvalues, which are the three compressional-wave, slow shear-wave and fast shear-wave velocities. They obviously depend on the direction of propagation. In our model, we account for the wave paths derived from the tomography study and consider a propagation cone with a maximum deviation angle of 40°with respect to the vertical axis. The elastic behavior of a medium with randomly oriented cracks has been exhaustively studied in the literature (Barthélémy, 2009;Budiansky and O'Connell, 1976;Henyey and Pomphrey, 1982;Hudson, 1980;Kachanov, 1993;Kachanov et al., 1994;Kemeny and Cook, 1986;Sayers and Kachanov, 1991;Schoenberg and Sayers, 1995;Walsh, 1965). However, natural fractures are not randomly distributed within the Soultz EGS site (Valley, 2007).
Their orientations have been driven by the complex geological history of the upper Rhine graben. The effective medium model is built assuming a preferential penny-shaped crack orientation in line with the Rhenish direction. In addition, the Soultz zone is within a horst area, which is bounded by subvertical N-S trending normal faults (Genter and Traineau, 1996). As a result, we develop an anisotropic effective medium model ( Fig. 4) with two sets of penny-shaped cracks with orientations defined from two angles in the spherical coordinate system: the h dip angle measured from the vertical axis and the u azimuth angle measured from the East axis, following the mathematical convention. Referring to both outcrop observations and borehole imagery, we consider the following crack orientations (Tab. 3): two North-South crack families with 155°dip, respectively eastward and westward. Two parameters characterize crack shape: length (two times the radius c i ) and aperture (denoted w in Fig. 4).

Inverse Modeling
As shown above, the simulated velocities depend on microstructural properties or parameters (crack density, crack radius and aperture) and matrix elastic parameters (bulk and shear moduli, denoted K 0 and G 0 , respectively). We now aim to solve the inverse problem. In other words, we want to estimate the microstructural parameters from the known compressional-and shear-wave velocities. From previous studies, we assume that water injection influences only the shape of the cracks, not their number. Indeed, Horálek et al. (2010) studied a pool of 45 micro-earthquakes with magnitudes between 1.4 and 2.9, which occurred during and after the 2003 massive fluid injection in the GPK3 borehole. They pointed out a massive dominance of double-couple events indicating that shear-slip is the main mechanism of the investigated micro-earthquakes. Therefore, microseismicity is mainly due to sliding on preexisting failure planes and tensile fractures, assumed to be created during injection phases, are rare.
Inverse problems can be solved using various optimization techniques (Tarantola, 2005). Within this framework, the leading idea consists in minimizing a given objective or cost function in order to find the unknown parameters. In the case under consideration, the objective function (J) is defined as: This two-term function quantifies the data mismatch in a least-square sense. The first term measures the misfit between the squared compressional-wave velocity obtained from tomography data (V tomo P ) and the corresponding numerical responses (V model P x ð Þ). The second term is very similar, but focuses on shear waves. Superscripts "tomo" and "model" refer to the tomography data and the corresponding numerical modeling results, respectively. Vector x encompasses the unknown parameters. It includes four components: the mean crack radius, the mean crack aperture and the matrix bulk and shear moduli. The minimization of the J objective function provides a x vector so that the difference between the simulated responses and data is as small as possible. To alleviate some of the drawbacks related to ill-posed problems, the unknown parameters are submitted to bound constraints (Tab. 4). The ones for matrix bulk and shear Effective medium modeling: parameters of the fractured model.  moduli are selected on the basis of the 1-D reference velocity model used in the tomography workflow (Tab. 2) as explained hereafter. Crack aperture is initialized to 10 À4 (corresponding to an aspect ratio of 10 À4 with a length of 1). It can vary between 10 À5 and 10 À3 , which is a reasonable range for crack aperture in granites (Schubnel et al., , 2006. Crack length boundaries are chosen regarding on the relationship with volumetric crack density expressed in Equation (3). Indeed a unit density can be more or less considered as a threshold for mechanical failure (Guéguen et al., 1997). In order to be far from this threshold, the upper crack length boundary is fixed to 3 in order to reach a maximal crack density of 0.5 according Equation (3). The minimization problem is then solved using the (Bound Optimization BY Quadratic Approximation BOBYQA) local algorithm. This derivative-free algorithm solves bound constrained optimization problems using a trust region methods that forms quadratic models for the objective function by interpolation. The interested reader can refer to Powell (2009) for more details. This method has been chosen principally because of its speed of convergence. At this stage, an essential issue is related to scale. It is commonplace to apply the EMT at the rock sample scale, less at the field scale. At the sample scale, the EMT relates microcracks and ultrasonic velocities. At the field scale, it links meso-or macro-fracturing to log and/or seismic data. In both cases, the scaling relationship between interrogating wavelength and REV size are similar. Furthermore, EMT methods allow to derive unrelaxed elastic moduli. These models are scale independent as long as the wavelength is larger than the REV scale. The tomography analysis used as a basis in this study resulted in passive seismic data in the frequency range 1-50 Hz with a prominent frequency around 10 Hz. Assuming a mean P-wave velocity of 6 km/s (for the whole ray paths, from the surface to the foci zone), and a mean frequency of 10 Hz corresponds to a wavelength of 600 m. Further, the horizontal dimension of the tomography mesh is 100 9 100 m. As a result, we consider that the REVof our EMT model is the same size as the tomography grid. This hypothesis is relevant because the REV is smaller than the wavelength. Thus, we claim that the fracturing scale of interest in this study is of the order of meters. We aim to characterize the meso-scale fracturing, which is the same scale than UBI observations (Valley, 2007).

Calibration of the Initial Crack Density and Matrix Parameters
Surface crack density is assumed to be constant within the REV. Its value could be derived from borehole data. Unfortunately, there is no UBI data available in well GPK2 below 3 860 m depth. Instead, we refer to the UBI data collected in deeper wells GPK3 and GPK4. Dezayes et al. (2000) estimated the cumulative fracture density in these two wells from 1 500 to 5 200 m depth. In the resolved tomography area, we inferred a 1-D crack density, often named P 10 , which is the number of fractures divided by length of scan-line. Thus, a constant P 10 crack density equal to 0.5 fractures.m À1 . P 10 is independent of the fracture length. In our modeling we need the surface crack density P 32 which corresponds to the area of fracture divided by the volume of rock mass. Even if we know that it is not obvious to link 3D crack density and 1D parameter like P10, especially in such anisotropic medium, we make the assumption to have a linear relationship between P 10 and P 32 and thus we put into the EMT model a constant P 32 of 1 m À1 . Another issue is the estimation of the crack-free matrix bulk and shear moduli, K 0 and G 0 . These two parameters are considered as unknown and are adjusted when minimizing the objective function. However, the inverse problem being ill-posed, the variations in these parameters have to be bounded to ensure regularization. To do so, we refer to the 1-D reference P-wave velocity model inferred from log data and calibration explosions within the framework of the tomography survey (Tab. 2). Then, the Vp/Vs ratio is assumed to be constant and set to 1.75. This provides the shear-wave velocities. The derivation of the bulk and shear moduli from these compressional-and shear-wave velocities makes it possible to estimate reasonable upper and lower bounds for these parameters. Thus, K 0 is considered to fall in the range of 47 to 53 GPa and G 0 of 25 to 35 GPa.
The primary calibration phase being completed, we move to the inversion process itself. The inversion is independently performed for every grid blocks populating the resolved sub-domain within the tomography grid. It entails the minimization of the J objective function by adjusting the unknown parameters: the matrix elastic moduli and above all, the crack lengths and apertures. Inversion results were not smoothed from one mesh to others. A sensitivity study performed to assess the reliability of our optimization results is presented in Appendix B.

SHEAR-WAVE SPLITTING METHODOLOGY
One of the main objectives of this paper is the comparison of the previous effective medium-based results with geophysical data such as the shear wave splitting.
Cracked isotropic media may exhibit anisotropic properties with cylindrical symmetry (Nur, 1971;Anderson et al., 1974). A lot of attention has been paid to such media, which are of great interest in seismic exploration and reservoir production (Al-Harrasi et al., 2011). The interested reader can refer to Crampin (1984) for more details.
Shear wave splitting is typically used as a tool for detecting anisotropy, in particular at reservoir scale (Verdon et al., 2009;Wuestefeld et al., 2011). When a polarized shear wave enters an anisotropic homogeneous medium, it splits into two orthogonally polarized shear waves with distinct velocities. There are two splitting parameters of interest: dt, the travel-time difference between the fast and slow waves and the orientation of the fast shear-wave polarization that is defined by the w angle in the plane perpendicular to the ray path. These two parameters make it possible to estimate the orientation and intensity of the anisotropy in the medium. Splitting has been attributed to the effects of extensive-dilatancy anisotropy, the distributions of stress-aligned fluid filled microcracks and preferentially oriented pore-space (Crampin and Booth, 1989). The pre-existing vertical cracks parallel to the direction of the maximum compressive horizontal stress (Sh max ) can open during injection while the others cracks remain closed. The polarization of the vertically fast shear-wave particle motion direction is parallel to the direction of the preferred orientation of the open cracks. In other words, the fast wave polarization is parallel to the direction of the maximum horizontal compression.
High-frequency shear waves split completely into two separate shear waves, which ease the determination of the splitting parameters. They are derived from two linear orthogonal components of the particle motion. However, the splitting parameters are even more accurately estimated from an elliptical particle motion due to the interference of fast and slow shear-waves. Parameter identification was performed using a version of the SPLIT shear code (Vecsey et al., 2008) that has been modified to handle the Soultz reservoir scale.
Shear-wave splitting is routinely observed in Soultz at the 6 surface 3-component stations. However, the records for only 2 of them, CARO and KIRK (Fig. 2), could be used to analyze the full particle motion ellipse as they were not affected by any disturbance due to converted arrivals.
The data processing was manually performed following in two steps (Fig. 5): signal pre-processing: we selected events with magnitude greater than 1 and high signal to noise ratio. Then, we chose a 0.5 s length window containing shear-wave and re-sampled the signal to 1 000 samples/s using the Wiggins (1976) interpolation method; splitting evaluation: we input station and event coordinates to rotate the signal to ray coordinate system (Sileny and Plomerova, 1996), set interactively the window length from the beginning of the S signal until the particle motion was properly constructed. We applied an Eigenvalue method as recommended by Vecsey et al. (2008) for direct shear-waves, and calculated error estimates using a bootstrap method. Finally, we only considered the results with dt and w errors less than 0.15 9 10 À2 s and 10°, respectively.

Crack Features Maps
Crack properties evolve with time because of hydraulic stimulation. Figures 6 and 7 display the results obtained for inversion of crack features (aperture and length). Figure 6 corresponds to North-South vertical cross-sections including well GPK2. Both crack lengths and crack apertures are presented in Figure 6. Besides, Figure 7 shows horizontal slices with only the evolution of it normalized, the reference state being 1. Therefore, values larger than 1 indicate crack extension while values lower than 1 are associated with crack shortening. While it is easy to interpret crack extension from a physical point of view, it could be difficult to understand the meaning of crack shortening. In this case, we have to consider that cracks are not stimulated on all their lengths. We may envision that a pore pressure drop within the crack network can contribute to close part of pre-existing cracks. Figure 6 shows that crack aperture is not an influential parameter during stimulation test: it remains more or less the same as the initial value of 10 À4 as shown by large green areas for all sets. This result is in agreement with theoretical studies (Simpson et al., 2001;de Dreuzy et al., 2002). On the contrary, crack length exhibits significant variations during hydraulic stimulation. Right after the beginning of injection, for the time periods identified by sets 1 and 2, large blue areas are detected on the crack length maps: they correspond to high crack length values. When the flow rate is set to 40 L/s, the crack length first decreases, especially close to the openhole section of GPK2 (set 3), before increasing again but in a moderate and diffuse way (sets 4 and 5). Again, the following flow rate increase to 50 L/s contributes to reduce crack length in the vicinity of GPK2 open-hole section (white cloud visible on set 6), similar to what is observed for set 3. For the following sets of the C phase, (sets 7 to 12), we observe higher values in the vicinity of the open-hole section (below 4.6 km depth). This general picture does not change within this period of time whatever the depth is. When hydraulic stimulation is over (sets 13 and 14), crack length tends to reach the reference state, even if there is still extension in the vicinity of the GPK2 well.
The horizontal cross-sections reported in Figure 7 evidence lateral variations in crack lengths at 3 depths: 4.2 (close-hole section), 4.6 (open-hole section) and 5.0 (below GPK2 borehole). It shows that a disturbed area, i.e. an area with larger cracks, is centered close to the injection well. The injection phases lead to an enlargement of this area. After shut-in, the area reduces once again around GPK2 (sets 13 and 14).
The relevance of our inversion results can be strengthened by the analysis of the J objective function evolution (Tab. 5). Whatever the target set, the final normalized objective Steps for the shear-wave splitting processing at stations KIRC and CARO. The example shown is an event that occurred on July 2, 2000, at 00h 39min UTC. a) the three components of the seismic signal (Vertical, North-South, East-West), b) magnification of the shear-waves on the horizontal components, c) rotation of the seismic signal to the ray coordinate LQT system (Llongitudinal, Qradial, normal to [L-T (transverse)] plane) and selection of the time window for splitting evaluation, d) left: construction of the elliptically polarized particle motion in the (L, Q, T) coordinates and rotation to the F (fast) and S (slow) axes, where the particle motion should be linearly polarized; right: diagram of the misfit function obtained by Eigenvalue method. The white circle in the blue minimum domain shows the splitting parameters: time delay dt of the slow split shear-wave and rotation angle W of the fast wave polarization in the (Q, T) plane. The fast polarization direction W is defined by two Euler angles, azimuth / and inclination h, e) histograms for splitting parameters W and dt. See Vecsey et al. (2008), for a complete description.
function is very small (mean value around 10 À5 ). In addition, we check that the objective function was decreasing monotonically and we observed a single valley with a single minimum, which makes the minimization process easier.

Extension and Opening Rates
In this study, the two crack parameters are the length and the aperture. We now calculate extension and opening rates. The e ext extension rate and the e op opening rate are defined as: L i is the optimized length, L ref the reference crack length (equal to 0.5 and corresponding to the smallest mean value obtained for set 3 between 4 and 5 km depth), x i the optimized aperture and x ref the reference crack aperture (equal to 5 9 10 À5 and corresponding to closed cracks). An extension rate of 100% corresponds to a crack length of 1 and reciprocally an opening rate of 100% corresponds to a crack aperture of 10 À4 . The reason why we focus on global extension and opening rates is that we want to compare a tomography set to another one. We first determine the rates for all depths investigated by tomography (Fig. 8a). Second, we focus on given levels in an attempt to capture the influence of the well (Fig. 8b, c): at the level of the open-hole section (4.5-4.9 km depth), above the open-hole section (4.0-4.4 km depth) and below (5.0 km depth).
As already suggested by Figure 6, no significant variation in the opening rate is evidenced when considering the whole depth range. However, if studying more specifically the area above close-hole, it turns out that sets 1 and 2 are associated with fracture opening (Fig. 8b). Then, after set 3, all the curves converge and the rate gets close to 100%. This result may be puzzling if considering only crack aperture. However, as suggested by lab experiments, elastic waves are very sensitive to variations in crack aspect ratio, not crack aperture only. Precisely, the extension rate shows even more variations with time. The same general trend is observed for the three depth sections with some differences for sets 2, 5 and 7 (Fig. 8c). Except the close-hole section, the extension rate is minimal for set 3. Again, this was already emphasized from the analysis of the crack length maps (Fig. 6, 7). A second minimal value is achieved for set 6. Remember that sets 3 and 6 correspond to the time period right after a sharp flow rate increase. Other maximal values are obtained with sets 4 and 9 (Fig. 8a). These sets actually correspond to the middle of the B and C injection phases. The extension rate is higher in the area around the open-hole section of GPK2 (violet/ purple curve in Fig. 8c) than in other areas, only except Horizontal cross sections at selected depths with the inverted crack lengths. Black dots are the projections of events used to obtain the velocity models. for set 3.This is in agreement with the results shown on the vertical cross-sections (Fig. 6). Maximum values for extension rates are reached for set 4. Second maximum values are reached for close-hole section during set 9 and for the open-hole section during set 12.

Shear-Wave Splitting
We ended up with about 15 events with high quality splitting evaluation for each tomography subset. We calculated w mean and an index of anisotropy Ia = (Rdt i /t si ) / N. dt i and t si are the significant shear wave splitting value and the shear wave splitting time delay along the ray of the fast shear-wave for the i th event. N is the number of events. We also checked that the centroids of the hypocenters of the selected events are very close (less than 100 m difference) for all subsets. Figure 5 explains the shear-wave processing used in this study. The anisotropy evidenced on the surface stations can result from various sources like the sedimentary layers above the granitic basement or the stimulation process performed a few years earlier in the granite. As our analysis is based only upon the relative variations in the splitting parameters at the two stations during the 2 000 stimulation, our results are independent of the aforementioned sources and they can be used to quantify the influence of stimulation on crack features. First, we observe two main features: the direction of the fast shear-wave polarization does not significantly change during stimulation and, its mean value is 155° (Fig. 9).
This direction is roughly parallel to the maximum horizontal stress orientation estimated to be 158°by Dorbath et al. (2010) from the inversion of well-constrained focal mechanisms.
The two Ia curves (Fig. 10) present similar behaviors: anisotropy decreases when the injection rate increases. This is particularly clear for set 3, when the water rate is increased from 30 to 40 L/s. The same behavior is also observed for set 6, but not as clearly, when the rate is set to 50 L/s. For these two sets, the seismic velocities came back to their initial values (Calò et al., 2011). Another common feature of the two curves is the continuous increase in Ia from set 10 till the end of seismic monitoring. This result could be related to the spreading of seismicity in an area getting wider with time at the end of stimulation. The increase is even reinforced after shut-in (Cuenot et al., 2008).
From a mechanical point of view, two trends can be distinguished from Ia curves: (1) between set 1 and sets 5/6, the main trend is a slight decrease with many variations and (2) beyond sets 5/6, there is a constant increase. Inspiring from results of Guéguen and Sarout (2011) concerning EMT in the case of vertical cracks, one could suggest that trend (1) corresponds to a re-saturation regime, the injected water replacing the fluid in place, perturbing the saturation state of fractures and also probably the crack features but not permanently; and trend (2) corresponds more to a "damage" regime, with permanent changes of crack features, implying a constant Ia increase.
Last, the differences in the absolute Ia values (Tab. 6) are easily explained by the locations of the seismic stations relatively to the reservoir. CARO is located just above the open hole section of GPK2, while KIRC is located to the northwest. Therefore, the rays travelling to KIRC cross a larger stimulated area.

DISCUSSION
If we focus on set 3 in Figures 5 and 6, an interesting local trend is pointed out. Crack lengths are larger around the close-hole section of GPK2 than in the vicinity of the open-hole section. However, microseismic events do not preferentially migrate upwards (Calò et al., 2011), usually meaning that there is no preferential flow direction. The EMT results could be explained by considering what is going on right after a sharp increase in the flow rate. During a very short period of time, a transient regime takes place. As described in Weidler et al. (2002), the reservoir is split into two sub-domains with distinct flow behaviors: a shallow part with a well-connected system where cracks are open on all their entire lengths and a deep part behaving like a closed system. However, effective medium modeling assumes a pure elastic behavior for the material, which is not suitable for describing a transient regime in a closed system. This is why our inversion process entails crack shortening in the small area close to the open-hole section (Fig. 8c). Shear-wave splitting is an useful method for characterizing anisotropy from travel-time differences. This makes it possible to estimate anisotropy indices. In this section, we compare these indices with the global extension rates introduced previously (Fig. 8a, 10, Tab. 6). These extension rates are the ones computed at all depths sampled by the velocity models, i.e., between 4.0 and 5.0 km depth.
The agreement is good for the first 9 sets, which correspond to about the first 4 days of the injection test. The lowest Ia values (set 3, Fig. 10) correspond to negative extension rate values (Fig. 8, crack shortening in fact). The maximum extension rate calculated for set 4 (Fig. 8a) is associated with a relatively high Ia value (Fig. 10). However, the comparison  Evolution of the index of anisotropy of the S waves for two seismic stations during GPK2 stimulation.
is poor for the final data sets. The index of anisotropy is then shown to depart from the crack extension rate. In short, it tends to increase for sets 10 to 14 while the extension rate ratio remains more or less the same or even slightly decreases. Shear wave splitting is mostly due to vertical cracks parallel to the direction of the maximum compressive horizontal stress. An increasing index of anisotropy (Fig. 10) is associated with a stronger splitting, hence to the extension of preferential cracks. Wuestefeld et al. (2011) also noted a temporal increase in splitting delays for steep arrivals during the stimulation of a gas reservoir. They explain this result as a "halo" of increased fracture density around the stimulated zone. Crack density is assumed to be constant within the EMT model, but it is clear that an increase in the number or the length of favorably oriented cracks will lead to the same anisotropy effect.
Two hypotheses of our EMT model are: cracks are oriented according to the regional geological direction; the variation in the extension rate is a macro parameter derived at the REV scale, i.e. an average for all crack populations embedded within the effective medium.
The influence of two factors could be strengthen due to hydraulic stimulation. During the first stimulation stages, the crack network is subjected to a huge fluid pressure change; injection flow rates impact crack properties, hence shear wave splitting. In such conditions, the system is said under hydraulic control. During the last stages of the stimulation process, the flow rate remains constant. This means that we potentially move from a hydraulic control to a geological control through regional stress orientation. In other words, the horizontal field stress could vary locally due to the hydraulic stimulation: this would explain the discrepancy observed between the extension rate and the time-shift between slow and fast shear arrivals.
An important issue to be discussed at this point is the assumption according to which surface crack density P 32 is constant during the whole injection experiment. The underlying question is about the capability of hydraulic stimulation to generate new cracks. Cornet et al. (2007) showed that the Soultz granite was in a critically stressed state from a mechanical point of view. Pore pressure variations of only 2 to 3 MPa induce microseismic activity and values hardly larger produce macroscopic failure. However, there are cracks in the granite both at the micro and meso-scales. As stated above, our EMT-based inversion method strictly focuses on the meso-scale. Thus, we only pay attention to macroscopic failure localization. The processes inducing microcracking such as thermal cooling are beyond the scope of this paper. This approach is consistent with previous studies (Zhao and Mizuno, 1999;Mishra and Zhao, 2003;Adelinet et al., 2011). The next step consists in identifying the mechanisms that can explain our results keeping in mind that crack density is constant. During sets 2 and 4, wellhead pressure slightly decreases after the sharp increases related to flow rate jumps (Fig. 1). This suggests an increase in the available storage volume. The injection flow rate may be not significant enough to fill up this additional space. During the same periods of time, our modeling shows higher crack extension rates ( Fig. 6-8). The induced dilation of cracks due to their extension could explain such a pressure reduction at the wellhead. This kind of process has been already observed during a hydraulic stimulation test in the Cornubian granite in England (Batchelor, 1982;Crampin and Booth, 1989). Finally the effective medium model built in this study seems to suit well to the very specific problem of geothermal stimulation of a fractured reservoir.

CONCLUSION
This paper introduces a novel approach to explain the behavior of the Soultz reservoir, which was submitted to a massive hydraulic stimulation. We solved an inverse problem using an effective medium model that makes the link between the available tomography data and the unknown parameters describing the state of fractures. The results obtained stressed that a strong increase in the injection rate leads to a decrease in crack length without influence on crack aperture. This provides evidence in favor of crack shortening. The equilibrium is achieved after some time with injection performed at the same flow rate. In such conditions, crack length keeps from fluctuating.
Extension rates deduced from the effective medium based inversion process correlate with time shifts recorded between the slow and fast shear-wave arrivals. The good agreement between the two independent data sets suggest that the processes occurring at the meso-scale (extension/ shortening of cracks) are quite similar to those occurring at larger scales even though the corresponding crack networks are different. In addition, in this paper, we distinguished clearly the effect of crack extension from the effect of crack closure. This helps better understand the effect of stimulation.
Finally, effective medium modeling is a very useful tool to explain reservoir and geophysical data in terms of microstructural evolutions. This is all the more essential than the maximization of productivity depends on reservoir structure, whatever the nature of the target resource (thermal energy or fossil resources). The Soultz reservoir is an ideal study case because of its homogeneous petrophysical properties and the large number of data collected at different scales. To date, an important issue is to extend the workflow proposed in this paper for the Soultz reservoir to more complicated field cases.

APPENDIX A
Evolution of the shear-wave seismic velocity at 4.6 km depth (Fig. A1) and corresponding vertical cross-sections (Fig. A2) during the 2000 stimulation in GPK2. Figures are chronologically ordered from set 1 to 14. The same results for compressional wave velocity were published in Calò et al. (2011). Evolution of the shear-wave seismic velocity at 4.6 km depth during the 2000 stimulation test. Images are in chronological order from set 1 to set 14. The black dots are the projections of the events used to obtain the shear-wave velocity models.

APPENDIX B
We repeated our inversion process several times to assess the reliability of our results and to evaluate the sensitivity of the numerical responses obtained with respect to the effective medium approach. The data used for this purpose correspond to the 5 th set and the 4.6 km depth. Vertical sections of the shear-wave velocity models along traces A-B reported in Figure A1. The black dots are the projections of the events used to obtain the shear-wave velocity models. Blue lines are the projections of GPK2. The thick part of the borehole trajectory corresponds to the open-hole section of GPK2 well.
In a first step, we investigated the influence of crack orientation and direction of wave propagation on the inversion of crack length (Fig. B1). We considered as the true model the one obtained with parameters described in this paper (Fig. B1a): two NS crack families with dip 65°and direction of wave propagation in agreement with tomography data (maximum deviation angle of 40°with respect to the vertical axis). Then, two crack distributions are considered (Fig. B1b): an isotropic one and another one for NS cracks with a dip 45°. The inverted crack length inferred from the isotropic distribution is quite close from the true model. Conversely, the single family of cracks with 45°dip yields a crack length distribution quite different. Finally, two wave propagation directions are tested (Fig. B1c): a vertical and an oblique one (45°). The inverted crack length obtained when considering a vertical propagation is quite similar to the true model although the oblique propagation leads to very different results. This sensitivity study emphasizes the fact that it is essential to properly characterize the true geological medium in order to input realistic parameters in the effective medium modeling.
In a second step, we studied the sensitivity and the robustness of the optimization process (Fig. B2). We considered as the true model the compressional wave velocities derived from the tomography study (Calò et al., 2011) displayed in Figure B2a. Then, we carried out various tests to evaluate the influence of the number of optimization parameters: optimization with one parameter (crack radius) in Figure B2b, 2 parameters (crack radius and aperture) in Figure B2c and 4 parameters (crack radius, crack aperture, bulk and shear matrix moduli) in Figure B2d. The results obtained showed that the number of optimization parameters is not a crucial parameter for the achievement of the optimization process. Whatever this number, the general trend is preserved, even if extreme compressional wave velocity values are not reached. Sensitivity study on direct parameters using data extracted from set 5 at 4.6 km depth. Comparison is made on the inverted crack length. Black dots are the projections of events used to obtain velocity models. a) True model: geological distribution of cracks and direction of propagation in agreement with tomography input data. b) Effect of the crack orientation on the inverted length, isotropic distribution (left picture) and cracks with dip 45°(right picture). c) Effect of the direction of the wave propagation: purely vertical (left picture) and oblique propagation (right picture).

Figure B2
Sensitivity study on inversion parameters using data extracted from set 5 at 4.6 km depth. Comparison is made on the compression-wave velocity (expressed in km/s). Different inversion parameters are tested: crack main radius (R1), crack aperture (R3) and bulk and shear moduli of the matrix (matrix). a) The true model is the compression-wave seismic velocity model calculated from the tomography study (Calò et al., 2011); b) compression-wave velocity model inverted with only R1 as optimisation parameter; c) P-wave velocity model inverted with R1 and R3 as optimisation parameters; d) compression-wave velocity model inverted with R1, R3 and matrix moduli as optimisation parameters.