Integrated tomographic methods for seismic imaging and monitoring of volcanic caldera structures and geothermal areas

    Integrated tomographic methods for seismic imaging and monitoring of volcanic caldera structures and geothermal areas O. Amoroso, G. Festa, P.P. Bruno, L. D’Auria, G. De Landro, V. Di Fiore, S. Gammaldi, S. Maraio, M. Pilz, P. Roux, G. Russo, V. Serlenga, M. Serra, H. Woith, A. Zollo PII: DOI: Reference:

S0926-9851(17)31040-6 doi:10.1016/j.jappgeo.2017.11.012 APPGEO 3374

To appear in:

Journal of Applied Geophysics

Received date: Revised date: Accepted date:

11 May 2017 14 October 2017 24 November 2017

Please cite this article as: Amoroso, O., Festa, G., Bruno, P.P., D’Auria, L., De Landro, G., Di Fiore, V., Gammaldi, S., Maraio, S., Pilz, M., Roux, P., Russo, G., Serlenga, V., Serra, M., Woith, H., Zollo, A., Integrated tomographic methods for seismic imaging and monitoring of volcanic caldera structures and geothermal areas, Journal of Applied Geophysics (2017), doi:10.1016/j.jappgeo.2017.11.012

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Integrated tomographic methods for seismic imaging and monitoring of volcanic caldera structures and geothermal areas

IP

T

O. Amoroso1, G. Festa1, P. P. Bruno2, L. D’Auria3,**, G. De Landro1, V. Di Fiore4, S. Gammaldi1, S. Maraio5, M. Pilz6, P. Roux7, G. Russo1, V. Serlenga1,**, M. Serra1, H. Woith8, A. Zollo1

Department of Physics “Ettore Pancini” University of Naples ‘Federico II, Italy.

2

Petroleum Geoscience Dept, The Petroleum Institute, Abu Dhab

3

Istituto Nazionale di Geofisica e Vulcanologia, Osservatorio Vesuviano, Italy

4

IAMC Napoli, Italy

5

Centro di GeoTecnologie - Università di Siena, Italy

NU

SC R

1

6

Swiss Seismological Service, Swiss Federal Institute of Technology, Sonneggstr. 5, 8092 Zurich, Switzerland 7

GFZ German Research Center for Geosciences, Helmholtzstr. 5, 14467 Potsdam, Germany

D

8

MA

ISTerre, Université Grenoble Alpes Joseph Fourier, Rue de la Piscine 1381,F-38058, Grenoble, France

*

CE P

TE

Now at Instituto Volcanológico de Canarias, Tenerife, Spain. ** Now at Consiglio Nazionale delle Ricerche, Istituto di Metodologie per l’Analisi Ambientale, Tito, Italy

Corresponding author: Ortensia Amoroso ([email protected])

Abstract

AC

Keywords Seismic tomography, volcano seismology, seismic attenuation, body waves, surface waves

In this paper we present innovative methodologies for seismic monitoring of volcanic structures in space and time (4D) which can possibly evolve toward an unrest stage. They are based on repeated phase and amplitude measurements done on active and/or passive seismic data including shots, vibrations, earthquakes and ambient noise in order to characterize the structure of the volcano and track its evolution through time. The characterization of the medium properties is performed through the reconstruction of an image of the elastic and anelastic properties of the propagation medium crossed by seismic waves. This study focuses on the application of specific tomographic inversion methods to obtain high quality tomographic images. The resolution of the tomographic models is influenced by the number and spatial distribution of data. The expected resolution thus guides the setup of, for example, active seismic surveys. To recognize and monitor changes in the properties of the propagation medium without performing an active survey we identify a fast proxy based on the time evolution of the Vp/Vs ratio. The advantages and limitations of the methods are discussed through synthetic tests, resolution analysis and case studies in volcanic areas such as the Campi Flegrei (southern Italy) and The Geysers geothermal area (California).

1

ACCEPTED MANUSCRIPT 1. Introduction Volcanic eruptions have always attracted interest from, media, and scientific community. Scientific studies have led to remarkable progress in the modeling of physical phenomena ruling volcano

T

dynamics, but much remains to be done to reach their complete understanding, and most of all,

IP

identify their precursory.

One of the most powerful geophysical prospecting methods is seismic tomography. It is an

SC R

inference technique, which exploits the information contained in the seismic records i.e., arrival times and amplitude of seismic phases to obtain 2D and 3D models of physical parameters such as seismic phase velocity or anelastic attenuation (Lees 2007). Tomography requires an inverse

NU

problem solution to match the model with the observations (Nolet 2008, Rawlinson et al 2010). From the pioneering work of Aki and Lee (1976), methodological advancements in this field allow

MA

for the imaging of tectonic and volcanic areas with resolution down to the space scale of the order of tens of meters (De Landro et al , 2017).

It is possible to reconstruct the subsurface structure using artificial seismic sources produced, for

D

example, by underground explosions (active seismic) or natural sources such as earthquakes or

TE

noise (passive seismic). In the first case, an array of seismic sensors is deployed at surface, close to the region to be studied. The acquisition layout (the relative receiver-source distribution) is

CE P

designed according to the desired resolution. In passive seismic, both the possible scattered deployment of sensors and the source location are unknowns for the inverse problem. The aim of this paper is to describe methods of seismic tomography used for imaging volcanic structures, with particular attention to inversion strategies. Several case studies of local tomography

AC

are presented to provide an overview of the different types of tomographic techniques that are carried out. In particular, we present the application of tomographic techniques to active seismic data sets relative to the volcanic area of Campi Flegrei (Southern Italy) and a passive seismic data set recorded at The Geysers geothermal area (California).

2. Seismic tomography methods for imaging volcanic structures Seismic tomography represents an excellent tool to image crustal heterogeneities (i.e. the position, extension and intensity of anomalies) using either shots and/or local earthquake data (e.g. Zollo et al 2002; Husen et al 2004). Very striking results for the P- and S-wave velocity models can be obtained by performing a simultaneous inversion of P- and S-wave arrival times to estimate the event location coordinates, their origin time, and P- and S-wave velocities at the nodes of a discretized 3D volume (Vanorio et al 2005) as long as a good approximation of the real medium is available as a starting model in the inversion. A useful strategy to minimize the influence of the 2

ACCEPTED MANUSCRIPT starting model is adopting a multi-scale approach (e.g. Amoroso et al 2014; De Landro et al 2017). While velocities provide information on the elastic properties of rock, the lateral variation of the anelastic parameter Q can be retrieved by the use of attenuation tomography (De Siena et al 2009;

T

Hauksson and Shearer 2006; Serlenga et al 2016). Possible time variations of seismic properties can

IP

be observed by repeating the tomography in time (i.e. 4D tomography - Patanè et al 2006; Koulakov et al 2013). Finally P- and S- wave models with different resolution can be constrained using

SC R

dispersion curves measured on surface waves either obtained by analyzing the signals generated by shots or by cross-correlation of ambient noise.

NU

2.1 3D velocity tomography

The most advanced tomographic inversions use an iterative scheme operating a linearized delay-

MA

time inversion to estimate both velocity models and earthquake locations (e.g. Latorre et al., 2004 and references therein). First arrival travel times of wavefronts are computed through a finitedifference solution of the eikonal equation (Podvin and Lecomte, 1991) in a finer grid of nodes,

D

consisting of constant slowness nodes computed by tri-linear interpolation from the inversion grid.

TE

For each source-receiver pair, travel times are recalculated by numerical integration of the slowness on the inversion grid along the rays traced in the finite-difference travel time field (Latorre et al.,

CE P

2004). Simultaneously, for each node of the inversion grid, travel time partial derivatives are computed for the P slowness, hypocenter location and origin time. The parameters are inverted using the LSQR method (Paige and Saunders, 1982). Model roughness is controlled requiring that the Laplacian of the slowness must vanish during the inversion procedure (Benz et al., 1996;

AC

Menke, 1989). The velocity model is parameterized by a nodal representation, described by a tridimensional grid.

2.2 3D attenuation tomography using t* and dt* If the velocity model is known and the t* parameter is estimated along the corresponding path, the 3-D attenuation quality factor Q of the medium along this path can be solved by an appropriate inversion method (Lees and Lindley, 1994; Haberland and Rietbrock, 2001; Rietbrock, 2001; Hansen et al., 2004; Martinez-Arevalo et al., 2005; Chiarabba et al., 2009; Bisrat et al., 2013; Amoroso et al., 2017). Since the velocity model is fixed, ray paths do not change in the inversion procedure and the problem is fully linear from a mathematical point of view. Nevertheless, the uncertainties on data, in addition to the mixed-determined nature of the problem (Menke, 1984) ask for the data to be inverted by means of a linearized, iterative, perturbative scheme requiring a starting attenuation model as well as some form of regularization of the inverse problem and constraints. In this model a theoretical calculation of t* is performed to evaluate the residue between 3

