Magnetotelluric study of the Remiremont-Epinal-Rambervillers zone of migrating seismicity, Vosges (France)

– The magnetotelluric method has been used to image the deep electrical structure of the Remiremont-Epinal-Rambervillers region in the French Vosges Massif, which has presented a significant seismic activity in the past decades. Several earthquakes of moderate magnitude (up to 5.1) occurred in this area with a systematic migration along a nearly N-S direction. Inversion of the magnetotelluric data reveals zones of high electrical conductivity. A large conductive body presents a significative spatial correlation with the region that was most recently affected by earthquakes. This conductive body is interpreted as a consequence of the presence of a fluid filled basement fault network in proximity to the zone affected by the last seismic crisis, where fluid pressure diffusion takes place for several years after the main shock and participates in maintaining a microseismic activity. Etude magnétotellurique des zones de migration de sismicité de Remiremont-Epinal-Rambervillers, Vosges (France) Mots-clés. – Magnétotellurique, Fluides, Réseau de failles, Sismicité, Massif des Vosges. Résumé. – La méthode magnétotellurique a été utilisée pour imager la structure électrique profonde de la région de Remiremont-Epinal-Rambervillers dans les Vosges, où a été enregistrée une importante activité sismique pendant les trois dernières décennies. Plusieurs séismes de magnitude modérée (jusqu’à 5.1) ont eu lieu dans cette région avec une migration systématique des événements sismiques le long d’un axe orienté environ Nord-Sud. L’inversion des données magnétotelluriques révèle plusieurs zones de conductivité électrique élevée. Un corps électriquement conducteur de grande extension présente une corrélation spatiale significative avec la région affectée par les tremblements de terre les plus récents. Ce corps conducteur est interprété comme la signature électrique de la présence d’un réseau de fractures saturé de fluides conducteurs, à proximité de la zone affectée par cette dernière crise sismique, et où la diffusion de la pression de fluide prend place pendant plusieurs années après le choc principal et participe ainsi au maintien d’une activité microsismique.


INTRODUCTION
Fluid flow and fluid pressure appear to play an important role in the dynamics of fault zones. The presence of high-pressure fluids at depth could cause fault failure during the seismic cycle [Byerlee, 1993;Rice, 1992]. Also it could play an important role in fault propagation and trigger shallow earthquakes [Mekkawi et al., 2003]. The increase in pore pressure can drive a fault system to failure. Stein [1999] suggested that stress changes as low as 0.1 bar are sufficient to induce failure in critically stressed zones. High-pressure fluid migration has also been proposed to explain the migration of aftershocks [Miller et al., 2004], or the overpressure found in accretionary wedge décollement several years after a large earthquake [Bourlange and Henry, 2007].
A series of recent magnetotelluric studies at the San Andreas fault (SAF) and other active areas imaged electrical conductivity anomalies in the crust which are related to earthquakes distribution [Becken et al., 2011;Yamaguchi et al., 2010;Wannamaker et al., 2010;Yoshimura et al., 2009;Brasse et al., 2009;Brasse and Eydam, 2008;Jiracek et al., 2007;Mekkawi and Saleh, 2007]. They have suggested that these conductive anomalies could result from the presence of crustal fluids and/or ore minerals (e.g. graphite, salt, sulphides) accumulated in the fractures system. The San Andreas fault has revealed a pattern of migrating earthquakes explained by fluid migration near Parkfield [Johnson and McEvilly, 1995]. In this region, the magnetotelluric technique has been used to image subsurface electrical resistivity. It has revealed a low resistivity zone extending to a depth of 2-3 km of the SAF, interpreted as a zone of fractured rock filled with fluids [Becken and Ritter, 2012;Becken et al., 2008;Unsworth et al., 2004Unsworth et al., , 2000Unsworth et al., , 1997Bedrosian et al., 2004].
The Remiremont-Epinal-Rambervillers region in the French Vosges Massif is a zone of significant seismic activity [Audin et al., 2002]. A systematic migration of the seismic events has been observed along a nearly N-S direction at a rate of about 5-10 km/yr. The migration pattern has been interpreted as the result of the propagation of transient pore pressure changes along a zone of low permeability. We present here the results of a magnetotelluric survey that was conducted in the Remiremont-Epinal-Rambervillers area in order to test the presence of fluids at depth in the region.

