The Solfatara is one of the major volcanoes of the Phlegrean Fields (Campi Flegrei) volcanic complex, and it is located in a densely populated area a few kilometres west of the city of Naples. It is an active resurgent caldera that has been characterized by a rich history of surface–ground deformation and soil diffuse degassing and fumarolic emissions, which are indications of the top of a hydrothermal plume. A seismic survey was completed in May 2009 for the characterization of the main subsurface features of the Solfatara. Using the complete data set, we have carried out surface wave inversion with high spatial resolution. A classical minimization of a least-squares objective function was first computed to retrieve the dispersion curves of the surface waves. Then, the fitting procedure between the data and a three-sediment-layer forward model was carried out (to a depth of 7 m), using an improved version of the neighbourhood algorithm. The inversion results indicate a NE-SW fault, which is not visible at the surface. This was confirmed by a temperature survey conducted in 2010. A passive seismic experiment localized the ambient noise sources that correlate well with the areas of high CO2 flux and high soil temperatures. Finally, considering that the intrinsic attenuation is proportional to the frequency, a centroid analysis provides an overview of the attenuation of the seismic waves, which is closely linked to the petrophysical properties of the rock. These different approaches that merge complete active and passive seismic data with soil temperature and CO2 flux maps confirm the presence of the hydrothermal system plume. Some properties of the top of the plume are indicated and localized.
The Solfatara is the volcano in the Phlegrean Fields caldera volcanic complex where most of the volcanic-hydrothermal activity is at present concentrated. It is located in a densely populated area that is a few kilometres west of the city of Naples (Fig. 1a). The Phlegrean Fields magmatic system is still active, as was shown by the 1538 A.D. eruption at Monte Nuevo, along with its rich history of surface–ground deformation, known as bradyseism.
Since the Monte Nuovo eruption, the caldera floor has experienced long-term subsidence, with an average rate of 1.3 cmyr-’1 (Berrino et al. 1984; Dvorak & Mastrolorenzo 1991). The area has also sustained two recent unrest episodes, in 1970–1972 and 1982–1984, which were characterized by large ground deformation, changes in the geochemical signals at the surface and some seismicity (Madonia et al. 2008). The ground deformation changes during these two seismic episodes were remarkable, with the uplift amplitudes reaching 1.7 m in 1972 and 1.8 m in 1984; these were then followed by long-term subsidence (Todesco et al. 2010).
The volcanic features of the Solfatara crater consist of soil diffuse degassing, both inside and outside of the crater, and fumaroles that are mainly located on the eastern part of the crater. As indicated by Chiodini et al. (2001), about 3000 tons of a gaseous mixture of steam and CO2 are released per day at the Solfatara, and the near-surface condensation of steam produces a thermal power in the order of 130 MW per day. The Solfatara volcano is thus an area where degassing is concentrated, which indicates that it might correspond to the top of a hydrothermal plume (Chiodini et al. 2001; Bruno et al. 2007).
The geophysical structure below the Solfatara area has been investigated using different methods that have included gravity and electromagnetic surveys (Berrino et al. 1998; Bruno et al. 2007), and passive and active seismic tomography. Using seismic reflection data, Zollo et al. (2008) identified a 7-km-deep reflector over a 1-km-thick low-velocity layer that they defined as a partial melting zone that feeds the Phlegrean Fields volcano. At shallower depths, De Siena et al. (2010) used passive seismic attenuation tomography to investigate the hydrothermal system just beneath the Solfatara. A high seismic attenuation area with a low Vs/Vp ratio was interpreted as a CO2-rich structure that is located between 0 and 2500 m deep below the Solfatara. This thus confirmed that a hydrothermal plume stands at the centre of the Solfatara Crater.
The main issue of this study is to determine the main subsurface features of the Solfatara hydrothermal plume at high spatial resolution, using seismic experiments. A preliminary seismic study indicated important lateral variations in the elastic properties of the medium. For instance, Rayleigh wave velocities were measured at between 40 and 130 ms-’1 inside the crater (for a study area of ca. 150 m × 200 m and at frequencies ranging from 8 to 24 Hz). A second issue of this study is to determine the relationship between the seismic properties over the area and the spatial variations in the temperature, liquid and gas contents of the soil.
From 2009 April 25–27, 67 seismic sensors were deployed with a grid spacing of 20–25 m, which covered the 150 m × 200 m Fangaia area (Fig. 1b). Using a low-frequency broadband vibrometer collocated with each seismometer, an exhaustive surface wave data set was produced, from which a 3-D study of the seismic velocities was performed.
It should be noted that the Fangaia area comprises very different kinds of soils, from soft areas around the two bubbling mud pools in the centre, to the dry and solid clay in the southern part. It also shows a large range of surface temperatures, from boiling point to air temperature.
The seismometers were grouped into eight different seismic arrays, with each array consisting of up to nine vertical seismometers with their own acquisition system and global positioning system synchronization. During the day time, active experiments were performed using the vibrometer source. During the night time, seismic ambient noise produced by the hydrothermal activity inside the Solfatara was recorded. Temperature measurements were also performed on the same grid at 15 cm in depth, which were complemented by more recent temperature and CO2 flux measurements, which were carried out in 2010.
During the active experiments, a low-frequency vibrometer source was successively deployed at every point of the grid, whereas the 67 seismometers were recording. The signal S0(t) injected into the vibrometer source was a 60-s-long sweep with a frequency ranging from 10 to 50 Hz. For each one of the (67 × 67) source–receiver pairs (source s; receiver r), the received signal Sr,s(t) was cross-correlated with the input signal to produce the point-to-point impulse response Gr,s(t):
Before describing the inversion process itself, it is of note that the application of the reciprocity theorem between the source–receiver pairs provides a quality check of the data set, to efficiently select the best impulse responses for the dispersion curve estimation.
The reciprocity theorem states that the vertical response to a vertical force should be identical when the source and the receiver are exchanged. Considering the example of grid points 33 and 34 (see Fig. 1b), the signal recorded on the seismometer at point 34 for the vibrometer source at point 33 perfectly overlaps with the signal recorded on the seismometer at point 33 when the vibrometer source was set at point 34 (Fig. 2a). Before comparison, the two impulse responses, G33,34(t) and G34,33(t), were pre-filtered between 8 and 24 Hz; that is, the frequency band for which the signal-to-noise ratio (SNR) was maximum.
A global study of all of the source–receiver pairs (s, r) was carried out following the same procedure. This showed that 91 per cent of the pairs with s and r coming from the same sub-array (Fig. 1b) satisfied the reciprocity theorem, with a 0.9 correlation coefficient between Gs,r(t)and Gr,s(t) for a maximum allowed time shift of 0.01 s. In contrast, when s and r were chosen from among different sub-arrays, only 57 per cent of the pairs satisfied the reciprocity theorem with the same conditions. The non-applicability of the reciprocity theorem is mainly due to clock errors between the seismic arrays. Another reason is the larger distances (and thus weaker SNRs) between the source–receiver pairs when working with different sub-arrays. Note that in the frequency of interest, the signal attenuation in the Fangaia area was always greater than 40 dB for distances greater than 50 m. The inefficient ground coupling of the vibrometer source is also a possible explanation for some observed weak correlation coefficients.
The surface wave travel times were calculated for the source–receiver pairs that satisfied the reciprocity theorem criterion defined above, as well as the criterion of a SNR > 8. From the 67 × 66 = 4422 source–receiver impulse responses, 3020 signal pairs were extracted, which correspond to 1510 source–receiver pairs. The impulse response shown in Fig. 2a consists of several wave packets, which confirms the complexity of the medium.
In Figs 2b and 3, some of the frequency spectra are plotted in the 8–24 Hz frequency band for source–receiver pairs at different locations in the Solfatara. Again, the gaps and peaks in the spectra amplitude modulation are a clear indication of the spatial heterogeneity of the propagation medium.
The travel times obtained at each source–receiver pair were then inverted to produce local dispersion curves for the surface wave group velocity at every point of a 205 m × 160 m rectangular surface grid that covers the entire area, with a cell size of 5 m × 8 m chosen as twice the minimum wavelength (around 3 m). To do so, minimization of a least-squares objective function is used, taking advantage of the linear relationship between the slowness and the travel times :
The spatial correlation length L is assumed to be stationary in the whole medium, and it is estimated at 30 m (ca. twice the maximum wavelength). This weighting allows us to stabilize the inversion result, although it decreases the spatial resolution. In eq. (3), r is an empirical factor that controls the importance of the model smoothing relative to the data covariance matrix that is defined from weighting coefficients that depend on the quality of the empirical data, which were evaluated using the reciprocity theorem efficiency. Different inversions were computed for different values of r, with the data filtered between 8 and 24 Hz. The optimal factor r= 0.08 corresponds to the maximum curvature of the classical L-shaped misfit curve.
Note that the straight ray approximation provided 70 per cent of reconstructed variance for travel times. Curved ray inversion or, even better, finite-frequency tomography (based on the computation of the point-to-point Fresnel zone) could be applied to this data set to further refine the group velocity estimation.
Fig. 4 shows the group velocity map that was obtained at 15 Hz. There are high lateral variations between 40 and 130 ms-’1 with error bars on the order of 5–20 ms-’1. This is typical of hydrothermal systems (Vandemeulebrouck et al. 2010), where the water saturation, porosity and temperature strongly affect the surface wave velocity. Similar group velocity maps were then estimated at eight frequency values, ranging from 12 to 26 Hz, with a constant 2 Hz step. At each frequency value f, the source–receiver impulse responses Gr,s(t) were filtered between f-’ 5 and f+ 5 Hz, prior to applying the spatial inversion of the travel times, as described above.
In Fig. 4, the estimated dispersion curves are plotted at four particular grid points that were chosen in areas of interest. Point 1 is located in a mud area (hot water). Very weak velocities (50 ms-’1) can be noted, and especially a drop around 15 Hz that might be due to the presence of water. At points 2 and 4, where the soil is dry and hot and there is an elevated CO2 flux (>2500 gm-’2d-’1), the velocities are higher, they decrease with frequency, and they are spatially homogeneous. The velocities at point 4 have a maximum around 150 ms-’1; that is three times the minimum velocity at point 1, which is only 50 m away. Point 3 is located in a dry area that has no specific geological features visible at the surface, but which has unusual weak and constant velocities.
The last step of the inversion consists of matching the dispersion curves with a 3-D shear-velocity model. At each cell of the surface grid, the dispersion curves were inverted to obtain the corresponding 1-D depth-dependent shear-velocity profiles. A simple three-layer parameterization was sufficient to image the shallow structures of the Solfatara. Each layer was defined by its depth and the elastic parameters: P-wave velocity, S-wave velocity and density. For each layer, both the P-wave and S-wave velocities were inverted. The medium density was fixed at 1500 kgm-’3, and the maximum penetration depth was fixed at 7 m, according to the common law of the third of the maximum wavelength (around 15 m). The fitting procedure between the data and the three-layer forward model was carried out using an improved version of the neighbourhood algorithm (NA) (Wathelet 2008). In contrast to basic Monte Carlo sampling, the NA attempts to guide the random generation of new samples according to the results obtained for the previous samples. Voronoi cells were used to model the cost functions across the parameter space. As compared to classical NA (Sambridge 1999), this improved version better explores the non-uniqueness of the problem, due to continuous scaling of the parameter space. For each cell, 150000 models were computed. Since surface wave inversion is well known for possible non-unicity of the solution, we checked the model convergence to its optimum and we carefully explored the space parameter (not shown here). Conclusion is that only S-wave velocities were satisfactorily constrained through the inversion process down to a depth of 7 m.
In Fig. 5, the S-wave velocity map is plotted for three depths: 1.0, 3.5 and 6.5 m. This clearly shows the seismic heterogeneities of the area and the complexity of the petrophysical properties. In the (E) area in Fig. 5, there is a high-velocity zone that corresponds to a dry and compact zone of the Solfatara basin that has low attenuation (see Fig. 6). In the (G) area in Fig. 5, close to one of the two Fangaia mud pools, the presence of water at the top surface can be seen by the thin low-velocity zone of ca. 50 ms-’1, above a 100 ms-’1 layer. This structure suggests that there is a condensation layer in the first meter of the soil, as also observed on the temperature profiles of Chiodini et al. (2005). Other geophysical surveys will now be carried out to interpret the spatial anomalies in the velocity maps.
The relationship between attenuation and shear velocity
As with velocity, the attenuation of seismic waves in a volcanic structure is a relevant parameter for the identification of the petrophysical properties (e.g. permeability and porosity) and conditions (e.g. saturation and fluid liquid/gas ratio) of the porous medium. Indeed, according to White et al. (1975) and Brajanovski et al. (2006), fluid flow waves at shallow depth are responsible for dispersion and attenuation at low frequencies. Considering that the intrinsic attenuation is proportional to the frequency (Johnston et al. 1981), high frequencies are more attenuated than low-frequency waves. Quan & Harris (1997) proposed an approach to evaluate the attenuation through calculation of the frequency centroid downshift during wave propagation in the medium. Hence, for the Solfatara survey, a study of attenuation was carried out based on the centroid estimation, following the laboratory studies and centroid definition described by Lei & Xue (2009). For every source–receiver pair (i, j), with i and j such that the source–receiver distance does not exceed two grid points, the signal Si,j(t) is transformed into the corresponding spectrum , from which an estimate of the centroid fi,jc is obtained from the Hedlin formula:
In Fig. 6, the value of the calculated centroid is represented by the coloured lines and the dots between i and j. Following Quan & Harris (1997), with the assumption of a frequency-independent Q model, the centroid map can be used as observed data, to reconstruct the attenuation distribution. If the attenuation is high, the centroid should be weak, and vice versa. The centroid spatial distribution over the area shows large variations, which ranges from 17 to 22 Hz. This clearly highlights a strong attenuation area (F) at the centre of the Solfatara. The use of this centroid-based method provides several advantages compared to other attenuation estimation methods, as it is relatively insensitive to geometric spreading, reflection and transmission effects, as well as source and receiver coupling, radiation patterns and instrumental responses (Quan & Harris 1997).
The comparison between the centroid map and the velocity map for the shear waves shows a good match between the S-velocities and attenuation (which are directly linked to the centroid values). Two areas with opposite features are clearly seen: the centre of the Solfatara (Fig. 6, area F) has strong attenuation and low surface velocity, whereas the large southern area has low attenuation (high centroid) and high velocity (Fig. 6, areas B and E). As the central Solfatara area corresponds to low surface temperatures, we interpret this area as a mechanically soft and liquid-saturated part of the hydrothermal system, as compared to the drier and hotter zones in the southeastern part of the study area. It should also be noted that the (F) zone has the lowest elevation of the area, and it usually drains the run-off waters from the southern part of the crater. More generally, the spatial variation of the attenuation into the Solfatara crater is surprisingly large for such a small area of a few hundred meters in size that has comparable soil-matrix properties, which indicates that fluid properties and temperature have major roles in the elastic properties. These local effects need to be taken into account when dealing with larger structures, such as the Phlegrean Fields.
Connections with the geological structures
The velocity maps are next compared with the ground temperatures measured at 15 cm in depth, or more precisely, with the horizontal variations (gradient) of the ground temperatures (Fig. 7). At depth z= 6.5 m, a drop in the velocity values can be seen in Fig. 5c along a marked NE-SW direction; this same trend can also be seen in the horizontal gradient of the surface temperature (Fig. 7). This thus indicates a clear geological boundary; that is a fractured zone. The small shift between the location of this fault from when observed at the surface through the temperature to the velocity drop at z= 6.5 m might represent the fault dip at depth.
Localization of the ambient noise sources and relationship with CO2 degassing and surface temperature
During the 2009 acquisition, the ambient noise was recorded over two nights. Pre-processing was applied, which consisted of the elimination of high-amplitude seismic events by truncating the recording amplitude at three times the standard deviation of the seismic noise. Then, frequency equalization was performed, to whiten the noise spectrum in the 5–25 Hz frequency interval (Bensen et al. 2007). This pre-processing was performed to improve the noise properties, such that the sensor-to-sensor travel times can be extracted from the time-averaged noise correlation function (Gouedard et al. 2008). In practice, the noise correlation function was averaged over the two nights for every receiver pair. Unfortunately, despite some interesting comparisons on a small set of receiver pairs (mostly at short distances), the passive data set cannot be used to construct a 3-D shear velocity structure that could be compared to the active data result. In general, because of both strong attenuation and medium heterogeneity, the SNR of the cross-correlated functions was too weak to provide an accurate travel-time measurement.
However, the travel times of the dominant wave packet could be extracted from the correlation function for receiver pairs exhibiting a SNR > 8. These travel times were then matched to the theoretical surface wave travel times that were calculated for a surrogate point source in the Solfatara and the group velocity map obtained from the active data inversion (Fig. 4). Using a grid search for the surrogate point source, the noise-source probability map was constructed (Fig. 8b) as described in Cros et al. (2011).
In 2010, a temperature and CO2 flux survey was performed at the Solfatara. When the noise-source map is compared with the CO2 flux and surface temperature maps, a good match is seen between the maxima of these three parameters (Fig. 8). The CO2 flux shows three maxima (Fig. 8a, areas A, B, C) that are linked to high probabilities of noise-source location (Fig. 8b) and high temperatures (Fig. 8c). The temperature map also shows another area (Fig. 8c, area D) with a very high temperature (around 100 °C); this is linked to secondary noise sources (Fig. 8b) and CO2 maxima (Fig. 8a). Thus, the dominant ambient noise probably comes from CO2 degassing and high temperatures, as observed in Fig. 8, area (A).
Starting from the dispersion curves between the station pairs, the inversion method follows a classical minimization of a least-squares objective function, to retrieve the Rayleigh-wave group velocity map. We have improved the robustness of the inversion by the introduction of weighting coefficients that depend on the quality of the data, which was evaluated using the validation of the reciprocity theorem. An improved version of the NA is used finally, to estimate the S-wave velocity maps that show important spatial variations of the medium properties. A NE-SW fault is seen, which was not visible at the surface. Using passive seismic measurements, the source of the elevated ambient seismic noise in the area of Solfatara Crater is seen to be located at the eastern part of the area, which has a very high CO2 flux and high surface temperatures that are close to boiling. Several studies have shown that the physical processes that can generate powerful acoustic noise in multiphase hydrothermal systems like the Solfatara are related to the presence of bubbles and boiling, and the propagation of bubbles in conduits and cavitations (Kieffer & Sturtevant 1984; Kedar & Kanamori 1996; Vandemeulebrouck et al. 2005). In the case of the Solfatara, the bubbles probably contain both CO2 and steam, and the steam-to-liquid phase change that takes place close to the ground surface (Chiodini et al. 2005) will generate the ambient seismic noise. As observed in the mud pools, the CO2 bubbles can percolate the liquid layer formed by condensation, and can also produce remarkable acoustic intensity. The large amounts of heat at the surface come from the latent heat exchanged during the phase transition. This area represents the top of a large gaseous hydrothermal rising plume in the Solfatara crater, which probably contains a few others. The western part of the Solfatara corresponds to the descending part of the plume, with lower temperatures and lower CO2 flux. The seismic velocity and the attenuation tomography describe the surface boundary of a hydrothermal plume, which appears as a very heterogeneous medium in terms of its elastic wave properties, even on a small scale. High-resolution active and passive seismic measurements thus appear to bring additional constraints to gas flux and ground temperature surveys, to aid in the deciphering of the structures of complex multiphase hydrothermal systems.
Thanks to B. Dupuy, J. Grangeon, M. Radiguet and B. Vial for their work during the field survey, and G. Cougoulat for his technical assistance. We are grateful to G. Angarano and the staff of Vulcano Solfatara who kindly allowed us to make measurements at Solfatara in excellent conditions. This study was funded by CNRS-INSU-ST and BQR ISTerre.