ACCEPTED MANUSCRIPT observed and expected t* values: Δt*=t*obs-t*cal. This condition defines the set of equations to be inverted, which are coupled with smoothing and preconditioning conditions to avoid instabilities to emerge in the solution. Finally the inversion of the linear system Δt*=H is performed using the

T

LSQR method, where the term H describes the data-kernel, is the perturbation of the reciprocal

IP

of the quality factor with respect to the starting attenuation model. Once the Q attenuation model has been obtained, the Root Mean Square (RMS) of residuals is evaluated, and the procedure is

SC R

reiterated until this value is smaller than a certain threshold.

Tomographic attenuation images may be obtained also by differential attenuation measurements (Serlenga et al., 2016), i.e. dt*. The problem can be faced using a similar approach to that of Teng

NU

(1968), where a thorough dissertation of the problem may be found. For a single source and a set of receivers (j is the reference station, i=1,N the remaining stations), the data vector is (1)

dti*  ti*  t *j .

(2)

MA

dt1*................dtn* ,

TE

D

where

Let us develop the problem only for the quantity dt*1 for sake of simplicity. With fixed velocity and

CE P

attenuation models, a residual quantity, dt*1, is obtained as

dt1*  (t1*OBS  t1*CALC )  (t*jOBS  t*jCALC )   t1*   t*j

(3)

AC

in which the term dt*1 may be also defined as the double-difference contribution (Waldhauser and Ellsworth, 2000). In the (3), each term in brackets, t*, represents the residual between observed and computed t* in fixed velocity and attenuation models. In a 3D tomographic grid, each t* residual can be expressed as

 ti*   cube

ti*

l ,n, p

l ,n, p .

(4)

In the above equation the subscript i represents the receiver index. The tern of indexes l,n,p, instead, refers to the nodes of the tomographic grid. Let us assume to have a single source attenuation tomographic problem, for which seismic signals have been recorded at N stations. Parameterizing the investigated medium in a 3D tomographic grid described by M nodes, constituting the set of parameters to be determined, the resulting system of double-difference equations could be written in a matrix form: 4

ACCEPTED MANUSCRIPT *m t1*m t j  1 m m dt1* ... ..... ... .....  ... ………… ..... ... ..... ... *1 *m dtn* tn*1 t j tn*m t j    m 1 1 m m *1 t1*1 t j  1 1

SC R

IP

T

(5)

The inversion scheme, as well as the adopted numerical algorithm, is analogous to the one

NU

described in the first part of this paragraph, with the exception that absolute t* values are substituted by differential values.

MA

2.3 4D imaging

2.3.1 Temporal variation of Vp/Vs ratio

D

The Vp/Vs ratio is a quantity directly correlated with the presence of fluids within the crust (Thurber

TE

et al. 1995; Amoroso et al 2017). The analysis of Vp/Vs ratio temporal variations is thus a technique for imaging the large-scale medium properties. It provides complementary information with respect

CE P

to a tomographic inversion procedure and it has the advantage that the results can be computed as soon as the seismogram is available. The Vp/Vs ratio value is computed from the time difference between S- and P-wave arrival times

AC

(indicated by tp and ts respectively) at a given station, divided by the P-wave travel-time (ttp) from the hypocenter to the station (Lucente at al., 2010; Wadati, 1933; Kisslinger and Engdahl, 1973; Chiarabba et al., 2009b) using the following formula

VP t -t =1+ S P VS ttP

(6)

The above relation holds for a uniform Vp/Vs ratio along the travel path (Scholz et al., 1973; Whitcomb et al., 1973). By considering the Vp/Vs ratio temporal variation for several events that occurred in the same source region, and recorded at different pairs of nearby stations, we are able to reveal possible changes in crustal volume of events occurrence (Figure 1). The basic idea is to analyze the graphs of VP/VS ratio as a function of time for each couple of stations, and to identify both spatial and temporal changes by comparing them to each other.

5

MA

NU

SC R

IP

T

ACCEPTED MANUSCRIPT

Figure 1: Schematic representation of the rays path for a given sources-stations geometry. The common

D

area of all rays is the source region (red dashed line).

TE

2.3.2 4D seismic tomography

Time-lapse tomography is a technique that consists in applying the three-dimensional tomography

CE P

in different time-windows, referred to as epochs. To properly define the epochs, a preliminary 3D reference velocity model and related resolution have to be obtained. These results represent the benchmark for the resolution analysis in the individual epochs. The resolution can be assessed

AC

through the ray coverage, the derivative weight sum (DWS - Toomey and Foulger, 1989), the checkerboard test, and the resolution matrix, represented by maps of its diagonal elements and the spreading function (Sj - Michelini and McEvilly, 1991). For each of these, it is possible to define a threshold value that can be used to delimit the model regions that are considered well-resolved. The threshold should be set up specifically for each tomography, for example implementing synthetic tests. The adaptive selection process of the temporal duration of an epoch consists in selecting specific time windows and using the data set recorded within this time interval to achieve a 3D tomography. The solution is evaluated and compared with the one obtained in the 3D reference model. If the resolution is too low with respect to the resolution of the benchmark, the length of the time window has to be increased including more data. In order to have consistent model resolution between the different epochs, the epochs may have different durations.

6

ACCEPTED MANUSCRIPT 2.4 Surface wave tomography from phase and group velocity analysis using dense array records The estimation of a 3D model for the investigated area is achieved collecting the 1D models coming

T

from the joint inversion of the dispersion curves in each selected subdomain, where the 1D

IP

approximation holds. The local 1D model can be obtained from a joint inversion of group and phase

SC R

velocities (e.g. using the Geopsy software, Wathelet et al. 2004). The dispersion curves depend on the number of layers, the density and the velocity values for the P- and S-waves in each layer. To properly set the number of layers, recursive test can be performed increasing the number of layers and evaluating the corresponding misfit reduction. An AIC criterion can be applied to define the

NU

optimal number of interfaces required to fit the dispersion curves. For a fixed number of layer the solution is sought using a global exploration method such at Neighborhood algorithm (Sambridge,

MA

1999), which divides the parameter space into Voronoi cells and refines the sampling in the regions nF



of the space where the cost function, defined as

i 0

 xdi  xci 2 i2nF

where xdi is the measured velocity

D

at a given frequency fi, xci is the computed velocity at a given frequency fi, σi is the uncertainty on

TE

the measurement of the considered frequency and nF is the number of frequency samples, is smaller

AC

2.5 Ambient noise

CE P

than elsewhere.

In contrast to global tomographic studies we are not seeking for small perturbations to an established reference model but we intend to locate velocity differences and corresponding velocity variations on a small local scale. Opposed to active seismic studies in volcanic environments (De Luca et al., 1997), several approaches have been proposed to extract information about the 3D velocity distribution of seismic waves in shallow structures from temporary recording of seismic noise (Chávez-García and Luzón, 2005, Brenguier et al., 2007, Picozzi et al., 2009) although the use of high frequency seismic noise interferometry implies by far not only a change in scale but also requires a sufficient understanding of the origin of seismic noise and of the spatial and temporal distribution of its sources. Generally, all techniques for retrieving surface wave dispersion curves are based on phase-coherency measurements between pairs (at least two) of signals. Aki (1957) proposed the spatial autocorrelation (SPAC) method, which has been generalized in the extended spatial autocorrelation (ESAC) method by Ohori et al. (2002). If the microtremor wavefield is stochastic and stationary in both time and space (Ohori et al., 2002), 7

ACCEPTED MANUSCRIPT the azimuthally averaged correlation function for a single angular frequency ω0 can be calculated by means of

 0  r  c(0 ) 

  r, 0   J0 

T

(7)

IP

Here, c(ωo) represents the phase velocity, r is the interstation distance and Jo is the zero order

SC R

Bessel function. Equation (7) can be applied to averaged correlation functions calculated for a set of narrow frequency bands.

In turn, the calculated travel times for the surface wave for individual frequencies were

NU

simultaneously inverted by a rapid one-step tomographic algorithm. In its general form, the travel time between a source and a receiver along a ray path element for a continuous slowness s (reciprocal of velocity) is written in an integral form. However, although the ray path is velocity

