Here, we apply our method to the red zone19,20 of Campi Flegrei, Italy, to produce a high-resolution amplification map of the region, to show the method's ability to generate high-resolution ShakeMaps after a seismic event, and to prove that our GMM outperforms existing regional GMMs based on VS30 information.
The red zone of Campi Flegrei is defined in the National Emergency Planning for Volcanic Risk in the Phlegraean area and was identified in the Decree of the President of the Council of Ministers of 24 June 2016. The red zone, encompassing an area of ~130 km², is characterised by elevated risk of pyroclastic flows and is inhabited by around 500,000 individuals. In the event of an emergency, the population is to be evacuated as a precautionary measure. The red zone is also the Campi Flegrei area most affected by local earthquakes, and the Italian area with the highest concentration of EQN participants in 2024. From April 2024 to June 2024, the number citizens in the red zone with the EQN app installed was between 7000 and 9000. The red zone is also monitored by the seismic stations of the Italian Seismic Network, of the Italian Strong Motion Network and of the Irpinia Seismic Network (see Supplementary Fig. 1 and Supplementary Dataset 1). Figure 1 depicts the active seismic stations and active EQN smartphones on 8 June 2024. While the density of stations is relatively high with respect to the European density, many highly populated areas of the red zone are only covered by smartphones. By integrating the accelerometric measurements from citizens' smartphones, we improve ShakeMaps right where people live.
We considered four seismic events (detailed in Supplementary Dataset 2) that occurred between 27 April 2024 and 8 June 2024 during the Campi Flegrei seismic sequence started in 2023. These events were characterised by a similar coverage (in terms of network extension) of the red zone by stations and smartphones (see Fig. 2). Similar coverage is important to ensure that site amplification quantified at both the smartphones and the stations has the same reference (average amplification of the covered area). In addition, the events ranged in duration magnitude (Md) from 3.6 to 3.9. The relatively low magnitudes allow the assumption of radial decay without directivity effects. The number of stations in the red zone that measured each event ranged from 27 to 29, while the number of smartphones ranged from 56 to 441. The reason for this high variability is that only smartphones that are being charged perform seismic monitoring, and many more smartphones are charged at night. For stations, we considered PGA measurements, while for smartphones the peak smartphone acceleration (PSmA) as defined in the Methods section. Measurements collected during the four events are shown in Fig. 2.
Our method is motivated by the fact that stations and smartphones measure different physical quantities. Stations directly measure the ground acceleration, while the acceleration measured by smartphones is linked to the ground acceleration through a transfer function that depends on the characteristics of the building in which the smartphone is located, the exact location of the smartphone within the building, and the characteristics of the object over which the smartphone is placed. Evaluating which factors are precisely controlling this bias between smartphone and seismological stations records is challenging. Moreover, when a citizen joins the EQN initiative, all of the above information is not available and does not have to be provided by the citizen to protect his/her privacy.
In the context of a citizen science initiative, the fact that we do not have access to citizens' smartphones, and that smartphones and stations are never co-located, makes it difficult to experimentally assess the bias of smartphone measurements. To overcome this problem, we adopted a statistical modelling approach. For each of the four seismic events, we fitted a spatial statistical model (Eq. 1 in the Methods section) on PSmA measurements. Then, we used the model to predict the PSmA at the station locations, where PGAs are available. The spatial model averages the building and smartphone-specific characteristics, allowing the average PGA-PSmA bias function to emerge.
Using linear regression, we assessed the relationship between the station PGAs and the predicted/filtered PSmA at the station locations. The regression lines depicted in Fig. 3, none of which are parallel to the bisector line, indicate that the average bias is amplitude dependent. On average, peak smartphone accelerations (PSmAs) are about 1.67 times higher than PGAs when PGA is 0.1g, and about 3 times higher when PGA is 0.01g. This bias prevents the direct integration of smartphone measurements into ShakeMaps, justifying the alternative strategy outlined in the following section.
Although PSmAs exhibit a bias with respect to PGAs, we contend that the information from smartphones can be incorporated into ShakeMaps indirectly by modelling the spatial cross-correlation between the two types of measurements. Specifically, for each seismic event, we model the cross-correlation between PGAs and PSmAs residuals after an event-specific isotropic decay has been independently removed (see Supplementary Fig. 2). Station and smartphone residuals when considering the four events together are shown in Fig. 4. Smartphone residuals show higher variability at the small spatial scale because they are affected by both site amplification and smartphone- and building-specific characteristics. Station residuals depend only on the site amplification.
To demonstrate that station and smartphone residuals are spatially cross-correlated, we used the models of Eqs. (1) and (2) to derive a station-only-based log-amplification map and a smartphone-only-based map. Supplementary Fig. 3 depicts both maps, along with the corresponding standard deviation maps and the maps of the areas where the log-amplification is significantly positive or negative considering a 95% confidence interval (Eq. 10). The pixel-by-pixel spatial correlation between the two maps is 0.84. This result is significant as it demonstrates that station and smartphone residuals (when derived and modelled independently) capture similarly the local site amplification. The area where the station-only-based log-amplification map is significantly positive or negative corresponds to around 38.9 km² (29.8% of the red zone area) while the area for the smartphone-only-based map is around 64.5 km² (49.0%).
We then used the model in Eq. (11) to obtain the high-resolution log-amplification map based on the data fusion of station and smartphone measurements. The back-transformed (exponential of the log-amplification) map is depicted in Fig. 5, while the uncertainty map of the log-amplification is shown in Supplementary Fig. 4. We stress that the uncertainty map is not a secondary product as it enables the identification of areas where the amplification map can be relied upon, and of areas where the amplification significantly differs from the average amplification of the region. The level of uncertainty depends on the number of stations and smartphone measurements, their spatial distribution across the red zone, the number of seismic events used to derive the amplification map, and the errors and biases in the smartphone measurements.
Using the data fusion model, the area where the log-amplification is significantly positive or negative is 79.0 km² (61.0% of the red zone area). Combining station and smartphone data increased the area by more than double compared to the map based on station data alone. Additionally, we assessed the power of the estimated log-amplification map in reducing the station residual variability using Eqs. (12) and (13). The standard deviation of the average station residual once the isotropic decay is removed is 0.6389 (0.2775 on the log scale). Correcting the residual by the log-amplification map at the station decreased the standard deviation to 0.1827 (0.0793 on the log scale), corresponding to a 71.4% reduction.
The log-amplification map and its uncertainty are also provided in Supplementary Dataset 3. The estimated model parameters for the four events are given in Supplementary Table 1.
Considering the high-resolution amplification map of Fig. 5, and the population and building exact locations from the European buildings dataset, we estimated the exposure distribution to site amplification for 498,817 people, for 8151 EQN users and for 334,837 buildings in the red zone. We found that all distributions (see Supplementary Fig. 5) are skewed toward amplifications lower than 1. This implies that people and buildings tend to be located in areas of the red zone where the PGA is deamplified rather than amplified (relative to the PGA estimated by an isotropic decay centred on the epicentre). About 19.4% of the population is exposed to an amplification factor of more than 1.5 while 5.8% to a factor of more than 2. About 24% of the buildings are exposed to an amplification factor of more than 1.5 and about 5.4% to a factor of more than 2. The EQN users and the population have a very similar exposure distribution, which means that the EQN users were a representative sample of the population (in terms of site amplification exposure).
The high-resolution amplification map issued from the fusion of station and smartphone data was used to produce a high-resolution PGA ShakeMap for the Md4.2 earthquake that occurred on 27 September 2023 with epicentre in the red zone (event detail in Supplementary Dataset 2). This was done using the model in Eq. (14). The event is the largest (in terms of magnitude) of the ongoing 2023-2025 Campi Flegrei sequence that was detected by at least 100 smartphones of the EQN initiative. Figure 6 depicts the PGA ShakeMap and its uncertainty derived using Eqs. (16) and (17), respectively. To better appreciate its high spatial resolution, the ShakeMap is restricted to an area of the red zone where the estimated PGA is significantly higher than zero. For a comparison, Supplementary Fig. 6 shows the ShakeMap provided by the National Institute of Geophysics and Volcanology (INGV). The PGA decay of the INGV map is nearly isotropic. On the contrary, Eq. (14) produced a ShakeMap that exhibits small-scale variability thanks to the amplification map and thanks to the event-specific smartphone measurements.
As a further relevant application, we used the high-resolution amplification map to improve the recently published regional GMM with site classes based on VS30. We considered a total of 661 station PGA measurements from 34 seismic events that occurred between 5 February 2023 and 18 June 2024, with an epicentre in the red zone and a magnitude between Md2.9 and Md4.4 (details in Supplementary Dataset 4). We then compared three GMMs: the regional GMM, the same GMM recalibrated on the 34 events, and our GMM of Eq. (19) based on the amplification map in Fig. 5 and calibrated on the same events. The standard deviation of log-PGA residuals is around 0.410 for the original GMM, around 0.397 for the recalibrated GMM and around 0.357 for our GMM. This is a 13% and 10% reduction, respectively, which is significant given the high impact of random variability in GMMs on probabilistic seismic hazard assessments. Additionally, the reduction in the standard deviation of the residuals is a measure of the GMM enhancement and an implicit validation of the amplification map. The estimated GMM coefficients are given in Table 1.
We then used our calibrated GMM to simulate the PGA ShakeMap assuming a Md5.0 event with epicentre in 40.82° latitude and 14.15° longitude. This is the same epicentre of the Md4.6 event of 13 March 2025. The ShakeMap, which is depicted in Fig. 7, retains the high spatial resolution and the spatial patterns of the amplification map. The PGA ShakeMaps obtained from the other GMMs are depicted in Supplementary Figs. 7 and 8, and they exhibit a radial decay with interaction of the VS30 site classes.