Geology and tectonics of the Vosges massif and Rhine graben
The Remiremont-Epinal-Rambervillers region is located on the eastern edge of the Paris sedimentary basin close to the unconformity on the Vosges massif ( fig. 1). Further east, NNE striking normal faults separate the Vosges massif of the Rhine graben. In the Hercynian Vosges massif, the Saxo-Thuringian Vosges to the north is separated from the Moldanubian Vosges to the south by the WSW-ENE Vittel fault zone. Another important tectonic feature is the NNE-SSW Sainte-Marie-aux-Mines fault [Dercourt, 2002]. These two faults are considered inactive [Fluck et al., 1991].
The Rhine graben is a part of the Cenozoic rift system of Central Europe. This graben separates the Vosges massif from the Black Forest massif. The Oligocene rifting of the Upper Rhine graben, which was associated with an E-W extension, terminated about 20 Ma ago at the end of the Aquitanian [Brun et al., 1992;Ziegler, 1992;Edel et al., 2006]. The Alpine convergence between Africa and Europe led to the rotation of the main stress direction to NW-SE in the stable western Europe at the end of the Miocene [Illies, 1974[Illies, , 1978. This direction is determined by focal mechanisms showing that the Rhine Graben is dominated by left-lateral shear on N-S to NNE-SSW faults [Abhorner, 1975;Bonjer et al., 1984;Delouis et al., 1993;Plenefisch and Bonjer, 1997;Haessler and Hoang-Trong, 1985]. It is also documented by in-situ stress measurements [Illies and Greiner, 1979] and displacements measured along recent faults [Villemin and Bergerat, 1987;Larroque and Laurent, 1988;Lacombe et al., 1993]. The world stress map [Muller et al., 2000] indicates a roughly NW maximum horizontal stress ( H max ) and SE minimum horizontal stress ( H min ).

Remiremont-Epinal-Rambervillers seismicity
Historical seismicity is poorly documented in this region, except for the earthquake that occurred in 1682 near the city of Remiremont [SIRENE, 1998]. This event is one of the most damaging earthquakes known to have occurred in France (Io = VIII, MSK). Since 1980, several earthquakes have been recorded in the region by the National Seismological Network RéNaSS ( fig. 2). The local magnitudes (Ml) of these earthquakes range between 1 and 5.4. Two moderate seismic events occurred in the region: a Ml = 4.8 earthquake occurred 6 km NW of Remiremont in December 1984, and a Ml = 5.4 magnitude earthquake occurred 4 km SE of Rambervillers in February 2003[Cara et al., 2005. Both earthquakes were followed by several smaller magnitude earthquakes. A previous seismic activity occurred in 1973-1974 to the east of the city of Epinal [Audin et al., 2002]. These last authors noticed that the epicenters of the Epinal and Remiremont seismic events were located on a NNE trending strip extending from the epicentral area of the 1682 historical event. This directional trend is confirmed by the 2003 Rambervillers seismic events that occurred further north. The 1984 main earthquake focal mechanism suggests a N-S fault plane with a sinistral lateral slip. The focal mechanisms of the two major shocks of the 1973-1974 Epinal crisis are not well constrained. However, previous authors supposed focal mechanisms similar to that of the 1984 event [Haessler and Hoang-Trong, 1985]. The 2003 Rambervillers earthquake focal mechanism indicates a normal faulting along a NW-SE striking fault dipping 40 o NE [Cara et al., 2005;Jacques et al., 2003]. These focal mechanisms are consistent with a NW maximum horizontal stress, perpendicular to the Alpine arc.

PHYSICS OF THE MAGNETOTELLURIC METHOD
The magnetotelluric method (MT) is an electromagnetic geophysical method which uses measurements of the natural magnetic (H) and electric (E) fields at the surface of the earth to determine the distribution of electric resistivity within the Earth. The invention of the MT technique is credited jointly to Tikhonov [1950] and Cagniard [1953]. Cagniard was the first to actually perform scalar field measurements. The tensor MT method was proposed by Cantwell [1960]. This approach defines and measures a tensor relationship between the electric and magnetic fields at the earth's surface.
Overview of the basic aspects and the application of the magnetotelluric methods in geophysical studies are described by several authors [Vozoff, 1972[Vozoff, , 1986[Vozoff, , 1991Kaufman and Keller, 1981;Yungul, 1996;Sharma, 1997]. It is usually assumed that the medium is linear, isotropic and homogeneous. At a given frequency two orthogonal measurements of magnetic fields (H x ) and (H y ), are related to electrical fields (E x ) and (E y ) on the same axes through a tensor impedance defined by: The complex elements of the impedance tensor Z in equations (1) and (2)  (1) is the ratio of orthogonal components of the electric and magnetic fields at a given frequency, and T = 2p w / is the period of the wave. The geomagnetic deep sounding method (GDS) measures the horizontal and vertical components of the magnetic fields. The ratio between the horizontal and vertical components is most sensitive to the lateral variations in conductivity [Parkinson, 1959;Schmucker, 1970]. The induced vertical component (Hz) field is related to the inducing field by: where H x and H y are horizontal magnetic fields in the north (x) and east (y) directions, A and B are complex functions of frequency. The real parts of these functions are of great interest. They allow calculating induction arrows (vectors), which point towards more resistive areas of the substratum. Their direction is therefore perpendicular to the geologic strike. Length and angle of induction arrows are expressed as T(Re, Im) and q(Re, Im): Two-dimensional MT fields can split into two independent modes TE (transverse electric) and TM (transverse magnetic). In the first mode, the electric currents flow in the geologic strike direction whereas in the TM mode, the currents cross the geologic strike at right angle. The apparent resistivity and phase computed for each mode can vary considerably with site location [Vozoff, 1986].
The apparent resistivities determined by MT method can be affected by frequency-independent offset, called static shift [Sternberg et al., 1998]. It occurs when electrodes are implanted on or near anomalous body (3D structures, near surface inhomogeneities, variations in geology, fluid content for example). If the static shift is too severe, corrections have to be performed on one or both TE and TM modes. Active induction technique (transient electromagnetic sounding and audio-magnetotelluric for example) can be used to correct the static shift [Schnegg, 1998].

MAGNETOTELLURIC DATA ACQUISITION
A first magnetotelluric (MT) study was conducted in order to determine the location of the active faults associated with the 1984-85 seismic events [Mekkawi, 2003]. A WNW-ESE magnetotelluric profile (C-D on figure 2) was carried out across faults documented on geological maps [Durand et al., 1988;Hameurt et al., 1985;Hameurt and Vincent, 1979;Minoux, 1978;Chrétien et al., 1974;Desprez et al., 1971]. The results of this previous study suggested the presence of a broad conducting body to the NE of the profile. This was evidenced for example by the real induction arrows for the period 12.4 s that point at all sites toward the SW (fig. 3).
The MT profile of the present study (A-B on figure 2) was carried out on a direction perpendicular to this previous profile (NE-SW) in order to image the region affected by the seismic cluster events, also corresponding to the direction where the broad conductive body to the NE was expected. The direction chosen for the profile approximately follows the seismicity pattern observed in the Remiremont-Epinal-Rambervillers seismic zone. The data collection was achieved during two measurement campaigns. The profile length is about 41.5 km. Nineteen observation stations have been realized, the distance between two stations varies between 1.0 and 4.0 km ( fig. 2). Two acquisition stations constructed at the University of Neuchâtel were used. ECA CM16 coils were used to measure short period (0.0046-10 s) magnetic field variations (AMT) [Andrieux et al., 1974], whereas longer periods (MT) were measured with ECA CM11E coils (1-560 s). Electric field variations were measured using non-polarisable Ag-AgCl electrodes immersed in saturated KCl solution [Filloux, 1987]. Telluric lines were deployed N-S and E-W over 50 m. Magnetotelluric data were logged continuously during about 40 hours at each location.

Data treatment and impedance tensor analysis
The MT data consists of five components (E x , E y , H x , H y and H z ). Due to water infiltration in one vertical coil, the vertical magnetic component (H z ) could not be measured for all locations. The MT response function estimation was carried out with the single station least-square scheme using data dispersion as a selection criterion [Schnegg, 1990] to determine the apparent resistivity and phase for each relative surface position and frequency.
The induction arrows are shown for 11 stations ( fig. 4). For periods below 1 s, the orientation can change from one station to the next. This probably reflects important heterogeneities in the resistivity structure of the shallow underground. For periods higher than 10 s, the real induction arrows are more homogeneously orientated SW, yielding a regional strike close to N030 o E. This orientation is similar to the strike of the earthquakes distribution ( fig. 2 In a first approach, we propose to use the simplifying assumption of a two dimensional structure. The apparent resistivity and phase for Z xy and Z yx (TM and TE) impedance tensor are shown in the figures 5a and 5b. An increase of the apparent resistivity can be seen at periods 0.0046-10 s, with accompanying phase drop. This may reflect the presence of water-saturated sediments above a cristalline substratum, especially on the northern part of the profile.

Depth inversion of the data and model description
The inversion of the data was done with the 2D magnetotelluric inversion program REBOCC . The data set consists of TE and TM apparent resistivity and phase responses at 41 periods from 0.0046 to 560 s. The initial model was chosen as a 10 ⍀m half-space with 31 layers (plus 10 air layers for TE mode). Several sensitivity tests were carried out. We present the result of the full TM and TE inversion. The static-shift distortion was automatically adjusted in the inversion process. A subset selection of the data was used (one data every 3 periods) which is usually enough to get a reliable inversion .
The model presents several zones of distinct apparent resistivities ( fig. 6). The first 3 km present a moderate to low resistivity. The electrical resistivity of this top layer decreases globally toward the north. This feature is interpreted as reflecting the sedimentary layers of the Paris basin covering the Hercynian basement. In this region (from top to bottom) the basin is composed of 150 m of Triassic sandstones and conglomerates (Buntsandstein),~100 m of limestones and mudstones (Muschelkalk) and~100 m of mudstones (Keuper). The cover thickness never exceeds 500 m. The Buntsandstein formation constitutes a good aquifer, which is a possible explanation to the low resistivity values of the near surface layer especially in the northern part.
Below this near-surface moderately conductive layer, the model reveals the presence of four bodies with distinct apparent resistivity signatures (a, b, c, and d). From south to north, the upper crust presents a high resistivity (~5000 ⍀m), at least 10 km extent (a: from station #1 to #5), a~5 km extent zone with a moderate resistivity of about 500 ⍀m (b: from stations #12 to #14), a~5 km extent zone (c) of high resistivity (higher than 1000 ⍀m) between stations #14 and #16, and finally a broad conductive zone (d) with a resistivity lower than 100 ⍀m (with a minimum of 10 ⍀m) of about 10 to 15 km extent.

INTERPRETATION AND RELATION TO SEISMICITY
An episod of seismicity was observed in the Epinal area during 1973-1974[Audin et al., 2002. The earthquakes that have affected the region between 1980 and 2005 have been projected on the magnetotelluric derived electrical resistivity profile ( fig. 6). All the earthquakes that occurred on a band 0.4 o longitude large (~30 km large) centered on the profile are plotted. The uncertainty on earthquakes depth determination is in the order of 5 km. All the earthquakes but one occurred in the first 10 kilometers, corresponding to the expected depth of the brittle layer in this normal continental crust. The 1984 Remiremont earthquake (M 4.9) occurred in the transition region between the high resistive body to the south (a) and the smaller moderately resistive body (b). The more recent 2003 Rambervillers earthquake (M 5.4) occurred in the vicinity of the broad conductive body (d) on the north. The magnetotelluric image shows that in 2006 a large conductive body was present south of the 2003 Rambervillers earthquake hypocenter. We think that this spatial relation between the seismic swarm and the electrical conductive body results from the presence in the underground of a network of fractures filled with an electrical conductive fluid, where seismic slip occurred and which are significantly connected one to each other to yield a large regional electrical conductivity anomaly. A similar interpretation was proposed at the San Andreas fault to explain the presence of a similar low resistivity zone .
The microseismic activity observed in the Rambervillers region has continued since the main shock of February 2003 and until 2011 at least. Such a long period of activity cannot result of only mechanical adaptation that usually follows earthquakes. Fluid diffusion in the fracture network is probably also involved and responsible of this long-term seismic activity.

Fault zone porosity considerations
In the broad anomaly of high conductivity (d), the lowest values of the calculated resistivity are about 3 ⍀m at about 10 km depth. With the same magnetotelluric method, a similar resistivity of 3 ⍀m at 3 km depth was determined below the San Andrea fault [Unsworth et al., 1997] and a value of 30 ⍀m at 3 km depth was determined at the Kalabsha fault in Upper Egypt [Mekkawi et al., 2005]. Both publications estimate the value of the fracture porosity in the fault zone from the apparent resistivities obtained in the fault zone using Archie's law, with the assumption that the electrical conductivity is the consequence of interconnected pore fluids. Archie's law expresses as: where r is the effective resistivity of the rock, f is the pore fluid resistivity, is the porosity, and m is an empirical constant, varying accordingly to the shape and interconnectivity of the pores. A higher degree of pore interconnection has been characterized by a smaller Archie's law exponent. It commonly lies between 1 and 2. An intermediate Archie's law exponent of maximum 1.5 is more likely to occur in the middle crust. It is known from boreholes that highly saline fluids exist to at least mid-crustal depths. The German Continental Deep Drilling Program (KTB) sampled fluids with dissolved loads in excess of 68 gL -1 at 4000 m depth [Huenges et al., 1997]. This translates to electrical resistivity of~0.025 ⍀m at 260 o C, the temperature at a depth about 9 km at the KTB site [Simpson and Bahr, 2005]. The KTB site is situated 600 km east from Remiremont, in a geological setting similar to that of the present study area, consisting of Hercynian Moldanubian unit basement. Using KTB fluid resistivity as an estimate, the porosity in the Remiremont fault can be estimated to range between 1% (for m=1) and maximum 4% (for m=1.5). A higher porosity range (8-30%) was proposed for the San Andreas fault, where a f value of 0.26 ⍀m was based on the chemical composition of fluids sampled at a nearby well [Unsworth et al., 1997]. The 1-4% porosity we estimate for the Epinal-Rambervillers region is very similar to the 1% porosity measured at 9 km depth in the KTB hole.

Fluid diffusion in a faults network
We proposed in the section 6 that the broad conductive anomaly could be the indication of a network of connected fractures. Any increase of the fluid pressure in such a network should induce fluid migration along the pressure gradient through diffusion processes. Similar explanations have been proposed for the seismic crisis monitored in the Aigion region of the gulf of Corinth in western-central Greece in the period [Bourouis and Cornet, 2009. Seismic tomography methods applied to the 2001 seismic crisis revealed a low-angle seismogenic zone dipping 20 o toward the north hypothesized as a possible highly fractured zone with possible deep fluids circulation [Gautier et al., 2006]. The 2003The -2004 seismic crisis involve deformation processes with short (few days) and long (several months) time characteristics that could correspond respectively to mostly slip along pre-existing structures for short-lived crisis, and fluid diffusion for long term behaviour [Bourouis and Cornet, 2009]. Fluids originated from the seismogenic zone have migrated upward, in combination with regional downgoing flux of meteoric water. These authors proposed that the upward migration of fluid results from transient pressurization generated by the deformation process. Even if the seismic activity evolution is not as well constrained in the Remiremont area, similar fluid diffusion processes probably occur.

Coupling of fluid diffusion with fault slip
Diffusion of fluids in a fracture network would change the local effective stress and then possibly promote slips on preexisting fractures that were initially in a state of stress close to failure. Such a diffusion of fluids necessitates the presence of a fluid overpressure and a sufficient amount of connectivity of the fracture network. The main seismic shock could potentially have broken some sort of hydraulic barriers, whether at the pore or fracture scale, thus creating some hydraulic pathways. This transient fluid overpressure would then be able to diffuse in the fracture network, and potentially generate new fault slips. This proposed mechanism is expected to be transient and to attenuate with time. Some works have also explored the possibility of self-propagating solitary wave of fluid pressure in systems coupling elastic deformation with fracture permeability modification such as the increase of the fracture permeability with any increase of the fluid pressure in the fracture, due to its dilatancy, in the décollement zone of accretionary wedges for example [Bourlange and Henry, 2007]. The Remiremont seismic swarm could reflect the occurrence of such mechanisms though mechanical and hydraulic property data for the surrounding and fault materials would be required to discuss further.

Fluid migration rate hypothesis
The present magnetotelluric profile can be considered as an instantaneous snapshot of the undergound apparent resistivity. Thus it is not possible to estimate any migration rate of the conductivity anomaly from the data of this study only. An estimate of the seismicity migration rate was previously obtained from the study of relocated seismic events for the Remiremont region [Audin et al., 2002]. They proposed a migration rate in the order of 5-10 km/yr over 30 km-long distance between 1980 and 1987 along some kind of fault zone. The seismicity pattern was proposed to reflect a zone of low permeability that would allow the propagation of transient pore pressure changes. After the main 1984 shock, a bilateral northward and southward earthquake pattern was also observed over a distance of 1-2 km in two months. From these migration rate, an estimate of the fault permeability in the order of 10 -13 to 10 -16 m 2 , was proposed [Audin et al., 2002] using characteristic time of diffusive process [Townend and Zoback, 2000]. Some permeability estimate were also obtained from the study of the micro-seismic events that occurred during a hydraulic fracturing experiment in Soultz-sous-Forêts geothermal site (150 km NE of Remiremont) [Bourouis and Bernard, 2007]. At this site, the Hercynian formations were affected by micro-seismic events on a characteristic horizontal distance. These events were considered as the trace of the fluid pore pressure diffusion in the fractured granitic medium. A permeability value of 10 -16 m 2 was estimated [Shapiro, 2000;Shapiro et al., 1999] which can be compared to the lower bound of the estimate proposed for the Remiremont region.
A new MT profile along the same transect as the profile of this study of the Rambervillers swarm would allow to quantify changes of the electrical resistivity structure over years and thus allow to determine any migration of the broad electrical conductive body (d). It would make possible to compare with the fluid migration proposed for the Corinth seismic swarms, the Remiremont earthquake and the Soultz-sous-Forêts micro-seismic events. This is however out of reach of the MT data acquired till now.

CONCLUSION
The electrical resistivity structure has been imaged in the Remiremont-Epinal-Rambervillers region in 2006 using the magnetotelluric geophysical method. The resistivity structure appears correlated to the distribution of earthquake hypocenters known in this region. A broad anomaly of high electrical conductivity has been revealed south of the 2003 Rambervillers earthquake. This region has been affected by replicas earthquakes over several years after the Rambervillers main shock. The presence of this anomaly and of the following long-term seismic activity is interpreted as the indication of the presence of a network of faults in the basement presumably dilated by fluid pressure, and siege of fluid pressure migration by diffusion processes along this network.