MA

dependent, meaning that travel time inversion is a non-linear problem, deviations of the paths from a straight line will either be of the same order as the dimension of the blocks size or smaller will be

D

less than a quarter of the wavelength. Considering the error level of the input data for shallow seismic surveys, we are sure that a bias of a few percent can be tolerated keeping the solution linear

TE

(Kugler et al., 2007, Pilz et al., 2012).In turn, the studied medium can be subdivided into smaller blocks, and the problem can be expressed in a simple discrete matrix form (Pilz et al., 2012, 2013).

CE P

Starting from a homogeneous 3D velocity model, we adopt an iterative procedure based on singular value decomposition for minimizing the misfit between the observed and theoretical travel times. In particular, as higher frequency sample shallow blocks whereas lower frequencies sample deeper

AC

blocks a further constraint on the solution can be introduced by adding a weighting matrix. According to Yanovskaya and Ditmar (1990), for a 3D problem, the weights can be calculated as the product of two functions with one weight depending on the horizontal properties, i.e. the number, the length and the orientation of the ray paths and the other one depending on the depth resolution, i.e. the frequency. A good resolution, meaning a large relative value for the horizontal weight, is achieved if many rays with different azimuths cross the cell. In a simplified way, the horizontal weights are computed by multiplying the ellipticity for the number of rays crossing each cell. The vertical weights account for the different penetration depths of the different frequencies of the surface waves. Consequently, the vertical weights are based on the analytical solution of displacement components for the fundamental mode in a half space (Aki and Richards, 1980). Additionally, for stabilizing the iteration process and for reducing the risk of divergence, an adaptive bi-weight estimation (Tukey, 1974; Arai and Tokimatsu, 2004) can be applied. The solution is constrained to smoothly vary in the horizontal domain, meaning that the slowness of each cell is related to the value of the slowness of all the surrounding cells. As the vertical weights 8

ACCEPTED MANUSCRIPT are based on a continuous functional form, this ensures that the slowness also smoothly varies vertically. For each cell, the individual weights and the corresponding surface wave velocity values are

T

updated after each iteration until a reasonable compromise between the reduction of the root-mean-

IP

square error between the observations and the predictions and the norm of the solution is reached. The final 3D velocity model is obtained from the inversion of the cross-correlations based on a

SC R

singular value decomposition technique (e.g., Arai and Tokimatsu, 2004; Pilz et al., 2012).

NU

3. Application on real datasets

3.1 Campi Flegrei caldera and the Solfatara Volcano

MA

Campi Flegrei is a volcanic polygenetic complex located along the western coast of Southern Italy, 15 km west of Naples (De Natale et al., 1987). The Campi Flegrei volcanic field is a partly submerged nested calderic system whose current shape is directly linked to its eruptive history. The

D

greatest events occurred about 39,000 (Campanian Ignimbrite) and 15,000 years ago (Neapolitan

TE

Yellow Tuff) (Deino et al., 2004; De Vivo et al., 2001). The current activity of this volcanic system consists mainly in widespread fumaroles inside and at the border of caldera, thermal spring

CE P

activities as well as ground deformation episodes. In this regard, the subsidence trend of the last five centuries has been sometimes interrupted by unrest episodes characterized by short-term ground uplift. Actually, two large bradyseismic crises occurred in the last 40 years: the former between 1969 and 1972, the latter between 1982 and 1984 (Beaducel et al., 2004; Battaglia et al., 2006; Del

AC

Gaudio et al., 2010) . Both of them were accompanied by a shallow seismicity in the range 1-4.2 and brought to a cumulative uplift of about 3.5 m in a 10-15 km wide circular area (del Pezzo et al., 1987, De Natale et al., 1995, de Lorenzo et al., 2001, De Natale et al., 2006, Zollo et al., 2008, Amoruso et al., 2014, Carlino et al., 2015). As a consequence of these episodes, the risk of an imminent eruption was perceived by the population living in the active portion of the supervolcano (about 350,000 people) (Di Renzo et al, 2011; Capuano et al., 2013). On these grounds, in recent years the interest of the scientific community greatly increased in order to gain insight into the interior of the volcano. To this purpose, in the recent past several subsurface images of the Campi Flegrei caldera were obtained by means of different seismic methodologies (de Lorenzo et al., 2001; Zollo et al., 2003; Judenherc and Zollo, 2004; Zollo et al., 2008; Battaglia et al., 2008; Dello Iacono et al., 2009, De Siena et al., 2010, Serlenga et al., 2016) Among the most studied areas of this volcanic field, there is the Solfatara volcano, which is one of the craters of the Campi Flegrei caldera with a diameter of about 0.6 km. It is believed to be the top of an hydrothermal plume 9

ACCEPTED MANUSCRIPT (Chiodini et al. 2001), being the shallow hydrothermal system of the Campi Flegrei caldera fed by the interaction between the meteoric water and heat from depth. Solfatara often experiences episodes of degassing, with very intense emission of hydrothermal–magmatic gases and high heat

T

flow, ascending along faults bounding the crater (Chiodini et al. 2001).

IP

The shallow part of this volcanic structure was investigated through the seismic campaign RICEN (Repeated InduCed Earthquake and Noise), recording both active seismic data and ambient noise.

SC R

Active seismic data were recorded using a single 6400 kg IVI-MINIVIB® vibroseis truck, which delivers a maximum theoretical peak force of ~27 kN at each sweep. For the two 2D seismic profiles a nominal source spacing of 4 m was used, with vibration points located halfway between

NU

geophones (Figure 2, left panel). At each vibration point, two, 15s long, 5-150 Hz sweeps were stacked and recorded respectively by a 216-channel and 240-channel, with spacing between the

MA

individual sensors of 2 m. The source-receiver configuration allowed sampling a wide range of offsets, from a minimum of only 0.5 m to a maximum of 453 m. 3D data were acquired on a regular grid, composed of a 240-channel array. For the regular 3D array we used 10 parallel sub-arrays of

D

24 channels each, oriented roughly with a NE strike. Inline spacing of geophones was 5 m, while

TE

cross-line spacing was 10 m. This yielded a 120x100 m2 wide grid of geophones (Figure 2, right panel). Vibration points are located along parallel lines located midway between two parallel

CE P

geophone lines. Nominal source spacing was 20 m, within the rectangular grid of four adjacent geophones. Ambient noise data were continuously acquired for three days at 50 stations on a larger

AC

area, spanning the entire crater.

Figure 2. Acquisition layout of seismic campaign RICEN. Left: location of the two linear arrays used for the 2D seismic reflection-refraction experiments. Data were collected in May 2014 along the NNE striking array and in November 2014 along the WNW trending array. Right : aerial view of the Solfatara crater with location of the recording array. Specifically the 2D array for the 3D tomography is located at the NE boundary of the mud region, named “Fangaia”.

10

ACCEPTED MANUSCRIPT A data set with more than 100,000 seismograms was collected and preliminary used to infer a 1D model of the elastic and anelastic properties of Solfatara. Starting from an initial dataset of 21,315 automatically determined P-wave first arrival time picks, 17,847 data have been manually selected

T

to be used to obtain a 3D P-wave tomographic image.

IP

In order to retrieve a P wave attenuation model, a further manual selection of the dataset has been performed by discarding signals with a low signal-to-noise ratio. To remove the source effect from

SC R

the signal, the velocity seismograms have been cross-correlated with the sweep (Brittle et al., 2001). The cross-correlated traces have been cut in a window of 0.128 s around the P-wave arrival time and zero-padded up to a total duration of 0.256 s. Then, the velocity signals have been integrated in

NU

the time domain and the displacement spectra have been computed. The natural logarithm of the spectral amplitudes as a function of the frequency has been fit with a straight line in the frequency

3.1.1 Resolution test and ray coverage

MA

range 40-125 Hz in order to obtain for each source-receiver couple the t* value.

D

In order to evaluate the spatial resolution of a 3D model for the source-receiver geometry of the

TE

RICEN experiment at Solfatara volcano, we performed a numerical computation of the resolution matrix through spike tests. For this analysis, a 1D velocity model for the area was inferred from the

CE P

previous study of Bruno et al. (2007) up to 70 m depth. The velocity model is specified at the nodes of a cubic lattice of 10 m spacing. We added a velocity perturbation of 400 m/s to each grid node and computed the synthetic travel times that were used as input data for the tomographic inversion. In figure 3 we show the diagonal elements of the resolution matrix (RDE) overlapped to the map of

