CC-BY 4.0
1. Introduction
The characterization of the seismic hazard in volcanic regions stands as a huge challenge with real both scientific and societal issues. It implies understanding physical and chemical processes leading to seismicity, potential fluid transfers and volcanic activity. Unlike major high magnitude seismic events, earthquake swarms are characterized by the lack of a dominant main shock. Such swarm-like seismicity can be triggered by numerous physical mechanisms in various contexts, especially in active volcanic areas where they are typically associated with either hydrothermal fluids, magma movements or volcanic eruptions (Mesimeri et al., 2023). Earthquake swarms also occur in area without active volcanism in continental rifts (i.e. Rio Grande, West Bohemia, and recently in the French Massif Central; Ibs-von Seht et al., 2008) where seismicity is restricted to deep-reaching zones of weakness and a heterogeneous stress field that could favor fluid circulation in crustal layers. If the seismicity of intracontinental regions is generally low compared to active margins for example, seismicity is less constrained due to the relative rarity of occurrences and the absence of dense sustainable monitoring networks (Larroque et al., 2021). Regarding the long-term history and phenomena in such areas (significant long-lived structural heritage, tectonics and erosion, diffuse seismicity and very rare strong earthquakes), seismic hazard assessment, and especially the characterization of associated sources and mechanisms, remains a particularly difficult exercise.
The French Massif Central (FMC) remains a large seismogenic area in memory of its past and long-term volcano-tectonic activity (Battaglia and Douchain, 2017). For example, at the local scale of the Monts Dore volcanic province, several seismic events of low magnitude are regularly recorded, and an important seismic swarm recently occurred in 2021–2022, remarkable both in terms of duration and number of events detected. To provide a detailed spatial characterization of this seismic swarm event, the local scientific and technical staff of the Observatoire de Physique du Globe de Clermont-Ferrand (OPGC) and the Laboratoire Magmas et Volcans (LMV) in collaboration with the Bureau Central Sismologique Français–Réseau National de Surveillance Sismique (BCSF-RéNaSS) deployed an operational response on the field during the seismic crisis in order to provide at the first order structural insights around this phenomenon. The relocated seismicity has evidenced a seismic swarm clustered between 3 and 6 km-depth aligned following two fault planes-oriented along main regional tectonic axes. Magnetic anomalies mapping, 2D Electrical Resistivity Tomography (ERT) and magnetotelluric (MT) soundings, highly sensitive to fluids contents, fracturing and alteration processes, were carried out to image any structural contacts that could link fluids circulation or storage at various depths with structural weakness paths. Complementary soil CO2 and self-potential (SP) measurements have also been carried out to highlight the possible circulation of fluids close to the surface. Details on the spatio-temporal behavior of the seismic swarm as well as dynamic insights from complementary temporal surveys are the focus of a companion paper (Boudoire, Gailler, et al., 2025).
2. Context of the study
The recent seismicity in the French Massif Central is mainly distributed in four main areas (Battaglia and Douchain, 2017): (1) the Monts Dore sector, (2) the Saint-Flour one, (3) the region of Ambert (including a seismic swarm in 2013 and more recently in 2024) and to the northwest of Clermont-Ferrand on a large area encompassing the Combrailles zone (Figure 1a). It is moreover important to stress a virtual absence of seismicity under most of the volcanic provinces (Cantal, Velay, Devès, Chaîne des Puys), with the exception of the Monts Dore area. Actually, this volcanic province is also the site of the last eruptive activity known in mainland France around 6700 years ago (Pavin’s Group of Volcanoes; Villemant et al., 2016). Seismic swarms in the Monts Dore volcanic province (less than one hundred seismic events with an occurrence of a few days to a few weeks) were documented in 1979, 1980, 1983, 1984, 2008 and 2013 (ibid.). A new very recent seismic swarm has occurred in 2021–2022, highly remarkable by the number of earthquakes detected (more than 400 events on the BCSF-Rénass catalog; https://api.franceseisme.fr/fr/search; BCSF-Rénass, 2022) and for the longevity of seismic activity (i.e. more than one year). In more detail, this latter is concentrated on an 8-km2 area with the first events initially located 2 km west of the Beaune le Froid village and 2 km north-west of the Chambon-sur-Lac town (Figure 1b).
(a) General context of the Monts Dore volcanic province (colored areas) in the FMC and main tectonic features. The rose diagram indicates the main orientations of the four families of emblematic faults. (b) Zoom on the simplified geological and structural scheme of the study Monts Dore volcanic province (after Chèvremont et al., 2023).
Regarding its overall context, the Monts Dore volcanic province covers an area of approximately 500 km2 with about 220 km3 of emitted products (Brousse, 1971). The volcanic system is divided into two complex units: the Guéry stratovolcano (or “Mont Dore”; 3.1 to 1.8 Ma) and the Sancy stratovolcano (1.1 to 0.24 Ma; Cantagrel et al., 1999; Chèvremont et al., 2023). It is also the site of the last eruptive activity known in mainland France around 6700 years ago (Pavin Volcanoes Group; Villemant et al., 2016). Thermomineral spring waters are currently widespread in the area (Boudoire, Pasdeloup, et al., 2023; Bräuer et al., 2017; Ricci et al., 2024) that make it of peculiar interest for geothermal prospects (Pauwels et al., 1997; Varet et al., 1980). From a geological and tectonic point of view, the Monts Dore volcanic province is marked by sedimentary and volcanic formations outcropping on a highly fractured Variscan-age plutonic basement. It is located between the two regional major tectonic accidents of the Limagne gap to the east (Oligocene) and the Sillon Houiller to the west. According to Merle, Aumar, et al. (2023), the Variscan faults (>300 Ma) can be divided, based on their orientations, into four families named after emblematic faults of the FMC (Figure 1a):
- The Oligocene N20 direction (Sillon Houiller Fault family), which acted in distension, and mainly conditioned the regional tectonics in horsts and grabens (Hulot, 1988; Varet et al., 1980).
- The NE–SW direction (Cévennes Fault family), a Variscan intracontinental strike-slip fault that was reactivated during the Tertiary.
- The Pleistocene N160-170E direction (Villefort Fault family), which played a distensional role feeding the late phases of volcanism (Monts Dore, Chaîne des Puys; Hulot, 1988; Varet et al., 1980). It replays the Hercynian and Oligocene faults in shear and reactivates them until the present time as suggested by thermal and mineral activity along all these weak directions (Besson, 1978; Varet et al., 1980).
- The late-Hercynian N110-140E direction (Armorican Shear Zone family) exhibits varying dips depending on their behavior, ranging from high-angle (sub-vertical) for strike-slip faults to steeper dips for normal faults.
These brittle directions have been characterized by recent reactivation (i.e. Tertiary and Quaternary), accommodating extension and uplift processes (Merle, Aumar, et al., 2023; Varet et al., 1980). The reactivation of these structures likely continues at present, considering the ongoing uplift of the FMC (Delfau and Lenôtre, 1992) and the presence of thermal activity emerging notably through these fractures (Besson, 1978; Varet et al., 1980). While these major discontinuities have been highlighted at a regional scale, it is noteworthy that they are also found at smaller scales (Aumar, 2022; Merle, Aumar, et al., 2023), consisting of multiple faults maintaining a general orientation rather than a single linear fault (Merle, Aumar, et al., 2023). Accordingly, around the seismic swarm area, several interesting structures have been identified in the literature, which may belong to the families defined by Merle, Aumar, et al. (ibid.) (Figure 1b):
- The Jassat strike-slip fault (Chèvremont et al., 2023) that is part of the Armorican Shear Zone family, whose faults network forms the boundary between the Monts Dore and the Plateau des Dômes (Merle, Aumar, et al., 2023).
- Closer to the swarm, around Chambon Lake and in the western part of the study area, a fault network was defined by direct observation or gravity measurements (Glangeaud et al., 1965) with a consistent orientation to the Sillon Houiller Fault family. A similar oriented faults network is also observed to the east in the substrate outcropping at Saint Nectaire (Veyre-Monton geological map edited by the BRGM; Chèvremont et al., 2023; Figure 1b).
Furthermore, the presence of numerous recent volcanic vents (cones, diatremes, maars; Chèvremont et al., 2023) evidenced along a ∼N–S alignment in the center of our study area suggests an organization of eruptive vents along Variscan fault structures, as proposed by Merle, Aumar, et al. (2023) in the Chaîne des Puys. A late feeding of Monts Dore volcanic activity along the N170 direction is also suggested by Varet et al. (1980).
3. Methods
In such a complex tectonic province and well-developed hydrothermal system, several hypotheses can be proposed to explain the seismic activity: (1) a tectonic readjustment along the previously described faulting directions and/or (2) a disturbance of an underlying hydrothermal system by endogenous (arrival of deep fluids) and/or exogenous (drainage) phenomena. To better characterize its origin and address potential issues raised by a renewal of the seismicity in local volcanic provinces as a whole, several actions have been initiated under the aegis of the Pôle Régional d’Observation de l’activité Volcano-tectonique d’Auvergne et d’Ardèche (PROVA2; Boudoire, Gailler, et al., 2025). Initially set up for structural investigations, the team carried out a series of geophysical and geochemical measurements at the scale of the area considered in response to this exceptional seismic activity (Figure 2).
(a) Diagram summarizing the main goal of each approach and their complementarity in the study context. (b) Location of the geophysical measurements superimposed to the relocated seismic events (zoom on area located in Figure 1).
We mainly focus on contrasts in magnetization deduced from magnetic anomalies and electrical resistivity at various scales since those parameters depend on numerous factors (mechanical heterogeneities, alteration, fluids storage and circulation) that accumulate in time, especially in long-dormant volcanic provinces. Beyond the overall area concerned, several key sites were also investigated in more detail. The purpose is to better characterize the overall context at various scales (both laterally and in depth), in order to bring additional constraints on the origin of the seismic crisis, and its potential for coupling tectonics with hydrothermal circulation at various depths.
3.1. Seismic data
The Monts Dore region experienced a significant seismic episode starting from late 2021. BCSF-Rénass detected the first events on November 16 using RESIF-Epos (now Epos-France) seismological stations. Seismic activity intensified from March 19, 2022, reaching 441 earthquakes detected and localized by BCSF-Rénass by year-end, with a peak of about 30 events recorded on July 24 and on November 19, 2022 (Figure 3a). In order to better characterize the spatial distribution of seismic events and thus the evolution of the swarm, three additional temporary seismic stations were deployed from April 4, 2022, together with a network of seismic nodes.
(a) Number of earthquakes per month (BCSF-RéNaSS vs. Fine-tuning localizations workflow Grunberg and Lambotte, 2024) with visualization of installations of temporary nodes and additional seismic stations. (b) Relocation of seismic events using the double-difference method, highlighting principal directions that distinguish two main families of seismic events. (c) Time evolution of the hypocenters related to the seismic swarm (relocated data BCSF-RéNaSS).
To enhance the quality of this seismicity catalog, BCSF-RéNaSS incorporated artificial intelligence tools, specifically an automatic phase picking method based on deep learning. The Seisbench/Phasenet model (Woollam et al., 2022) was employed, trained on manually phase picked from the BCSF-RéNaSS catalog covering 2016 to early 2024. This approach significantly improved the detection and localization of seismic events, particularly for low-magnitude earthquakes previously undetectable using SeiScomP picker (SeisComP, 2008).
For earthquakes previously localized, manual picks were supplemented with automatic picks obtained through deep learning, using waveforms from all Epos-France stations within a 150 km radius. Data from a temporary network specifically deployed to monitor this seismic swarm was also incorporated. This network included six short-period seismic nodes, operational between 30/06/2022 and 08/08/2022, along with three broadband stations whose data are retrieved since 13/04/2022 (XX.CHBN), 06/16/2022 (XX.BNLF), and 06/22/2022 (XX.VRNF).
The catalog construction followed several key steps. First, picks closely spaced in time were clustered to reduce redundancy, prioritizing manual picks. Next, seismic phase association was performed to form events, combining the HDBSCAN algorithm (McInnes et al., 2017), to regroup picks close in time and space, with PyOcto (Münchmeyer, 2024), which discarded picks not conforming to typical travel-time curves. The third step involved event localization using the NonLinLoc algorithm (Lomax et al., 2000) with a five-layer 1D regional velocity model (Mazabraud et al., 2005) suited for FMC. Finally, when possible, moment magnitude (Mw) was estimated through waveform spectral fitting using a modified version of SourceSpec (Satriano, 2021). In order to further refine earthquake locations, a double-difference relocation was conducted using the HypoDD algorithm (Waldhauser and Ellsworth, 2000). This method incorporated both catalog-based travel times and differential times derived from waveform cross-correlation, minimizing relative travel-time residuals between nearby events. By improving the relative positioning of clustered seismicity, this approach enabled a more accurate reconstruction of fault structures.
The combination of deep learning for automatic phase picking, traditional seismological methods, and data from temporary seismic stations significantly improved the quality and completeness of this seismicity catalog. The refined catalog now includes 1616 earthquakes, compared to the 441 previously identified by BCSF-Rénass from November 2021 to the end of 2022. Additional details regarding the temporary stations are described in the companion paper (Boudoire, Gailler, et al., 2025).
Final relocation of seismic events (as discussed hereafter) using the double-difference method, highlighting principal directions that distinguish two main families of seismic events, as well a time evolution of the hypocenters related to the seismic swarm are presented in Figures 3b and 3c respectively.
3.2. Associated geophysical and geochemical one-off surveys
3.2.1. Magnetic field survey
A ground magnetic survey was performed in 2022 (Figure 4a) to derive a high-resolution magnetic anomaly map at the scale of the seismic area (i.e. about 25 km of cumulated profiles; Figure 2b) using a portable Overhauser proton-precession magnetometer (GEM GSM19) in walking mode (0.5 s sampling interval), with a sensor mounted on a backpack at a height of about 1.60 m above the ground. The absolute instrumental accuracy is commonly specified at ±0.1 nT, and the position was determined simultaneously while surveying with an integrated GPS with a fairly precise location of the magnetic measurements (i.e. around a few meters) even without any post-processing. To provide a more homogeneous and detailed map at the scale of the overall area of interest, we also used data extracted from the Bureau des Recherches Géologiques et Minières (BRGM) helicopter-borne geophysical campaign—HADGR211022 project, carried out in 2023 at an altitude of 80 m above the ground (Figure 4b). The overall dataset of this latter is available at https://data.gouv.fr/fr/datasets/donnees-geophysiques-aeroportees/ where more details can be found (so called E-sector in the data statement). This large-scale survey is conveniently oriented N–S, which is complementary to our further E–W oriented ground survey (Figures 2b, 3 and 4a). Note that stitching together such overlapping surveys acquired at different times, altitudes and spatial resolutions should be carefully considered. For this, we use a dedicated algorithm based on a Rectangular Harmonic Analysis (RHA) described in Thebault and Gailler (2023), with which field measurements are automatically corrected for the global main field (CHAOS-7 model; Finlay et al., 2020) to derive an interpolated homogeneous magnetic anomaly grid leveled at a constant altitude above the ground, i.e. 80 m (Figure 4c). The final interpretative dataset (Figure 4c) was interpolated using the minimum curvature method in order to derive the smoothest possible grid that fits the given data values. The method is particularily suitable (1) for data from both regional and local surveys, and (2) is more efficient when data are randomly distributed, as our present case. It was defined by how closely the solution surface matches the initial values, for a final interpolation versus extrapolation of 20 m in cell size and a blanking distance of 260 respectively (Figure 4d). Here again, the goal was to derive a first order and qualitative regional map.
(a) Total magnetic field extracted from the 2023 BRGM helicopter-borne survey (HADGR211022 project). (b) Detailed magnetic ground measurements performed in 2022. (c) Combined map of magnetic anomaly computed at an altitude of 80 m above the ground. (d) Reduced to the pole magnetic anomaly computed using the averaged magnetic field vector at the time of the survey. White text indicates village names. (e) Analytic signal computed from the RTP map. Seismic events (black spheres) and the main reported fault (dashed line) are shown in background.
A Reduction To the Pole (RTP; Figure 4d) of the magnetic anomalies was lastly applied to reduce the dipolar appearance of the anomalies and shift them to lie above their sources (Baranov, 1957). This transformation requires the directions of the apparent magnetization and the ambient field to be specified (Inclination-I and Declination-D), here the induced magnetization vector collinear with the ambient field (i.e. the averaged magnetic field vector at the time of the survey: I = −61.3° D = 1.8°) was used. Note that such an assumption remains acceptable since the study area does not span the Brunhes–Matuyama reversal as detailed hereafter. In order to reveal major magnetization contrasts we also computed the analytic signal (Figure 4e), a post-processing that transforms the shape of a magnetic anomaly to a positive one centered above the source (Nabighian, 1984; Nabighian, 1972). The absolute value of the analytic signal is obtained by taking the square root of the squared sum of the vertical and two horizontal derivatives of the magnetic field. The analytic signal was computed on the RTP grid and upward continued at 100 m in order to reduce the noisy very short wavelengths of the magnetic field commonly observed in such a volcanic and anthropic area. Independent of the direction of both the ambient magnetic field and the source magnetization, the RTP exhibits maxima over magnetization contrasts and is helpful for highlighting mechanical lithologic variations. Magnetic data treatment, RTP and analytic signal maps were computed using Oasis Montaj software (©Geosoft).
3.2.2. ERT-SP-CO2 reference profile
For a more detailed investigation, a 5 km long reference profile (red line in Figure 2b) was conducted in the northern part of the seismic swarm including Self-Potential (SP, i.e. passive measurements of disturbances in electrical currents relevant in mapping subsurface preferential flow paths) as together with soil CO2 measurements (50 and 100 m spacing respectively). Soil CO2 measurements were performed using a MEGA multi-analyzer designed at the Istituto Nazionale di Geofisica e Vulcanologia (INGV)—Sezione di Palermo, equipped with two Gascard NDIR CO2 sensors (Boudoire, Grassa, et al., 2020). The “dynamic concentration” method was applied for this survey (Gurrieri and Valenza, 1988). An Electrical Resistivity Tomography (ERT) profile was also acquired along a collocated but more restricted profile in an area where soil CO2 and magnetic anomalies had been previously detected. The measurements were performed using an ABEM Terrameter LS2 resistivity meter with 64 electrodes 10 m spaced, and a roll along made up of a half of the total cable length to increase the investigated profile length (i.e. 630 + 315 m long). We used a Schlumberger field protocol which has the advantage of being sensitive to both horizontal and vertical structures, down to a depth of about 100 m below the surface in its central part, given the 10 m spacing between electrodes. After an initial data analysis and filtering of residual spurious measurements, we computed a 2D inversion model of the electrical resistivity distribution at depth using RES2DINV commercial software (Loke and Barker, 1996).
3.2.3. Magnetotelluric soundings
For greater depth investigation (i.e. deeper than 8 km), magnetotelluric soundings were carried out at several key sites (Figure 2) over 4 to 5 days each, with two main goals: (i) to characterize the resistivity of the formations at the depth of the seismic swarm and to image any potential presence of crustal fluids (Patel et al., 2020; Pearce et al., 2020), (ii) to record any temporal variation that could be linked with the seismic activity. Since we focus here on structural insights, we do not discuss this latter point in the present study. The survey was conducted from December 2022 to February 2023 using a broadband Phoenix Geophysics instrument (MTU-5C). The station was successively installed in strategic areas regarding the seismic swarm: three above the main area of seismicity, one in its northern part, and the last one outside (south-west) where previous MT measurements were available in Vallée de Chaudefour (Figure 2b; Hulot, 1988) in order to ensure the consistency of resistivity measurements through time.
Data processing was carried out in a single station using so-called “robust” algorithms (Egbert and Booker, 1986). Static-shift effects that commonly affect MT (Pellerin and Hohmann, 1990) data were corrected by conducting simultaneous and collocated electromagnetic soundings (TDEM for Time Domain ElectroMagnetic). A careful analysis of the phase tensor was also carried out as discussed hereafter, allowing to graphically determine the strike axis of the main regional conductive structure. In terms of modelling, 2D inversions were then performed to model the best resistivity distribution along the main area of interest, i.e. a profile selected perpendicularly to the main fault of the area (grey line in Figure 3a). In order to maximize the impedance tensor’s off-diagonal elements and minimize its diagonal elements, data were rotated to the strike axis (N140). The 2D inversions were run with the RLM2D code (Rodi and Mackie, 2001) that uses a finite-difference modeling algorithm based on the approach of Tikhonov regularization (Tikhonov and Arsenin, 1977). The model was segmented into 2D rectangular blocks measuring 25 by 25 meters at the surface and in the central part of the profile, from which the block dimensions increased laterally by a factor of 1.5 and vertically by a factor of 1.06 with depth. In terms of initial constraints, an a-priori model calculated using a constant conductivity value of 100 Ω⋅m has been fixed as an initial step, and resistivity bounds have been also set from 5000 to 1 Ω⋅m, which is consistent in such geological contexts. Several methodologies were tested in terms of data rotation, providing similar results so we consistently used the best fit on inversion using the regional strike rotation (i.e. N140).
4. Results
4.1. Details on the seismic swarm and tectonic implications
The 3D geometric processing of the relocated events, the temporal analysis of their succession and the calculation of a few focal mechanisms, allow us to envisage the effect of a “flower structure” type tectonic system. Two main faults orientations associated with two main families of seismic events can be therefore distinguished (Figure 3a):
- A N173-85W dextral strike-slip fault including the deeper seismic events at a depth between 3900 m and 5600 m below the surface (Figure 3b). This fault corresponds to a shear zone approximately 350 m wide.
- A N13-68W normal fault located above the shear zone, including the shallower seismic events at a depth between 2600 m and 4400 m below the surface (Figure 3b).
Note that the dextral nature of the deep shear zone is not solely based on focal mechanism calculations, but is also supported by field evidence (in the form of cataclasites within the Variscan basement and the formation of Cenozoic sedimentary basins), as well as by rock mechanics calculations that are note detailed here. These calculations further allow for the determination of the maximum horizontal regional stress, which had previously been interpreted only in subsurface. The same applies to the normal faulting observed at the surface.
This tectonic system is in line with the regional fracturing of the Variscan basement, made up here of plutonic and metamorphic rocks. The shear zone thus appears as the conjugate fault of the Jassat fault (Vidal et al., 1996). Finally, it should be stressed that this seismic swarm confirms the current tectonic stress in the northern Massif Central, with a direction of around N148 for the maximum horizontal component, in agreement with the homogeneous maximum horizontal component at N150 recorded below a depth of 500 m in the Chassolles, Echassières, Mayet de Montagne and Auriat boreholes (ibid.).
4.2. Qualitative description of the magnetic anomaly map
Basically, spatial variations in magnetization mainly depend on rocks ages (i.e. if the concerned area spans a magnetic reversal), composition, fluids contents, mechanical heterogeneities (i.e. fracturing and destabilization structures) and associated alteration processes (e.g. Bouligand et al., 2020; Gailler and Kauahikaua, 2017) but are also largely affected by anthropic artifacts that should be considered. The combined airborne (Figure 4a) and ground (Figure 4b) magnetic anomaly map (Figure 4c) clearly highlights the strong effect of the main villages (especially Beaune-le-Froid and Bressouleille) that generate “anthropic” negative anomalies. However, an interesting long wavelength negative anomaly axis oriented about N170 appears in the central part of the study area, where the seismic swarm occurred. In long-dormant and complex volcanic provinces, numerous phenomena could generate such negative magnetic anomalies. Regarding the low elevation of the combined magnetic dataset as well as the overall small wavelength of the magnetic signal, we will consider that the sources of the associated anomalies are concentrated within the shallow pile of the study area, which does not span the negative reversal.
If thermomagnetic effects are to be proscribed in such an ancient volcanic context, stress variations (due for example to fracturing and hydrothermal activity) are commonly associated with piezomagnetic effects that could contribute to rocks demagnetization and therefore negative magnetic anomalies. Electrokinetic variations due to fluid circulation (e.g. Lénat et al., 2012) as well as slow natural or hydrothermal alteration processes would also locally affect the magnetization. The observation of an anomalous magnetic axis is even more striking on the RTP (Figure 4d) and analytic signal (Figure 4e) maps. There, the highest values of the analytical signal delimit a N140 to N170 corridor encompassing the seismicity. Discrete increase of the analytical signal is also associated with the family of seismic events aligned along the N13-68W. These maxima therefore suggest the presence of magnetization contrasts probably associated with lithological and/or mechanical heterogeneity at depth where the seismicity may (at least partly) cluster.
4.3. Combined analysis of the ERT-SP-CO2 reference profile
SP and soil CO2 measurements highlight two main trends along the northern reference profile (Figure 5a; red line in Figure 2b): (i) positive SP anomalies correlating to highest soil CO2 contents (well above 1.0%) and (ii) negative (or low) SP anomalies in areas where soil CO2 content is low (mostly lower than 0.2%). These trends can be attributed to two different behaviors in fluid circulation. In such a volcanic context, SP anomalies are mainly linked to electrokinetic phenomena, water–rock interactions leading to positive charges in the direction of flow (Revil, Naudet, et al., 2003; Revil, Pezard and Glover, 1999; Revil and Pezard, 1998). Therefore, negative SP anomaly, without notable soil CO2 degassing can be attributed to meteoric infiltration, the water buffering the CO2 content in the soil. A similar process was described by Boudoire, Finizola, et al. (2018) coupled with carbon isotopes fractionating by a Rayleigh condensation process. Such influence of meteoric water infiltration along the CO2 channels has been also referenced in different areas (e.g. Flechsig et al., 2008). Conversely, positive SP anomalies coupled with high soil CO2 content argue in favor of a rise of fluids (see Revil, Finizola, et al., 2023, and references therein for a review).
(a) SP and soil CO2 measurements along the 5 km long reference profile located in Figure 2b. (b) Resistivity distribution along the ERT profile, along with SP, soil CO2 and magnetic signals. Light shaded vertical rectangles highlight the main area of SP and soil CO2 positive anomalies. Dark shaded rectangles highlight the areas of negative RTP magnetic anomalies.
Electrical resistivity being highly sensitive to fluid circulation, the ERT profile can be used to reveal inherited weakness paths at depth that could favor an uprise fluid circulation highlighted by both soil CO2 and SP positive anomalies at the surface (Figure 5b). Similarly to the qualitative analysis of magnetic anomalies, the highly resistive layer observed at the surface (>500 Ω⋅m) can be attributed to massive, stable, and unaltered formations as detailed in various contexts (Portal et al., 2019). Two main conductive structures are evidenced deeper, in the western and central parts of the profile. This latter seems to rise toward the surface along two vertical conductive axes east of the profile. Note that these two axes are associated, at least in relative, to short wavelength but high amplitude negative magnetic anomalies (Figure 5b). Accordingly, the correlation between conductive and probably demagnetized axes rising toward the surface in areas where both SP and soil CO2 positive anomalies are measured could argue for a preferential rise of fluids from shallow weakness paths inherited from fracturing, alteration, and hidden structures at depth.
4.4. Deep information from magnetotelluric (MT) investigation
The analysis of the MT data using the phase tensor method (Figure 6a) reveals a N170 orientation at depth below the study area with high phase tensor value suggesting the presence of a conductor. This direction is fully consistent with the main trend evidenced by the magnetic surveys and by the family of the deepest seismic events (oriented N173-85W).
(a) Map showing phase tensor of the four MT stations located in the seismogenic area (black dots), and location of the 2D profile used for inversions (black line). (b) 2D cross-section along Peyre Lévade-Coujat soundings (West and East, respectively), and associated fit for each sounding. The black dots represent the relocated seismicity within 400 m around the profile. Northing (m) and Easting (m) WG84 UTM zone 31N and iso-altitude (m).
The resistivity distribution on the best 2D inversion result (i.e. RMS = 2) along a selected profile crossing the main regional inferred structures highlights three main features (Figure 6b): (i) In the central part of the model, a conductive sub-vertical body dips first to the west and then to the east as depth increases (indicated by iso-contour 24 Ω⋅m on Figure 6b). (ii) This conductive structure is encompassed by resistive borders (100–1000 Ω⋅m). (iii) A deep conductive body (with resistivity values lower than 10 Ω⋅m) is evidenced at depth (deeper than 6 km). Based on the geology of the area, even if the resistive surrounding could be interpreted as granitic formations that are more or less altered and water-saturated (i.e. between 300 and 1000 Ω⋅m according to Hulot (1988) and Dupis et al. (1980) respectively, who studied the area), the nature of the conductive bodies needs to be discussed in more detail. Another interesting feature is that the distribution of seismicity around the profile appears consistently located at the resistive/conductive contact, as it will be discussed hereafter.
5. Discussion and conclusions
5.1. Links between seismicity and local tectonics
The magnetic anomaly maps, especially the long and short wavelengths, observed in the analytic signal, clearly evidence some limits of mechanical heterogeneities (mostly N140-to-N170) in accordance with regional tectonics (Figure 4). The presence of mechanical and lithological weakness paths inherited from fracturing and alteration is confirmed at shallow depth by the combination of ERT, SP and soil CO2 measurements (Figure 5). The 2D magnetotelluric model shows the extension of these shallow heterogeneities in greater depth with a main conductive axis dipping west and extending down to 5–7 km-depth. Both geophysical and geochemical surveys thus evidence the presence of a significant anomalous structure that extends in depth.
We therefore derive a synthetic interpretation based on geological constraints and our multi-method approach in order to propose a structural scheme for this seismogenic area (Figure 7). This synthetic interpretation is refined by Boudoire, Gailler, et al. (2025) taking into account the temporal evolution of the processes involved in the 2021–2022 seismic swarm in the Monts Dore volcanic province. Our analysis of geological data indicates that the main shallow conductive structure could be associated with the Variscan Jassat strike-slip fault (Figure 7). Chèvremont et al. (2023) described this formation as composed of ultra-mylonite rocks with a N140-85SW foliation that is consistent with our geophysical observations. This is also consistent with the main long wavelength evidenced in the analytic signal, which could mark such a shear zone. It should be noted that Variscan faults are known to have a control on the recent volcanic activity in the Monts Dore and Chaîne des Puys volcanic provinces (Merle, Aumar, et al., 2023; Varet et al., 1980). Actually, the eastern border of the ultra-mylonite corridor is marked by numerous volcanic manifestations (i.e. springs, maars, flows and domes according to the geological map; Figure 1b). These findings suggest that, even nowadays, the Variscan faults, such as the Jassat fault in the present study (Figure 7b), are still major paths for channeling fluids from depth to the surface.
Combined synthetic schemes of structural insights derived from our multi-methods analysis in the area of the 2021–2022 seismic swarm. The rose diagram recalls the main orientations of the four families of emblematic faults.
5.2. A seismic swarm triggered by fluids circulation?
The nature of the deeper part of the conductive body evidenced at about 8–10 km in depth is more conspicuous. Hulot (1988) evidenced and interpreted a deep conductor at a depth of about 10 km below the Sancy volcano as an ancient magma chamber. The recent volcanic events (i.e. springs, maars, flows and domes; Figure 1b) along the Eastern border of the fault, such as the recent Tartaret eruption (about 15 ka; Figure 1), could support the idea of the existence of an ancient magma chamber, possibly a mushy zone area in relation with this deep conductor body. Further investigations are required to decipher the depths of the magma ponding zones and the eruptive dynamics related to this eruption. However, the current gas and water chemistry from thermomineral springs in the Monts Dore volcanic province does not support the idea of current magma storage and degassing at crustal depth (Granet et al., 1995) as it may be observed for other active volcanoes (Boudoire, Gailler, et al., 2025; Bräuer et al., 2017; Pauwels et al., 1997; Ricci et al., 2024).
Another explanation for the deep conductor would be the existence of a fluid reservoir at this depth. We observed that the seismic activity occurred at the contrast between high and low resistivity formations that is characteristic of fluid injection-related seismicity as proposed by Pearce et al. (2020, and references therein). These authors suggest that this phenomenon is linked to the migration of fluids in a permeable zone (mechanically weak, with low resistivity) adjacent to an impermeable zone (mechanically strong, with high resistivity), resulting in the accumulation of fluid pressure and in successive rock failure and/or fault reactivation. Based on the data presented here, we hypothesize that the impermeable granitic borders and permeable discontinuity of the Jassat strike-slip fault allow hydrothermal fluids to channel and accumulate pressure, leading to the fault reactivation by reducing the effective normal stress (ibid.).
This hypothesis is also sustained by studies in similar contexts such as the Nový Kostel region (Czech Republic) located in the Eger Graben (Fischer, Horálek, Hrubcová, et al., 2014), that has experienced numerous seismic swarms since the 1980’s. The most recent events show a clear trend of vertical migration of seismicity in time and space, interpreted as a progressive rupture along a major fault plane (Bachura et al., 2021; Fischer, Matyska, et al., 2017; Hainzl and Fischer, 2002) consistent with the regional stress field (Fischer, Horálek, Michálek, et al., 2010). To explain them, Fischer, Horálek, Hrubcová, et al. (2014) reject the tectonic hypothesis, because the accumulation of stress sufficient to cause such recurrent activity seems very unlikely regarding the tectonics of Central Europe. These authors propose that these seismic swarms result from the reactivation of faults, caused by the migration toward the surface of CO2-rich fluids coming from a deep mantellic reservoir where they have accumulated pressure. This model also explains periods of quiescence between swarms as reflecting the time needed for the “deep reservoir” to re-pressurize after a rupture. Such an interpretation is relevant in our present case, as the seismic swarms in the Monts Dore volcanic province are very similar to those in the Nový Kostel region, in terms of geodynamic context (Merle and Michon, 2001), recurrence, duration, magnitude, number of events, and behavior. It suggests that magmatic processes are still active in the mantle beneath the Monts Dore volcanic province as supported by the isotopic signature of the gas emissions in the area (Bräuer et al., 2017; Boudoire, Pasdeloup, et al., 2023; Boudoire, Gailler, et al., 2025).
5.3. Concluding remarks
According to a detailed seismic analysis combined with (1) geological information and (2) complementary geochemical and geophysical surveys, we presented a spatial definition of the 2021–2022 Monts Dore large seismic swarm recorded beneath the most recent stratovolcano of the French Massif Central.
We evidenced its link with the reactivation of pre-existing tectonic structures and suggested the influence of a circulation of fluids at depth as commonly inferred in such a complex structural context. However, the origin of this reactivation is still subject to question and requires further evaluation of time series of monitoring records from both permanent and temporary networks/surveys. The companion paper of Boudoire, Gailler, et al. (2025) aims at better deciphering the dynamic processes involved in this particular seismic zone, through highly complementary information derived from temporal surveys.
Declaration of interests
The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.
Open research
The ground geophysical data (i.e. magnetic field measurements, ERT profile and MT soundings) presented in this study are available at https://doi.org/10.25519/PROVA2.
The helicopter-borne magnetic data are extracted from the BRGM HADGR211022 project and available at https://data.gouv.fr/fr/datasets/donnees-geophysiques-aeroportees/ (ETALAB-Licence-Ouverte-v2 0 ID 65412d94ef687b2a614553f8).
Licensed commercial softwares have been used for geophysical data modeling and visualization: EmPower (2022), Geotools (2022), Geosoft (2022), Golden Software Surfer (2015), Matlab (2020) and Seequent Limited (2020).
Funding
This work was partly funded by the OPGC through the “Pôle Régional d’Observation de l’activité Volcano-tectonique d’Auvergne et d’Ardèche: The Auvergne-Ardèche regional volcano-tectonic monitoring center” (PROVA2), the FEDER Project “CAPRICE”, and the CNES in the framework of the projects Volcadrone and RMagEDrone “Multi-scale geophysical imaging of volcanic edifices”.
The HADGR211022 project was financed by ERDF funds for the Auvergne and by BRGM.
This is contribution no. 676 of the ClerVolc program of the International Research Center for Disaster Sciences and Sustainable Development of the University of Clermont Auvergne.
Acknowledgements
We are grateful to the Town halls of Chambon sur Lac, le Vernet Sainte-Marguerite and Murol as well as all the landowners. We also acknowledge Georeva and Phoenix Geophysics for their scientific and logistic support as well as the 2023–2024 Master 1 promotion, Adèle, Maxime, for their help during the ERT profile acquisition.
We would like to thank the tho anonymous referees for their constructive comments.