Anatomy of the Campi Flegrei caldera using Enhanced Seismic Tomography Models

Campi Flegrei caldera (Southern Italy) is a densely inhabited area and suffered several unrest episodes in the last centuries. The dynamic of the caldera is highly debated because of conflicting interpretations. Here we present a detailed reconstruction of the Campi Flegrei structure obtained using the microseismicity recorded during the 1984 unrest. Enhanced Seismic Tomography models obtained with these data allow us describing seismic velocities, attenuation, and scattering patterns. Results show: (1) a plumbing system with a diameter of 1 km located between 2.3 km and 4 km depth (2) a 0.5 km thick caprock located at 2 km depth interpreted as the main structure regulating the fluid interchange between deep and shallow sectors of the caldera, (3) the shape and volume of a shallow reservoir beneath the city of Pozzuoli; this reservoir played a key role during the 1982–1984 unrest, (4) several small reservoirs beneath the main craters of the caldera. All these features fit into the debated question on magmatic or hydrothermal mechanism driving the caldera deformation resulting of crucial importance to allow a better assessment of the hazard.


Anatomy of the Campi Flegrei caldera using Enhanced Seismic
Tomography Models

1) Maps with the final location of the seismicity.
Figure S1.1. Maps of the whole seismicity located with the final 3D seismic models and slices at different depth ranges.  Figure S1.2c shows the horizontal and vertical projections of the P-wave ray path coverage using the new preliminary locations. The average RMS of the relocated events on the final 3D WAM models is of 0.048s, resulting in an improvement of the hypocenter pattern that is twice with respect to that obtained with the previous models that reached a minimum value of 0.08s [2, 3, 4].

2) Set up of the 1D models Vp , Vs and
Figure S1.2. a) Layered Vp (black) and Vs (red) models used for the preliminary location and interpolated models (blue) used as starting velocities for the double difference tomography inversions. b) Layered Vp/Vs distribution in depth (black) and corresponding interpolated model (blue). c) Horizontal and vertical projection of the P-wave ray path coverage of the studied area.
3) Set up of the WAM models from the inversion of 13 perturbed grids. Figure S1. 3 shows the horizontal and vertical density sampling of the investigated volume obtained merging the 13 grids used for the 13 Double Difference tomographic inversions that have been generated to build the WAM. It is worth noting that the shift of the grid in depth produces small perturbations of the initial Vp and Vs models increasing the variability of the space of the models sampled during the inversions. The figure shows the difference in the density of points sampled with the WAM rather than using a grid alone.

4) Assessment of the Vp and Vs models.
To assess the reliability of our results and to show the resolution power of the data and method, we constructed a checkerboard model characterized by alternating positive and negative velocity changes of ±10 percent with respect to the initial 1-D velocity model. Each patch has size of 0.5 km along the three dimensions. The model was then used to calculate synthetic travel times using the same configuration of earthquakes and stations as in the real inversion. Travel time errors were simulated by adding Gaussian distributed noise with a standard deviation of 0.04s and 0.06s to the P and S arrival times, respectively, corresponding to a 2σ of the largest expected picking error for this dataset

5) Spike tests
To further assess the reliability of the structures observed in the models obtained using the real data, we performed two discrete spike tests or sparse checkerboard tests [8] that involve a sparse distribution of spikes (figure S1.5a and S1.5b). These tests were performed to evaluate 1) the magnitude of smearing for near features [8], and 2) the level of data noise needed to degrade the tomographic images with the consequence of limiting the interpretation of the models. Figure S1.5a shows the spike test performed with levels of noise of 0.04s and 0.06s for the P and S arrival times, respectively. The method allows to finely recover the spikes and, as evident, the smearing effects are very reduced. These levels of noise are 2σ of the largest expected picking error for this dataset [3]. Figure S1.5b shows that only with travel time errors about 8 times greater than the ones expected for the CF (4 if we consider a 2σ of the errors), the models start to be degraded to levels where the reconstruction of the spikes is poor, especially at depths grater than 2.5 km. However so big errors on the pickings (0.08s for P and 0.14s) would imply a miscomputation in the travel times of almost 10% of their values and an initial error of the event locations of more than 2 km. These conditions are unrealistic for the CF database. Therefore, this test demonstrates that with the used data it is difficult to introduce biases into the velocity models large enough to lead to a misinterpretation of the recovered structure. Figure S1.5a. Horizontal sections of the P-and S-wave discrete spike tests using a realistic level of noise of the data. Top left models are the pattern of the true model. Note that the layer at 3.5 km depth has a pattern inverted with respect to the other ones. The true model was built with this feature to further increase the complexity of the model. Figure S1.5b. Horizontal sections of the P-and S-wave discrete spike test using level of noise of the data four times than the expected 2 errors. Top left models are the pattern of the true model. Note that the layer at 3.5 km depth has a pattern inverted with respect to the other ones. The true model was built with this feature to further increase the complexity of the model 6) Restoration-Resolution test Finally, to assess the reliability of the deep structure observed, we built a synthetic model characterized by a P-and S-wave velocity anomaly 25% higher with respect to the initial 1-D velocity model. The anomalous body has the shape of a rectangular prism with horizontal section of 1x1 km 2 at 5 km depth flatting at 3 km depth into a shallow body 0.5 km thick elongating about 2 km S-W and 0.7 km large (figures S1.6a, b, and c). This model can be considered as an approximation of the velocity anomalous structure that has been observed in the inversions done with the real data, allowing us to perform a Restoration-Resolution test [9]. With the same configuration of earthquakes and stations as in the real data inversion, we calculated synthetic travel times, adding picking errors of…., and running an inversion procedure as for the real data. Results show that the data and method used are able to fairly recover such a complex structure down to 4.5 km depth for both Vp (figures S1.6d, e, and f) and Vs models (figures S1.6g, h, and i). Following all the tests and considerations done, we can affirm that the velocity models presented here are able to fairly reconstruct features with linear dimension of 0.5 km in the CF for depths ranging from 0 km to about 4.5 km. The possibility that these models contain artifacts that could be misinterpreted is considered low in the best resolved regions. . In this work we adopted a threshold of 100 for delimiting the best-resolved volume (figures S1.7a and S1.7b), which is twice that generally used in other attenuation studies of similar size and amount of data (e.g.

7) Assessment of the model
[10]). Furthermore we interpret only the part of the models falling into the regions where the DWS is greater than 500 to ensure the best reliability of our interpretation.  The reliability of the scattering model has been assessed using the data coverage for each grid cell, i.e. the number of possible scatters associated with each volume during the analysis. Figure S1.7c shows the number of scatters obtained for the envelopes filtered at the dominant frequencies around 12 Hz calculated at different depths delimiting the region of influence of the asperities. Figure S1.7c. Number of scatters calculated at different depths.