AC

Solfatara Volcano, and the station/source distribution. The resolution is closely related to the source-receiver geometry. In fact, the resolution is maximum in the area cover ed by the array. Figure 4 shows both RDE and the spread function Sj, which gives a measure of the order of magnitude of the off diagonal elements of the resolution matrix. For this analysis we did not consider data errors and we assumed that all stations recorded all events, i.e., we considered a total amount of 25,920 P-arrival times with respect to the 21,315 preliminary picks.

11

CE P

TE

D

MA

NU

SC R

IP

T

ACCEPTED MANUSCRIPT

Figure 3: Spatial resolution at surface, superimposed to the map of Solfatara and the acquisition geometry

AC

of the experiment .The black box delimits the same area of Figure 4.

12

AC

CE P

TE

D

MA

NU

SC R

IP

T

ACCEPTED MANUSCRIPT

Figure 4: Model resolution estimation. a) Map view showing RDE at several depth ranging between 0 and 70 meters.b) Map view showing the Sj spread function with depth ranging between 0 and 70 m.

According to this analysis, the model has a good resolution up to 40 m depth. In fact, RDE is larger than 0.7 down to 40 m depth in the central part of the investigated area, whereas, at larger depth, the number of well resolved nodes drastically drops down. The spread function has a pattern very 13

ACCEPTED MANUSCRIPT similar to that of the RDE. This behavior is indicative of the fact that the off-diagonal elements are small and there is little correlation between parameters. As a confirmation of the results, in Figure 5

MA

NU

SC R

IP

T

the plot of the ray coverage in the vertical and horizontal planes are shown.

TE

D

Figure 5: Ray coverage. a) Ray distribution in the XY plane, b) ray distribution in the XZ,plane, c) 1D velocity model used for ray tracing.

3.1.2 Synthetic test of 3D attenuation tomography of Pozzuoli bay area using dt* data

CE P

A three-dimensional Qp tomography of offshore part of Campi Flegrei caldera was already performed in the very recent past by inverting differential t* measurements, with the method described in section 2.2 (Serlenga et al., 2016), Nevertheless, a synthetic test is here described to

AC

further show the reliability of a three-dimensional attenuation tomography based on the inversion of dt* measurements. The adopted source-receiver configuration allowed investigating a tomographic volume of 13x13x5 km3. It is the same layout used for the active seismic experiment SERAPIS, in September 2001, in the Pozzuoli Bay, southern Italy (Judenherc and Zollo, 2004). The node spacing of the tomographic grid in the described synthetic test is 500x500x250 m3. Seismic rays were traced in a 1D velocity model obtained by averaging the depth-dependent P wave velocity model of Battaglia et al. (2008). Synthetic attenuation model is characterized by a low-Qp annular anomaly (Qp=100) located at depths between 500 m and 1.5 km in a homogeneous attenuation medium (Qp=500). In the tomographic volume 14,450 synthetic dt* are computed (Figure 6). The tomographic inversion of synthetic data is run using a homogeneous Qp=300 structure as starting model. The inversion is stopped at the 8th iteration and provides the tomographic images as shown in figure 7. The attenuation anomaly is fully recovered from a geometrical point of view. Nevertheless, due to the poor ray coverage, high relative errors (about 40 %) affect the volume around the annular anomaly at depth of 1.25 km and 1.5 km. 14

SC R

IP

T

ACCEPTED MANUSCRIPT

AC

CE P

TE

D

MA

NU

Figure 6: 3D attenuation tomography synthetic test. a) Attenuation model used in the synthetic test. The triangles represent the receivers, whereas the black points identify the shots. b) 1D P wave velocity model used for ray-tracing. c) Synthetic dt* data used for the inversion.

Figure 7:Results of 3D attenuation tomography synthetic test. a) Tomographic images at different depths retrieved after 8 iterations. Grey regions indicate areas not covered by rays. b) relative error of the tomographic results. c) RMS as a function of the number of iterations. A stable value is achieved after 8 iterations. d) Residuals after the first iteration (blue histogram) and after the 8th iteration (green histogram).

3.1.4 Surface wave tomography Phase and group velocity estimated on the array of Figure 2, right panel, are inverted to obtain a 3D model of the area, from combination of local 1D models. In each sub-array, the inversion of phase and group velocity using a global exploration technique (Wathelet et al. 2004) furnishes more models that almost equally fits the measured dispersion curves. Therefore, it is more reasonable to infer statistical properties from the family of solutions, having the similar misfits to the 15

ACCEPTED MANUSCRIPT experimental dispersion curves. At each depth, the final model was selected as the one resulting from the average of all the models with misfit less than 15% with respect to the minimum one (Figure 8). This choice will produce a final non-layered model, which presents a smooth variation

T

with depth. In this analysis only the S wave model is constrained by the dispersion curves, as these

IP

curves are poorly affected by variations of both P-velocity and density. The maximum penetration depth is about 12 m (Figure 8 right panel)

SC R

The 3D S-wave velocity model obtained from the collection of all the 1D models is represented in Figure 9, in slices at fixed depths. The most remarkable trend is the presence of low-velocity anomalies due to the water that permeates the westernmost area, outcropping to the surface in the

AC

CE P

TE

D

MA

NU

Fangaia.

Figure 8: Results of the inversion obtained using the software Geopsy. On the left and central panels, the retrieved dispersion curves for the group and phase velocities respectively, superimposed with the measured ones. On the right panel, the S-wave models explored by the inversion procedure and the average model(dotted black line). The colors refer to the misfit value .

16

NU

SC R

IP

T

ACCEPTED MANUSCRIPT

MA

Figure 9: Slices of the 3D model at different depths. Results beyond 12 m are not shown because of the lack of resolution.

D

3.1.5 Ambient noise tomography results

TE

Using the vertical cross-correlations and a fraction of two hours of the entire noise data set, the inversion results after 500 iterations are shown in Figure 10. At first glance one can notice the

CE P

occurrence of strong lateral velocity variations. The topmost layer just below the surface is characterized by S-wave velocities ranging from 120 m/s for the north-eastern part of the Solfatara crater close to the vegetation cover and associated background values of ground temperature and gas flux up to velocities of 200 m/s for the south-western part in agreement with the results of Serra

AC

et al. (2016). However, the penetration depth of their study was limited to the uppermost 10 to 15 m. A significant increase of the S-wave velocity starts immediately below the shallow deposits in the southern and central part at depths of around 10 and 20 m (Figure 10), probably at the lower boundary of the unconsolidated shallow deposits. The lower velocity area remains confined in the northern and south-eastern part of the crater. At a depth of around 35 m, the S-wave velocity in the southern area increases up to 850 m/s. At these depths, the volume is divided in a higher velocity part (corresponding to the fumaroles area) and a lower one (corresponding to the northern edge), with a transition central zone. High S-wave velocities and corresponding resistivity and gas flux values emphasize a highly resistive body indicating that resistivity changes at this depth are related to high temperature gas saturation (Byrdina et al., 2014, Pilz et al., 2017). This large velocity contrast and the measured velocity values might assure the hypothesis of the presence of a small volcanic edifice (most probably a tuff-cone) formed in a smaller crater in the eastern sector during the phreatomagmatic activity (Petrosino et al. 2012). Below 40 to 45 m, the S-wave velocity 17

ACCEPTED MANUSCRIPT increases up to 1200 m/s in the south-western area. The absolute depth of the velocity contrasts and the corresponding estimated S-wave velocities are robust and compatible with the findings of

MA

NU

SC R

IP

T

previous studies (Petrosino et al. 2006, 2012).

AC

CE P

TE

D

Figure 10: Horizontal cross sections of the Solfatara volcano in terms of S-wave velocity after inverting two hours of seismic noise.

18

ACCEPTED MANUSCRIPT 3.2 The Geysers geothermal area The Geysers geothermal field is located 120 km north of San Francisco. It has been actively exploited since the 1960's and is now the most productive geothermal area in the world. The

T

geothermal field is monitored by a dense seismic network maintained by the Lawrence Berkeley

IP

National Laboratory Calpine (BG) and by the Northern California Seismic Network (NCSN). The BG network consists of 32 three-component stations, 29 of which were used for the present study

SC R

(black triangles in the Figure 11).

The analysis is here focused on natural and induced seismicity recorded from August 2007 to October 2011. The whole dataset contains about 15,000 events with magnitude ranging between 1.0

AC

CE P

TE

D

MA

NU

and 4.5. We selected only 1320 events for which at least 20 high-quality P-wave picks are available.

Figure 11. Distribution of induced seismicity recorded at The Geysers. Black triangles indicate the seismic stations of the Lawrence Berkeley National Laboratory (LBNL) Geysers/Calpine seismic network used in this study, and gray triangles are additional stations from the Northern California Seismic Network (NCSN) in the region (modified after Convertito et al., 2012).

Several authors have noted a difference in the seismicity features along the north-west -south-east direction (Convertito et al., 2012; Picozzi et al. 2017). In particular, Beall and Wright (2010) identified a clear “M ≥ 4.0 dividing line”, which splits the whole field into two different seismic areas. The north-western area (ZONE1) contains all the epicenters of the earthquakes with 19

ACCEPTED MANUSCRIPT magnitude larger than 4.0, whereas the southeastern one (ZONE2) is characterized by lower magnitude earthquakes. In ZONE1 the normal steam dominated reservoir (at temperature ~240°C) is underlined by a high temperature steam dominated reservoir (HTR) at 260-360 °C (see figure 14,

T

right part). The top of the shallower steam reservoir is located at about 1 km depth. The top of

IP

felsite (granitic intrusion), which underlines much of the steam reservoir, is shallower in ZONE2

SC R

compared to ZONE1. 3.2.1 3D body-wave seismic tomography

The 3D P- and S-wave velocity models are obtained from the tomographic inversion of the first Pand S-wave arrival times. The selected stations and events distribution allowed us to investigate a

NU

volume of 36×25×5 km3. The velocity model is parametrized by trilinear interpolation on a tridimensional grid with node spacing of 1×1×0.5 km3. The inversion starts from the 1D velocity

MA

model, optimized for the area, which is also used for the initial earthquake locations (Emolo et al., 2012). The misfit function, defined as the sum of the squared time delays (RMS), is a posteriori analyzed to check the convergence, which is reached after 8 iterations, where the final RMS value is

D

equal to 0.1 s, with an RMS reduction of 50%.

TE

The tomographic P- and S-wave velocity models shown in Figure 12 a,b respectively, seem to delineate the main geological features of the field. For example, the iso-velocity curves at 5.4 km/s

CE P

for Vp and 3.2-3.4 km/s for Vs (see Figure 13) correspond to the top of the felsite (Thompson 1989). Most of the seismicity in ZONE1 is clustered between 2.5 and 4.5 km depth, whereas the seismicity in ZONE2 is shallower.

AC

The Vp/Vs ratio throughout the reservoir reveals strong deviations from the expected value of 1.73. A high Vp/Vs anomaly is present in ZONE1 where the highest temperature in the field is measured (Beall and Wright 2010) and the most of seismicity occurs (Figure 13).

20

CE P

TE

D

MA

NU

SC R

IP

T

ACCEPTED MANUSCRIPT

AC

Figure 12: P and S-wave velocity models. a), Map view of P-wave velocity model at 0, 1.5, 2, 3, 3.5 and 4 km depth from 3D travel time tomography. The seismic stations are imposed on the 0 km depth layers as open triangles, whereas black dots in all panels represent the earthquake locations. Regions not covered by ray-paths are in grey. The solid white contour lines identify the model regions where derivative weight sum (DWS) is greater than 5000. The images indicate the presence of a strong lateral variation of seismic velocity along the NW-SE direction. b) Map view of S-wave velocity model at 0, 1, 2, 3, 3.5 and 4 km depth. Again, the images indicate the presence of a strong lateral variation of seismic velocity along the NW-SE direction.

21

AC

CE P

TE

D

MA

NU

SC R

IP

T

ACCEPTED MANUSCRIPT

Figure 13:. Cross-section of P and S-wave velocity models and Vp/Vs ratio.The earthquakes are projected onto the NW-SE cross-section depicted by the black line in the first panel of figure 12.The blue solid line indicates the top of the steam reservoir, the orange line the top of felsite and the red line corresponds to the top of the high-temperature reservoir as reported by Beall and Wright (2010). a) Vertical cross-section along the NW-SE transect crossing the P-wave velocity model. b) Vertical cross-section along NW-SE ransect crossing the S-wave velocity model.c) Vp/Vs ratio images along the NW-SE transect.

3.2.2 4D velocity imaging Time variations of the seismic properties of the hosting medium can be observed by means of the variations of the Vp/Vs ratio evaluated at single stations and through time-lapse tomography (4D 22

ACCEPTED MANUSCRIPT tomography). The whole catalogue was divided in consecutive epochs of 6 months with an overlapping of 2 months for the analysis. A technique for the large scale medium properties identification was applied, which can provide

T

complementary information with respect to tomographic analysis The method allows to evaluate the

IP

temporal variations of the Vp/Vs ratio, using the arrival times of the P and S phases as seismological observables and then relate the ratio to the directional properties of the fluid diffusion process.

SC R

Taking into account separately the events belonging to the two zones (Fig. 15) and the catalogue subdivision defined above, the Vp/Vs ratio was evaluated for each event at each station. For each epoch an interpolated color map of the mean values of the ratio at each station is presented to show

AC

CE P

TE

D

MA

NU

its spatial variation.

Figure 14: Vp/Vsratio vs time epoch. The four panel show the Vp/Vs as a function of time epoch for different stations, grouped according to the similarity of their trends.

Vp/Vs ratios range from 1.65 to 1.85, in agreement with the work of Gunasekera et al. (2004). The assumption of this technique is that the Vp/Vs ratio is constant along the ray path, so that the value found at the stations may be considered the combination of the contributions of the anomalies at different depths. The plots in Figure 14 represent the 1D Vp/Vs ratio trends for some stations of the network that show similar patterns. Even these plots indicate a high variability of the Vp/Vs ratio in the different time intervals. For some stations it varies from values lower than 1.7 to 1.8, the former value being attributed to the depletion of pore liquid water and the replacement with steam, the latter being associated to the presence of liquid water (Gritto et al 2013).

23

AC

CE P

TE

D

MA

NU

SC R

IP

T

ACCEPTED MANUSCRIPT

Figure 15: Cross-sections of 4D seismic tomography results. Each panel shows the tomographic model, for each time epoch, projected onto the NW-SE cross-section plotted as a black line in the first panel of figure 19. The first and second columns represent the Vp and Vs the variation with respect to the 3D Vp and Vs model calculated on the entire period respectively (see Figure 13. ) The fifth column represents the Vp/Vs ratio model

The 4D tomography consists in applying the 3D tomography at consecutive epochs. The model inferred by considering the whole dataset was used as starting model in the inversion. The model parametrization is the same as for the initial model. We evaluated the results in terms of Vp and Vs, and their percentage variation with respect to the initial model (Figure 15). For each model the 24

ACCEPTED MANUSCRIPT DWS is computed to determine the well resolved regions. The 4D tomographic images show for all the epochs a variation that does not exceed 10%. The presence of some anomalies seems to be rather recurrent. This is the case for example of the + 5%

T

anomaly at the epochs from A and B and -3% anomaly at the epochs D, F, I and L. The observed

IP

changes in velocity might be correlated to the field operations, such as fluid injection or steam production, which are not constant during each year, but show seasonal variation (e.g. Convertito et

SC R

al., 2012). Finally, the inferred time variation in the velocity pattern (Vp,Vs, and Vp/Vs ratio) suggests

NU

a non- isotropic fluid diffusion in the whole geothermal field.

MA

4. Conclusions

In this paper we have presented four recent, integrated methods for high-resolution seismic tomographic imaging of volcanic structures. Methods differ in the input data type used and

D

inversion strategies. The first method is the body-wave velocity tomography performed with a

TE

linearized approach. We have also shown how to assess the resolution of the final model. The second method, attenuation tomography, uses t* absolute and dt* relative measurement, using

CE P

iterative inversion approaches. The last two methods proposed are the surface wave and shear wave tomography obtained from active shots and environmental noise. The advantages, limitations, and domains of applicability of these methods have been shown through practical examples applied to data provided by both seismic exploration campaigns and,

AC

synthetic passive seismic data. Acknowledgements

This paper was carried out within the framework of the MedSuv project, which received funding from the European Union Seventh Programme for research, technological development and demonstration, under grant agreement no. 308665. The authors wish to thank the editor and the reviewers, including L. De Siena, for their valuable comments, which have contributed to improve the quality of the manuscript. A Post-Doctoral Fellowship at the University of Naples, Federico II, supported this work.

Reference list Allen, R. (1978). Automatic phase pickers: their present and furure prospects. Bull. Seism. Soc. Am. 72, S225–S242.

25

ACCEPTED MANUSCRIPT Amoroso, O., A. Ascione, S. Mazzoli. J. Virieux and A. Zollo (2014), Seismic imaging of a fluid storage in the actively extending Apennine mountain belt, southern Italy, Geophys. Res. Lett., 41, 3802-3809, doi: 10.1002/2014GL060070.

T

Amoroso, O., G. Russo, G. De Landro, A. Zollo, S. Garambois, S. Mazzoli, M. Parente, and J.

IP

Virieux (2017), From Velocity and Attenuation Tomography to Rock Physical Modeling:

Geophys. Res. Lett., 44, doi:10.1002/2016GL072346.

SC R

Inferences on fluid-driven earthquake processes at the Irpinia fault system in Southern Italy,

Amoruso, A., Crescentini, L., Sabbetta, I., 2014. Paired deformation sources of the Campi Flegrei

119, 858-879, doi: 10.1002/2013JB010392.

NU

caldera (Italy) required by recent (1980-2010) deformation history, J. Geophys. Res. Solid Earth,

MA

Aki, K. (1957). Space and time spectra of stationary stochastic waves with special reference to microtremors. Bull. Earthquake Res. Inst. Univ. Tokyo,35, 415–456.

D

Aki, K. and W. Lee (1990). Determination of three-dimensional velocity anomalies under a seismic

Geophys. Res., 4381–99.

TE

array using first p arrival times from local earthquakes, a homogeneous initial model. Journal

Francisco, USA.

CE P

Aki, K., P. Richards (1980). Quantitative seismic theory and methods, W.H. Freeman, San

Arai, H., K. Tokimatsu (2004). S-wave velocity profiling by inversion of microtremor H/V

AC

spectrum. Bull. Seismol. Soc. Am., 94, 53-63. Battaglia, M., Troise, C., Obrizzo, F., Pingue, F., De Natale, G., 2006. Evidence for fluid migration as the source of deformation at Campi Flegrei caldera (Italy). Geophys. Res. Lett., 33, L01307, doi: 10.1029/2005GL024904. Battaglia, J, Zollo, A., Virieux, J., and Dello Iacono, D., 2008. Merging active and passive data sets in travel time tomography: The case study of Campi Flegrei caldera (southern Italy). Geophysical prospecting, 56, 555-573, doi: 10.1111/j.1365-2478.2007.00687.x. Beaducel, F. and De Natale, G., 2004. 3D modelling of Campi Flegrei ground deformations: role of Caldera boundary discontinuities. Pure Appl. Geophys., 161, doi:10.1007/s00024-004-2507-4.

26

ACCEPTED MANUSCRIPT Beall J. J., and M. C. Wright (2010). Southern extent of The Geysers high temperature reservoir based on seismic and geochemical evidence. In: GRC Transactions, Vol. 34. Benz, H., B. Chouet, P. Dawson, J. Lahr, R. Page, and J. Hole (1996). Three dimensional p and s wave velocity structure

T

of redoubt volcano, alaska. J. Geophys. Res. 101 (B4), 8111–8128.

IP

Bisrat., S. T., DeShon, H. R., Pesicek, J., Thurber, C., 2014. High-resolution 3-D P wave

SC R

attenuation structure of the New Madrid Seismic Zone using local earthquake tomography. J. Geophys. Res. Solid Earth, 119, doi:10.1002/2013JB010555.

Boschetti, F., Dentith, M.C. & List, R.D., 1996. Inversion of seismic refraction data using genetic

NU

algorithms, Geophysics, 61, 1715–1727.

MA

Boué, P., P. Roux, M. Campillo and X. Briand (2014). Phase velocity tomography of surface waves using ambient noise cross correlation and array processing. J. Geophys. Res., Solid Earth, 119(1), 519-529. Brenguier, F., N. M. Shapiro, M. Campillo, A. Nercessian, V. Ferrazzini (2007), 3D surface wave

D

tomography of the Piton de la Fournaise Volcano using seismic noise correlations, J. Geophys. Res.

TE

34, L02305, doi: 10.1029/2006GL028586

CE P

Byrdina, S., J. Vandemeulebrouck, C. Cardellini, A. Legaz, C. Camerlynck, G. Chiodini, T. Lebourg, M. Gresse, P. Bascou, G. motos, A. Carrier, S. Caliro (2014). Relations between electrical resistivity, carbon dioxide flux, and self-potential in the shallow hydrothermal system of Solfatara

AC

(Phlegrean Fields, Italy). Journal of Volcanology and Geothermal Research, 283, 172-182.

Brittle, K.F., Lines, L.R., and Dey, A.K. ,2001. Vibroseis deconvolution: a comparison of crosscorrelation and frequency-domain sweep deconvolution, Geophys. Prosp., 49, 675-686, doi: 10.1046/j.1365-2478.2001.00291.x Bruno, P.P.G., Ricciardi, G.P., Petrillo, Z., Di Fiore, V., Troiano, A. and Chiodini, G., 2007, Geophysical and hydrogeological experiments from a shallow hydrothermal system at Solfatara Volcano, Campi Flegrei, Italy: response to caldera unrest, J. Geophys. Res., 112, B06201, doi:10.1029/2006JB004383. Capuano, P., G. Russo, L. Civetta, G. Orsi, M. D’Antonio, and R. Moretti, R. (2013) The active portion of the Campi Flegrei caldera structure imaged by 3-D inversion of gravity data. Geochemistry, Geophysics, Geosystems, 14, 10, 4681-4697, doi: 10.1002/ggge.20276.

27

ACCEPTED MANUSCRIPT Carlino, S., Kilburn, C. R. J., Tramelli, A., Troise, C., Somma, R., De Natale, G., 2015. Tectonic stress and renewed uplift at Campi Flegrei caldera, southern Italy: New insights from caldera drilling. Earth and Planet. Science Lett., 420, 23-29, http://dx.doi.org/10.1016/j.epsl.2015.03.035.

T

Chávez-García, F. J., F. Luzón (2005), On the correlation of seismic microtremors, J. Geophys. Res.

IP

110, B11313, doi: 10.1029/2005JB003671

SC R

Chiao, L.-Y., and B.-Y. Kuo (2001). Multiscale seismic tomography, Geophys. J. Int., 145, 517527.

Chiarabba, C., D. Piccinini, and P. D. Gori (2009). Velocity and attenuation tomography of the

NU

umbria marche 1997 fault system: Evidence of a fluid-governed seismic sequence. Tectonophysics

MA

476 (1-2), 73–84.

Chiodini,G.,Frondini,F., Cardellini,C.,Granieri,D.,Marini,L.&Ventura, G., 2001.CO2 degassing and energy

releaseatSolfataravolcano,Campi

Italy,

J.

geophys.

Res.,

106,

D

doi:10.1029/2001JB000246

Flegrei,

TE

Cichowicz, A. (1993). An automatic S-phase picker, Bull. Seismol. Soc. Am. 83, 180–189.

CE P

Convertito V., N. Maercklin, N. Sharma and A. Zolllo (2012). From induced seismicity to direct tiime dependent seismic hazard . Bull Seismc Soc. Am. 102, 2563-2573, doi:10.1785/0120120036 Crampin, S. (1977). A review of the effects of anisotropic layering on the propagation of seismic

AC

waves, Geophys J. Roy. Astron. Soc. 49, 9–27, doi: 10.1111/j.1365-246X.1977.tb03698.x De Landro G., V. Serlenga, G. Russo, O. Amoroso, G. Festa, P. P. Bruno, M. Gresse, J. Vandemeulebrouck and A. Zollo (2017), 3D ultra-high resolution seismic imaging of shallow Solfatara crater in Campi Flegrei(Italy): New insights on deep hydrothermal fluid circulation processes, Accepted for publication on Scientific Reports. de Lorenzo, S., A. Zollo, A., and F. Mongelli (2001). Source parameters and three-dimensional attenuation structure from the inversion of microearthquake pulse width data: QP imaging and inferences on the thermal state of Campi Flegrei caldera (Southern Italy). J. Geophys. Res., 106, B8, 16,265-16,286, doi: 10.1029/2000JB900462. De Luca, G., R. Scarpa, E. Del Pezzo and M. Simini (1997). Shallow structure of Mt. Vesuvius volcano, Italy, from seismic array analysis. Geophysical Research Letters, 24(4), 481-484.

28

ACCEPTED MANUSCRIPT Del Gaudio, C., Aquino, I., Ricciardi, G.P., Ricco, C., Scandone, R., 2010. Unrest episodes at Campi Flegrei: A reconstruction of vertical ground movements during 1905-2009. J. Volc. and Geoth. Res., 195, 48-56, doi:10.1016/j.jvolgeores.2010.05.014.

T

Del Pezzo, E., G. De Natale, M. Martini, and A. Zollo, (1987). Source parameters

of

IP

microearthquakes at Phlegrean fields Flegrei (southern Italy) volcanic areas. Phys. Earth Planet.

SC R

Inter., 47, 25-42, doi: 10.1016/0031-9201(87)90064-1.

Dello Iacono, D., Zollo A., Vassallo, M., Vanorio, T. and Judenherc, S., 2009. Seismic images and rock properties of the very shallow structure of Campi Flegrei caldera (southern Italy). Bull. Volc.,

NU

71, 275-284, doi: 10.1007/s00445-008-0222-1.

De Natale, G., Zollo, A., Ferraro, A., Virieux, J., 1995. Accurate fault mechanism determinations

MA

for a 1984 earthquake swarm at Campi Flegrei caldera (Italy) during an unrest episode: implications for volcanological research. J. Geophys. Res.100 (B12), 24167–24185.

D

De Natale, G., G. Iannaccone, M. Martini, and A. Zollo (1987). Seismic sources and attenuation

doi:10.1007/BF00879360

TE

properties at the Campi Flegrei Volcanic area, Pure and Appl. Geophys., 125, 6, 883-917, doi:

CE P

De Natale, G., Troise, C., Pingue, F., Matrolorenzo, G., Pappalardo, L., Battaglia, M., Boschi, E., 2006. The Campi Flegrei caldera: unrest mechanisms and hazards. In: Troise, C., De Natale, G., Kilburn, C.R.J. (Eds.), Mechanisms of Activity and Unrest at Large Calderas. In: Spec. Pub.

AC

Geological Society London, vol.269, pp.25–45. De Siena, L., Del Pezzo, E. and Bianco, F., 2010. Seismic attenuation imaging of Campi Flegrei: Evidence of gas reservoirs, hydrotermal basins and feeding systems. J. Geophys. Res., 115, B09312, doi:10.1029/2009JB006938. De Siena, L., Del Pezzo, E., Bianco, F., Tramelli, A., (2009). Multiple resolution seismic attenuation imaging at Mt. Vesuvius, In Physics of the Earth and Planetary Interiors, Volume 173, Issues 1–2, Pages 17-32, ISSN 0031-9201, https://doi.org/10.1016/j.pepi.2008.10.015. De Vivo, B., G. Rolandi, P.B. Gans, A. Calvert, W.A. Bohrson, F.J. Spera, and H.E. Belkin (2001). New constraints on the pyroclastic eruptive history of the Campanian volcanic Plain (Italy). Mineral Petrol 73: 47-65.

29

ACCEPTED MANUSCRIPT Deino A.L., G. Orsi, S. de Vita, and M.Piochi (2004). The age of the Neapolitan Yellow Tuff caldera – forming eruption (Campi Flegrei Caldera-Italy) assessed by 40Ar/39Ar datig method, J. volc. and Geoth. Res., 133, 157-170, doi: 10.1016/S0377-0273(03)00396-2.

T

Di Renzo, V., I. Arienzo, L. Civetta, M. D’Antonio, S. Tonarini, M.A. Di Vito, and G. Orsi (2011).

IP

The magmatic feeding system of the Campi Flegrei caldera: Architecture and temporal evolution.

SC R

Chemical geology, 281, 227-241, doi:10.1016/j.chemgeo.2010.12.010

Emolo A., Maercklin N., Matrullo E., Orefice A., Amoroso O., Convertito V., Sharma N., Zollo A. “Analysis of induced seismicity at The Geysers geothermal field, California” AGU Fall Meeting,

NU

San Francisco, USA, 3 - 7 December, 2012.

Goldberg, X., 1989. Genetic Algorithm in Search, Optimization and Machine Learning, Addison-

MA

Wesley Pub Co., Boston.

Gritto, R., Seung Hoon, YOO, Steve P. JAR (2013). Three Dimensional Seismic Tomography at

D

The Geysers Geothermal Field, CA, USA, Thirti-Eighth Workshop on Geothermal Reservoir

TE

Engineering, Stanford University, California, February 11-13, SGP-TR-198. Gunasekera , R. C., G. R. Fougler and B.R. Julian (2003). Reservoir depletion at The Geysers

CE P

geohermal area, California, shown by four-dimensional seismic tomography. J.Geophys Res., 108(B3), 2134, doi:10.1029/2001JB000638. Haberland, C., and Rietbrock, A., 2001. Attenuation tomography in the western central Andes: a

AC

detailed insight into the structure of a magmatic arc. J. Geophys. Res., 106, B6, 11151-11167. Hansen S., Thurber, M., Mandernach, F., Haslinger, F., Doran, C., 2004. Seismic velocity and attenuation structure of the East rift Zone and south flank of Kilauea Volcano, Hawaii. Bull. Seism. Soc. Am., 94, 4, 1430-1440. Judenherc, S., and A. Zollo (2004). The Bay of Naples (southern Italy):Constraints on the volcanic structures

inferred

from

a

dense

seismicsurvey,

J.

Geophys.

Res.,

109,

B10312,

doi:10.1029/2003JB002876 Kisslinger, C., and Engdahl, E.R., 1973, The interpretation of the Wadati diagram with relaxed assumptions: Seismological Society of America Bulletin, v. 63, p. 1723–1736. Kugler, S., T. Bohlen, T. Forbriger, S. Bussat, G. Klein (2007): Scholte-wave tomography for shallow-water marine sediments, Geophys. J. Int. 168, 551-570

30

ACCEPTED MANUSCRIPT

Koulakov, I.

Evgeniy I. Gordeev, Nikolay L. Dobretsov, Valery A. Vernikovsky, Sergey

Senyukov, Andrey Jakovlev, Kayrly Jaxybulatov, Rapid changes in magma storage beneath the

T

Klyuchevskoy group of volcanoes inferred from time-dependent seismic tomography, In Journal of

IP

Volcanology and Geothermal Research, Volume 263, 2013, Pages 75-91, ISSN 0377-0273,

SC R

https://doi.org/10.1016/j.jvolgeores.2012.10.014.

Latorre, D., J. Virieux, T. Monfret, V. Monteiller, T. Vanorio, J.-L. Got, and H. Lyon-Caen (2004), A new seismic tomography of Aigion area (Gulf of Corinth, Greece) from the 1991 data set,

NU

Geophys. J. Int., 159, 1013-1031.

Lees, J. M., and Lindley, G. T., 1994. Three-dimensional attenuation tomography at Loma Prieta:

MA

Inversion of t* for Q. J. Geophys. Res., 99, B4, 6843-6863.

D

Lucente F.P., De Gori P., Margheriti L., Piccinini D., Di Bona M., Chiarabba C., Agostinetti N.P.,

TE

2010. Temporal variation of seismic velocity and anisotropy before the 2009 M-W 6.3 L'Aquila

CE P

earthquake, Italy. Geology 38, 1015–1018 Martinez-Arevalo, C., Patanè, D., Rietbrock, A., Ibanez, J. M., 2005. The intrusive process leading to the Mt. Etna 2001 flank eruption: Constraints from 3-D attenuation tomography. Geophys. Res.

AC

Lett., 32, L21309, doi:10.1029/2005GL023736. Matheney, M.,P., and Novack, R. L., 1995. Seismic attenuation values obtained from instantaneousfrequency matching and spectral ratios. Geophysical Journal International, 123, 1-15, doi: 10.1111/j.1365-246X.1995.tb06658.x. Menke, W. (1989). Geophysical Data Analysis: Discrete Inverse Theory. Academic Press. Michelini, A., & McEvilly, T. V. (1991). Seismological studies at Parkfield. I. Simultaneous inversion for velocity structure and hypocenters using cubic B-splines parameterization. Bulletin of the Seismological Society of America, 81(2), 524-552. Nolet G. (2008). A breviary of Seismic Tomography, Imaging the Interior of the Earth and Sun, Cambridge University Press.

31

ACCEPTED MANUSCRIPT Ohori, M., A. Nobata, K. Wakamatsu (2002). A comparison of ESAC and FK methods of estimating phase velocity using arbitrarily shaped microtremor arrays. Bull. Seismol. Soc. Am., 92, 2323-2332.

IP

squares. ACM Transactions on Mathematical Software 8 (1), 43–71.

T

Paige, C.,C., Saunders, M.A., 1982. Lsqr: An algorithm fo sparse linear equation and sparse least

SC R

D. Patanè, D., G. Barberi, O. Cocina, P. De Gori C. Chiarabba (2006). Time-Resolved Seismic Tomography Detects Magma Intrusions at Mount Etna, Science , Vol. 313, Issue 5788, pp. 821-

NU

823 DOI: 10.1126/science.1127724

MA

Petrosino, S., P. Cusano, G. Saccorotti (2006). Shallow shear-wave velocity structure of Solfatara volcano (CampiFlegrei, Italy), from inversion of Rayleigh-wave dispersion curves. Bollettino di Geofisica Teorica ed Applicata, 47, 89-103.

D

Petrosino, S., N. Damiano, P. Cusano, M. A. Di Vito, S. Vita, E. Del Pezzo (2012). Subsurface

TE

structure of the Solfatara volcano (Campi Flegrei caldera, Italy) as deduced from joint seismic‐

CE P

noise array, volcanological and morphostructural analysis. Geochem.Geophys.Geosys., 13, doi: 10.1029/2011GC004030. Picozzi M., A. Oth,

S. Parolai, D. Bindi, G. De Landro and O. Amoroso (2017). Accurate

AC

estimation of seismic source parameters of induced seismicity by a combined approach of generalized inversion and genetic algorithm: application to The Geysers geothermal area, California" Accepted for publication in Journal of Geophysical Research - Solid Earth Picozzi, M., S. Parolai, D. Bindi, A. Strollo (2009): Characterization of shallow geology by highfrequency seismic noise tomography, Geophys. J. Int., 176, 164-174

Pilz, M., S. Parolai, M. Picozzi, D. Bindi (2012).Three-dimensional shear wave velocity imaging by ambient seismic noise tomography. Geophys. J. Int., 189, 501-512. Pilz, M., S. Parolai, D. Bindi (2013). Three-dimensional passive imaging of complex seismic fault systems: evidence of surface traces of the Issyk-Ata fault (Kyrgyzstan). Geophys. J. Int., 194, 19551965.

32

ACCEPTED MANUSCRIPT Pilz, M., S. Parolai, H. Woith (2017). A 3-D algorithm based on the combined inversion of Rayleigh and Love waves for imaging and monitoring of shallow structures. Geophysical Journal

T

International, 209(1), 152-166.

IP

Podvin, P. and I. Lecomte (1991). Finite difference computation of travel-times in very constrasted

SC R

velocity models: A massively parallel approach and its associated tool. Geophys. J. Int, 271–284. Rawlinson, N., S. Pozgay and S. Fishwick (2010) Seismic tomography: A window into deep Earth. Physics of the Earth and Planetary Interiors 178 101–135.

NU

Rietbrock, A. 2001, P wave attenuation structure in the fault area of the 1995 Kobe earthquake, J. Geophys. Res., 106, 4141–4154, doi:10.1029/2000JB900234.

MA

Rowe, C., R. Aster, B. Brochers, and C. Young (2002), An automatic, adaptive algorithm for refining phase picks in large seismic data sets, Bull. Seism. Soc. Am., 92, 1660–1674.

D

Sambridge, M. (1999). Geophysical inversion with a neighbourhood algorithm—II. Appraising the

TE

ensemble Geophys. J. Int. 138 (3): 727-746 doi:10.1046/j.1365-246x.1999.00900.x Sambridge, M. & Drijkoningen, G., 1992. Genetic algorithms in seismic waveform inversion,

CE P

Geophys. J. Int., 109, 323–342.

Scales, J. A., 1987, Tomographic inversion via the conjugate gradient method, Geophysics, 52/2,

AC

179-185.

Scholz, C.H., Sykes, L.R., and Aggarwal, Y.P., 1973, Earthquake prediction: A physical basis: Science, v. 181, p. 803–809, doi: 10.1126/science.181.4102.803 Serlenga, V., S. de Lorenzo, G. Russo, O. Amoroso, S. Garambois, J. Virieux, and A. Zollo (2016), A three-dimensional QP imaging of the shallowest subsurface of Campi Flegrei offshore caldera, southern Italy, Geophys. Res. Lett., 43,11,209–11,218, doi:10.1002/2016GL071140 Serra, M., G. Festa, P. Roux, M. Gresse, J. Vandemeulebrouck, A. Zollo (2016). A strongly heterogeneous hydrothermal area imaged by surface waves: the case of Solfatara, Campi Flegrei, Italy. Geophys. J. Int., 205, 1813-1822. Teng, T, 1968. Attenuation of body waves and the Q structure of the mantle. Journal of geophysical research, 73, No. 6., 2195-2208, doi: 10.1029/JB073i006p02195.

33

ACCEPTED MANUSCRIPT Tomey and Fougler (1989). Tomographic inversion of local earthquake data from the HengillGrensdalur central volcano complex, Iceland, J. Geophys. Res., 94, 17,497–17,510. Thompson R.C. (1989). Structural stratigraphy and intrusive rocks at The Geysers geothermal field.

IP

T

Geothermal Resources Council, Transaction Vol.13.

Thurber, C. H., Atre, S. R., & Eberhart‐Phillips, D. (1995). Three‐dimensional Vp and Vp/Vs

SC R

structure at Loma Prieta, California, from local earthquake tomography. Geophysical Research

NU

Letters, 22(22), 3079-3082.

Tukey, J. W. (1974). Introduction to today´s data analysis. Proc. Conf. on critical evaluation of

MA

chemical and physical structural information, National Academy of Sciences, Washington, 3-14. Vanorio, T., J. Virieux, P. Capuano, and G. Russo (2005), Three‐dimensional seismic tomography from P wave and S wave microearthquake travel times and rock physics characterization of the

TE

10.1029/2004JB003102.

D

Campi Flegrei Caldera, Journal of Geophysical Research: Solid Earth, 110(B3), doi:

CE P

Wadati K. (1993). On the travel time of earthquake waves, Part II. Geophys. Mag. 7, 101-111. Waldhauser, F., Ellsworth, W.L., 2000. A double difference earthquake location algorithm: Method and Application to the Northern Hayward Fault, California. Bull. Seism. Soc. Am., 90, 6, 1353-

AC

1368, doi: 10.1785/0120000006.

Wathelet, M., Jongmans, D. & Ohrnberger, M., 2004. Surface-wave inversion using a direct search algorithm and its application to ambient vibration measurements, Near Surf. Geophys., 2, 211–221. Whitcomb, J.H., Garmany, J.D., and Anderson,D.L., 1973, Earthquake prediction: Variation of seismic velocities before the San Fernando earthquake: Science, v. 180, p. 632–635, doi: 126/science.180.4086.632. Whitley, D., 1994. A Genetic Algorithm Tutorial, Samizdat Press (http://samizdat.mines.edu/ga tutorial). Yanovskaya, T. B., P. G. Ditmar (1990). Smoothness criteria in surface wave tomography. Geophys. J. Int., 102, 63-72.

34

ACCEPTED MANUSCRIPT Zhang, J., Toksoz, M.N., 1998, Nonlinear refraction traveltime tomography, Geophysics, 63/5, 1726-1737. Zollo, A., L. D’Auria, R. D. Matteis, A. Herreo, J. Virieux, and P. Gasparini (2002). Bayesian

T

estimation of 2-d p-velocity models from active seismic arrival time data: imaging of the shallow

IP

structure of mt vesuvius (southern italy). Geophys. J. Int. 151, 566–582.

SC R

Zollo, A., N. Maerklin, M. Vassallo, D. Dello Iacono, J. Virieux, and P. Gasparini, (2008). Seismic reflections reveal a massive melt layer feeding Campi Flegrei caldera. Geophys. Res. Lett.,

AC

CE P

TE

D

MA

NU

35, L12306, doi:10.1029/2008GL034242.

35

ACCEPTED MANUSCRIPT Highlights

T

IP

SC R NU MA D TE

•

CE P

•

Innovative methodologies for seismic monitoring of volcanic structures in space and time. Elastic and anelastic imaging techniques to characterize medium properties using seismic waves from active and/or passive sources. Campi Flegrei caldera (southern Italy) and The Geysers geothermal area (California) as case studies to evaluate performance and applicability of tomographic methods.

AC

•

36