Spatial analysis for an evaluation of monitoring networks: examples from the Italian seismic and accelerometric networks

Abstract

In this work, we propose a statistical approach to evaluate the coverage of a network based on the spatial distribution of its nodes and the target information, including all those data related to the final objectives of the network itself. This statistical approach encompasses descriptive spatial statistics in combination with point pattern techniques. As case studies, we evaluate the spatial arrangements of the stations within the Italian National Seismic Network and the Italian Strong Motion Network. Seismic networks are essential tools for observing earthquakes and assessing seismic hazards, while strong motion (accelerometric) networks allow us to describe seismic shaking and to measure the expected effects on buildings and infrastructures. The capability of both networks is a function of an adequate number of optimally distributed stations. We compare the seismic network with the spatial distributions of historical and instrument seismicity and with the distribution of well-known seismogenic sources, and we compare the strong motion station distribution with seismic hazard maps and the population distribution. This simple and reliable methodological approach is able to provide quantitative information on the coverage of any type of network and is able to identify critical areas that require optimization and therefore address areas of future development.

1 Introduction

Country-scale seismic networks usually extend over an area of 10⁵–10⁶ km² (D’Alessandro et al. 2019) and, from planning to completion, can require several years to decades. In addition, the final spatial distribution of the nodes can be remarkably different from the planned configuration for a multitude of reasons.

As some examples, the Italian National Seismic Network, (D’Alessandro et al. 2011a, 2019) the Red Sísmica Nacional in Spain (D’Alessandro et al. 2013), the Romanian network, and the Global Seismographic Network (GSN) took on average 20–30 years (or more) to reach their present-day configuration, and they are still developing (Nacional and Rodríguez 1995; Gee and Leith 2011; Popa et al. 2015; Michelini et al. 2016; D’Alessandro et al. 2011b, 2012; D’Alessandro and Ruppert 2012; D’Alessandro and Stickney 2012). The development process for most seismic (or strong motion) networks is often analogous to that for these large networks. In other cases, a network could be the result of the joining of different subnetworks. This is the case for the Hellenic Unified Seismological Network, which originated from four distinct networks, each of which developed independently in time and with partial spatial overlap (D’Alessandro et al. 2011b). In most cases, the long span of time needed to establish a network is also ascribable to successive technological improvements. Therefore, it is always beneficial to verify whether the distribution of nodes within a given network satisfies its final needs and objectives.

In this work, we propose a statistical approach based on the spatial distribution of nodes within a large-scale network together with ancillary information related to the aims of the network itself (e.g., seismicity, seismogenic sources, and hazards). We aim to provide a simple and reliable methodological approach that can assess the degree of coverage of a network and address future optimizations.

2 Methodology

In this paper, data are analyzed by means of descriptive spatial statistics and point process methods. The methodologies are briefly described in the following paragraphs in terms of their formulas, meaning, and usage with particular focus on the proposed applications. The analysis is performed with R statistical software (R Development Core Team 2005), and the functions employed herein are in the packages raster (Hijmans 2018), spatstat (Baddeley and Turner 2005), and spatstat.local (Baddeley 2018).

2.1 Point process methods

A point process model $\mathcal {X}$ is a random finite subset of $W \subset \mathbb {R}^{d}$, where $|W|< \infty$. A spatial point pattern $\boldsymbol {\upsilon }=\{\textbf {u}_{1},\dots ,\textbf {u}_{n}\}$, which is a realization of $\mathcal {X}$, is an unordered set of points in the region W, where N(υ) = n is the number of points and d = 2.

For a spatial point process, the first-order intensity (λ(u)) is the expected number of points in a small region around a location (u) divided by its area (du) in the limit as du goes to 0:

$$ \lambda(\textbf{u}) = \lim\limits_{|d\textbf{u}| \rightarrow 0} \frac{\mathbb{E}(N(d\textbf{u}))}{|d\textbf{u}|} $$

(1)

The intensity (λ(u)) is assumed to be inhomogeneous, so it is not constant over the region; for spatial point patterns, the intensity is estimated non-parametrically to understand the spatial trend. The usual kernel estimator of the intensity function is computed as proposed by Baddeley et al. (2015):

$$ \hat{\lambda}(\textbf{u})=\frac{1}{e(\textbf{u})} \sum\limits_{i=1}^{n}k(\textbf{u}-\textbf{u}_{i},\boldsymbol{h}) $$

(2)

where e(u) is the edge correction and k(⋅) is a bivariate Gaussian density distribution function with a standard deviation (smoothing bandwidth) equal to h = (h_x,h_y). The bandwidth controls the degree of smoothing, and a larger set of values gives more smoothing. The estimation of the bandwidth vector h = (h_x,h_y) is computed with Scott’s rule (Scott 1992):

$$ h_{x}=n^{-1/6} \sqrt{var(X)} \qquad h_{y}=n^{-1/6} \sqrt{var(Y)} $$

(3)

In the proposed approach, the datasets concerning the nodes of the seismic and accelerometric networks and the instrumental and historical seismicity are treated as point patterns, and their spatial intensities are properly computed using the nonparametric estimator in Eq. 2.

Generally, it can be valuable to describe the first-order characteristics of a point pattern by studying its relationship with external variables (also called covariates). In particular, for a seismic network, the intensity is computed depending on an external covariate, namely, the distance to the nearest geological fault, D(u), to understand the spatial displacement of the nodes relative to a seismic source. An inhomogeneous Poisson model with a parametric log-linear form is assumed following (Baddeley et al. 2015):

$$ \lambda(\textbf{u})=exp(\beta_{0}+\beta_{1} D(\textbf{u})) $$

(4)

where D(u) is a spatial covariate and β₀ and β₁ are the parameters to be estimated. The parameters in Eq. 4 can be interpreted as follows: $\alpha = e^{\beta _{0}}$ gives the intensity at locations close to a geological fault (i.e., D(u) ≈ 0), while the intensity changes by a factor of $e^{\beta _{1}}$ for every 1 unit of distance away from the seismic source.

Moreover, we locally infer the spatially varying estimates of the parameters of the inhomogeneous Poisson process model in Eq. 4 following Baddeley (2017):

$$ \lambda(\textbf{u})=exp(\beta_{0}(\textbf{u})+\beta_{1}(\textbf{u}) D(\textbf{u})) $$

(5)

This approach, also called geographically weighted regression, has the potential to detect and model the gradual spatial variations in the parameters that govern the intensity of the point pattern as a function of the spatial covariate.

2.2 Descriptive spatial statistics

In the second step of the analysis, the intensities of the two networks are related to the ancillary information to assess their coherence. The estimated intensities of the point patterns, the hazard maps, and the population distribution are considered raster data. Pairs of maps are then compared to determine if the spatial information within the raster data is locally correlated. The local correlation coefficient for a pair of raster maps is computed considering a grid with a resolution of 5 km and a neighborhood of 5×5 cells (25×25 km²) around each cell. This resolution is calibrated according to the unique scale of each subsequent case study, but the resolution can be advantageously set to adapt to other applications.

In addition, we compute some descriptive statistics to obtain further information about the distributions of both networks analyzed herein with respect to the ancillary datasets.

3 Data

As case studies, we consider two different monitoring networks in Italy: one devoted to earthquake detection and characterization (i.e., the Italian National Seismic Network) and another devoted to measuring the effects of earthquakes (i.e., the Italian Strong Motion Network). For both networks, we also consider ancillary information that is strictly related to their respective objectives and is described in the following paragraphs.

3.1 The networks

The Italian National Seismic Network was developed immediately after the Irpinia seismic crisis in 1980. Initially, the National Seismic Network (NSN) consisted of only a few stations spread throughout Italy. In the 1990s, the ING (National Institute of Geophysics), later known as the INGV (National Institute of Geophysics and Volcanology), progressively upgraded the network; after approximately 30 years, these upgrades led to the present-day configuration, which consists of approximately 500 seismic stations (Michelini et al. 2016). All the collected data are managed in real time by the INGV seismology center in Rome. Over the years, the network grew in terms of its spatial coverage and sensor quality. However, the network did not develop with uniform coverage or with a focus on seismically active areas; rather, the evolution of the network was controlled mainly by logistic needs and the availability of financial resources. The result was a patchy network; and thus, the quality of earthquake locations and the magnitude of completeness are very heterogeneous (Marchetti et al. 2004; Amato and Mele 2008; Schorlemmer et al. 2010; D’Alessandro et al. 2011a; Chiarabba et al. 2015).

The effective development of the monitoring infrastructure of a network should strictly focus on the areas with the greatest potential seismic release (as determined from instrumental and historic seismicity) and on the areas where the seismogenic sources are well defined. A reasonable correlation between the seismicity rate and network coverage is currently lacking. This can clearly affect the quality of seismic monitoring and related seismological research because, as is well known, the quality of the estimated focal parameters depends not only on the quality of the recorded waveforms but also on the density and geometry of monitoring stations.

With similar criteria as the Italian National Seismic Network but a different purpose, the Italian Strong Motion Network (Rete Accelerometrica Nazionale, RAN) was developed by the Italian Department of Civil Protection (Dolce 2009). The network consists of 528 stations, each of which is equipped with an accelerometer; the main objective is to assess seismic responses on the basis of recorded ground motions. The strong motion network provides information about the extent of the area over which an earthquake is felt and about the post-event effects. The development of the RAN followed criteria that changed over time. The early sites were selected in areas with medium–high seismicity, and the stations were occasionally installed on outcropping, non-fractured rocky substrates to avoid amplification. Later stations were deployed according to a more uniform distribution criterion on a regular grid of approximately 20–30 km. However, after the 2009 L’Aquila seismic sequence, when the lack of information did not allow a proper assessment of the post-event situation, some relevant weaknesses of the network were realized. As a general rule, a strong motion network should be strengthened in areas with relatively strong expected shaking and in areas that are both vulnerable and exposed (i.e., urban areas in earthquake-prone zones).

Within the study window (W), 465 and 569 stations belong to the NSN and RAN, respectively. The data are available at http://cnt.rm.ingv.it/instruments (NSN) and http://www.protezionecivile.gov.it/jcms/it/ran.wp (RAN).

3.2 Ancillary information

The study area (W) comprises the Italian Peninsula and the island of Sicily. Because only a few stations in both networks (seismic and accelerometric) are located in offshore areas and on Sardinia, those areas are not included. Moreover, the seismic release and seismic hazard in these areas are negligible. The ancillary data are subdivided into two groups. The first group is related to the aims of the seismic network and includes information about the seismicity and seismogenic sources. The second group is related to the aims of the accelerometric network and includes seismic hazard maps and the distribution of the population in Italy. All the datasets are briefly described in this section.

The Italian seismic catalogue contains events since 1985. We consider earthquakes until the end of 2018 with magnitudes greater than 3. This threshold is necessary because the earlier configurations of the seismic network (i.e., with sparse distributions of stations) were not able to regularly detect small earthquakes. However, we can be fairly confident that earthquakes of M > 3 were all detected (Marchetti et al. 2004; Amato and Mele 2008). A subset of this catalogue containing 7141 events falling within the study area is thus analyzed. Furthermore, historical seismic events are selected from the Parametric Catalogue of Italian Earthquakes (Rovida et al. 2019, 2020), which contains 3268 events during the period 1000–2014 within W.

Additionally, information about seismic sources should be considered to evaluate the locations of the seismic stations. Thus, the composite seismogenic source (CSS) dataset extracted from the Italian Database of Seismogenic Sources (DISS) (Group et al. 2010) is considered. A composite source represents a complex fault system that includes an unspecified number of aligned individual seismogenic sources that cannot be separated spatially. In particular, we consider the projection of a CSS on the horizontal plane; for the sake of simplicity, in the text, the CSS is referred to as a “fault.” We assume that the probability of earthquake occurrence is uniformly distributed on each fault plane and consequently on its projection at the surface.

The seismic hazard map of Italy is a probabilistic model assessing the likelihood that a given ground acceleration threshold is overcome in a given time span. The hazard map is calculated taking into account the frequency-magnitude record of Italian earthquakes together with the attenuation laws of ground motion. The most common way to display a seismic hazard model is as a map relative to the 10% exceedance probability for the maximum expected ground acceleration in 50 years, which corresponds to a return period of 475 years (PCM 2006). For our analysis, we also consider the two available extreme hazard maps, which are relative to exceedance probabilities of 81% and 2% in 50 years, corresponding to return periods of 30 and 2475 years, respectively (PCM 2006).

Finally, we consider information about the distribution of the Italian population. The data come from the most recent official census preformed by the Italian Institute for Statistics (ISTAT) in 2011 (www.istat.it).

The seismic stations, accelerometric stations, and both earthquake catalogues (instrumental and historical seismicity) are treated as three spatial point patterns in W. The other data are treated as raster maps with a cell size resolution of 6 km.

4 Network spatial analysis

For simplicity, the study area is divided into labeled subregions that facilitate further analysis and interpretation of the results (Fig. 1). The boundaries between the zones should not be considered as lines separating areas with different characteristics but rather indefinite boundaries drawn simply to facilitate discussion of the results.

Fig. 1

Subdivision of the study area into subregions labeled with abbreviations. NW northwest, NE northeast, ENE eastern northeast, LNE lower northeast, UCW upper central west, UCE upper central east, LCW lower central west, LCE lower central east, SW southwest, SE southeast, C Calabria, and S Sicily. The plot refers to the Mercator projection UTM33 in km

4.1 Spatial distributions of the networks

Given the spatial point patterns of the two networks (seismic and accelerometric), we estimate the inhomogeneous intensity with a kernel estimator (Eq. 2). In our analysis, W is divided into a 6-km square grid. Integrating the intensity over the entire study area, the expected number of points is obtained; therefore, the values reported in the plots (Fig. 2) indicate the number of stations over 36 km².

Fig. 2

Kernel intensity estimation (with the smoothing bandwidth selected by Scott’s rule) for the a seismic monitoring network and b accelerometric monitoring network. The intensity values correspond to the number of stations in a 6×6 km cell. Black points denote stations. The plots refer to the Mercator projection UTM33 in km

The intensity maps are shown in Fig. 2. Among the seismic stations (Fig. 2a), the map indicates three areas with the highest intensity, namely, the SW, S, and UCE regions (see Fig. 1). The high concentrations of stations in the SW and S regions are ascribable mainly to the resolution necessary to monitor the peculiar characteristics of the seismicity in volcanic areas. Conversely, the high intensity in the UCE region is due to the local oversizing of the network thanks to several project-funded upgrades. The minimum values of the intensity are detected in the NW and SE regions. The intensities of the accelerometric network are more uniform than those of the seismic network at the country scale (Fig. 2b). Locally, there are areas with a higher concentration of stations along the Apennines Mountains, especially in the central and southern parts. In contrast, the minimum intensities are found in the NW, NE, and SE regions (Fig. 2b).

4.2 Seismic network versus seismogenic sources

Here, the spatial relationship between the faults and the nodes of the seismic network is explored. The distance from each node (i.e., seismic station) to the nearest fault is computed, and vice versa, and the corresponding cumulative density functions are shown in Figs. 3a and 3b, respectively. Among all the seismic stations, 33% are located on the surface projection of a fault, and 80% are less than 20 km from the nearest fault. Siting stations as close to the hypocenter as posiible is crucial for obtaining reliable depth estimations. In total, 50% of the faults are directly covered by at least one seismic station, and approximately 95% of the faults are monitored from a distance shorter than 20 km.

Fig. 3

a Cumulative density function of the distance to the nearest fault of each station. b Cumulative density function of the distance to the nearest station of each fault

As explained before, the seismic station intensity is computed as a function of an external covariate: the distance to the nearest fault. The main goal is to assess whether the stations are spatially arranged to occur more frequently near faults. The geological information in Fig. 4a is transformed into a spatial variable defined at all locations u ∈ W, namely, D(u), which is the distance to the nearest fault (Fig. 4b). The estimates of the parametric model in Eq. 4 are $\hat {\beta _{0}}=-19.90$ and $\hat {\beta _{1}}=-0.000018$. Generally, by increasing the distance to the nearest fault, the station intensity decreases. However, according to the spatially varying slope coefficient related to the model in Eq. 5, this reduction is higher in central and northern Italy (mainly the NW, LNE, UCW, and LCW regions) (Fig. 5). Moreover, in region C, the sign of the local slope coefficient is positive, indicating that the stations are not located on faults.

Fig. 4

a Polygons indicating the projection of the complex seismogenic sources in the Italian and Mediterranean areas onto the horizontal plane. b The distance in km to the nearest fault for each grid point in the study area

Fig. 5

Estimated components of the geographically weighted regression model in Eq. 5. a Fitted intensity of seismic stations according to the local likelihood fit of the log-linear model. b Spatially varying estimates of the intercept. c Spatially varying estimates of the slope coefficient related to the spatial covariate D(u), that is, the distance to the nearest fault

4.3 Seismic network versus seismicity

Here, the relationship between the spatial distribution of seismic stations and instrumental seismicity is evaluated. The areas with the highest number of events are in central Italy. Most of the seismicity relates to the three main sequences that occurred in 2009 (L’Aquila, LCW region), 2012 (Emilia, LNE region), and 2016–2018 (Amatrice-Norcia-Visso, UCE region) (Fig. 6a). High intensities also characterize the area around Mount Etna (S region), which is not related to a specific seismic sequence but rather to continuous background seismicity (Fig. 6a). The historical seismicity data indicate that most of the seismic activity in Italy has occurred along the Apennines and in Sicily (Fig. 7a).

Fig. 6

a Kernel intensity estimation for the instrumental seismicity between 1985 and 2018 for earthquakes with magnitudes greater than 3 and with the smoothing bandwidth selected by Scott’s rule. b Local correlation coefficients for the raster intensity objects: the seismic stations (Fig. 2a) and the instrumental seismicity (Fig. 6a). Each focal square area has dimensions of 5×5 grids, namely, 30×30 km²

Fig. 7

a Kernel intensity estimation of the historical seismicity with smoothing bandwidth selected with Scott’s rule. (b) Local correlation coefficients for the raster intensity objects: the seismic stations (Fig. 2a) and the historical seismicity (Fig. 7a). Each focal square area has dimensions of 5×5 grids, namely, 30×30 km²

Furthermore, the relationship between the spatial distribution of the stations with respect to the seismicity (both historical and instrumental) is explored at the local scale. The global and local correlation coefficients are computed to compare the pairs of the estimated intensities between Figs. 2a and 6a and between Figs. 2a and 7a. As expected, the overall correlation coefficients between the station intensity and the instrumental seismicity and historical seismicity are positive, where the Pearson correlation coefficients are 0.468 and 0.625, respectively.

Moreover, considering the estimated intensities as raster data, we check for spatial variations in the correlation. Around each cell in both rasters, we define a focal square area of 5×5 cells and record the correlation of the central cell with each of the surrounding cells in the focal area. In both cases, there is a generally positive correlation, while zones with a negative correlation occur only locally (Figs. 6b and 7b). Moreover, the areas identified as having a negative correlation are shared between both maps, where negative correlations characterize parts of the NW, LNE, UCW, LCW, C, and S regions. These negative correlations indicate discordance between the distributions of the stations and the seismicity (both instrumental and historical). It is not possible to discriminate whether these negative correlations are due to the lack of stations in a seismic area or poor seismicity in a well-covered area. However, in both cases, the seismic network in these areas requires an accurate re-assessment.

Furthermore, the cumulative number of events with increasing distance is calculated around each station (Fig. 8a). By selecting the stations with the highest number of events at a distance of 5 km (Fig. 8c) and the stations with the lowest number of events at a distance of 100 km (Fig. 8d), it is possible to distinguish which stations are located the best and worst with respect to the instrumental seismicity. A map of all the stations ranked according to the cumulative number of events at a distance of 50 km is also shown in Fig. 8e.

Fig. 8

a Cumulative number of earthquakes (instrumental seismicity from 1985 to 2018) for each seismic station with increasing distance and b the corresponding intensity. c The top 5% of all stations with the highest cumulative number of earthquakes at 5 km. d The top 5% of the stations with the lowest cumulative number of earthquakes at 100 km. e Stations classified according to the cumulative number of earthquakes at a distance of 50 km. The plots refer to the Mercator projection UTM33 in km

The seismic moment M₀ for historical earthquakes is calculated by converting the provided values of the moment magnitude M_w by means of the relation proposed by Hanks and Kanamori (1979):

$$ {M}_{0}=10^{(M_{w}+10.7)*1.5} $$

(6)

Then, similar to the instrumental seismicity, the cumulative seismic moment M₀ with increasing distance is computed around each station. The obtained curves (Fig. 9a) provide information about the coherence between the location of each station and the historical seismic release. Among the stations, we select the stations having the highest values of M₀ at 5 km and those having the lowest values at 100 km (Figs. 9c and 9d, respectively), which correspond to the “best” and “worst” placed seismic stations with respect to the historical seismic release. A map of all the stations ranked according to the cumulative M₀ at 50 km is also shown (Fig. 9e), providing a comprehensive explanatory framework for the effectiveness of the network.

Fig. 9

a Cumulative seismic moment release (M₀) of historical seismicity for each seismic station with increasing distance and b the corresponding intensity. c The top 5% of the stations with the highest value of the cumulative M₀ at 5 km. d The top 5% of the stations with the lowest value of the cumulative M₀ at 100 km. e Stations classified according to the cumulative M₀ at a distance of 50 km. The plots refer to the Mercator projection UTM33 in km

4.4 Strong-motion network versus seismic hazard maps and population

In this section, the relationship among the spatial distribution of strong motion stations, the hazard maps, and the population is evaluated. Seismic hazard maps provide information about the likelihood of ground shaking in a given area. Hazard maps are fundamental tools for estimating the likelihood of damage and losses in such areas and are the basis for quantitatively assessing seismic risk. The local correlation between the expected peak ground acceleration in 50 years and the accelerometric station intensity (Fig. 10a) is calculated for each 5×5 km² cell considering the surrounding 25 cells (30×30 km²). The results indicate a general positive correlation, whereas the correlation is negative only in restricted areas, namely, parts of the UCW, UCE, SE, C, and S regions (Fig. 10b). The local correlations with the hazard maps at different probability thresholds and return periods are entirely analogous.

Fig. 10

a The ground motion (peak ground acceleration) expected to be reached or exceeded with a 10% probability in 50 years. b Local correlation coefficients for the raster intensity objects: the accelerometric stations (Fig. 2b) and the seismic hazard map (a). The focal square area has dimensions of 5×5 grids, namely, 30×30 km²

One of the main goals of a strong motion network is to assess the earthquake intensity (i.e., the effects at the surface); and thus, the network should also be planned according to the distribution of the population. To evaluate the coherence of this relationship, we calculate the cumulative number of people living around each station with increasing distance (Fig. 11a). By selecting the stations with the highest number of people at a distance of 5 km (Fig. 11c) and the stations with the lowest number of people at a distance of 100 km (Fig. 11d), it is possible to distinguish which stations are placed the best and worst with respect to the population distribution. A map of all the stations ranked according to the cumulative number of people at 50 km is shown in Fig. 11e. Focusing on the main (i.e., most populated) cities (more than 100,000 inhabitants), which alone represent 19.2% of the whole population of Italy (Table 1), we calculate the cumulative number of strong motion stations with increasing distance (Fig. 12).

Fig. 11

a Cumulative number of people around each accelerometric station with increasing distance and b the corresponding intensity. c The top 5% of the stations with the highest value of the cumulative population at 5 km. d The top 5% of the stations with the lowest value of the cumulative population at 100 km. e Stations classified according to the cumulative population at 50 km

Table 1 List of the main Italian cities with more than 100,000 people and the counts of accelerometric stations at increasing distances (10, 20, and 30 km). PGA refers to the expected peak ground acceleration (in g) from the seismic hazard map and the corresponding classification (see Table 2)

Fig. 12

Cumulative number of accelerometric stations with increasing distance from the center of major cities (population greater than 100,000) listed in Table 1. The black solid line is the mean curve, the black dashed lines are drawn considering the standard deviation, and the red dashed lines correspond to regular grid distributions of 20, 30, and 40 km

5 Discussion

Assessments of the coverage of seismic and strong motion networks (and consequently their performance) are not only desirable but also of crucial importance to find areas where the networks need to be concentrated or, conversely, where the networks are redundant. In this work, we approach this topic by means of descriptive spatial statistics and point process methods. It should be noted that the analyses herein are limited to the mainland of Italy and the island of Sicily, as taking into account the stations located on the various islets along the coast would have led to misinterpretations of the results of the applied statistical methods.

The main purpose of a seismic network is to accurately characterize seismicity in terms of hypocenter locations and magnitude estimations. For this reason, we cross-checked the distributions of the networks with the distributions of faults and seismicity (both instrumental and historical). More than 95% of the stations are located with 50 km of a fault; this can be considered a valid distance for the adequate detection of seismicity exceeding M ≈ 2.0. In general, all the mapped faults are closely (< 25 km) monitored by at least one station; however, this distance is shorter than the longitudinal length of several faults, and therefore, the information about the coverage could be partial. Further insights are provided by the spatial correlations between the intensity of the seismic network (Fig. 2a) and the fault distribution (Fig. 4b). We found, as a general rule, a negative correlation: the greater the distance from the fault, the fewer seismic stations there are. At the local scale, the highest coherence between the station intensity and fault distance was detected in the ENE, UCE, and LCE regions; conversely, the C region showed disagreement (Fig. 5c). The positions of the stations, when compared with the distribution of instrumental seismicity, indicate truly scattered behavior (Fig. 8a and b). The better-placed stations (more events at a distance of 5 km) fall within the UCE region and in the vicinity of Mount Etna (S region), while the less coherent stations (fewer events at a distance of 100 km) fall in the NW, SE, and western S regions (Fig. 8c and d, respectively). Of course, instrumental seismicity alone does not provide a thorough framework; therefore, considerations of historical earthquakes are essential. Again, some nodes are better placed than others. With increasing values of cumulative M₀ (Fig. 9a and b), the better-placed stations (highest M₀ at 5 km) fall within the UCE and SW regions, and a few others are scattered throughout the whole territory (Fig. 9c). The stations exhibiting lower coherence with the historical seismic release (lowest M₀ at 100 km) fall within the NW and SE regions (Fig. 9a). The computed local correlations between the station intensity and distribution of instrumental seismicity indicate agreement between the two compared rasters. In fact, the map is not patchy as we could expect; rather, variations occur at a large scale. In general, the correlations are positive even though some small areas with negative correlations occur, namely, the UCW, NE, SE, and western S regions (Fig. 6b). Unfortunately, the main limitation of the proposed methodology is the impossibility of distinguishing between the two possible combinations, leading to the same conclusion: positive correlations in areas with either many stations and high seismicity or few stations and low seismicity and, conversely, negative correlations in areas with either few stations and high seismicity or many stations and low seismicity. In both cases, the seismic network should be better assessed in the latter areas. The local correlations with historical seismicity are negative in parts of the LCW and LCE regions and in minor areas within the UCW, LNE, and S regions. In contrast, the UCE, ENE, and eastern S regions appear to be the most reliable.

In summary, the Italian National Seismic Network shows generally good coverage over the whole territory in consideration of the seismicity and mapped seismogenic sources. The network should be re-assessed (and, when possible, opportunely refined) primarily in the SE and western S regions and secondarily in the NW and UCW regions.

We also cross-checked the distribution of the strong motion (accelerometric) network with the distributions of the seismic hazard and the population. The local correlations with the hazard map indicate a generally positive correspondence with local negative correlations (Fig. 10b). From this map, it is not possible to unequivocally identify whether the negative correlations are due to the lack of stations in hazardous areas or an abundance of stations in low-risk areas; however, some conclusions can be derived. As a general rule, the greater the expected ground acceleration, the higher the number of strong motion stations in a given area should be. Upon classifying the study area according to the ground acceleration classes established by the Italian regulation (Table 2), an increasing proportion between the percentage of stations and the percentage of the study area is actually observed (Table 2). The distribution of the population is extremely uneven: peaks of people are crowded into very restricted areas, while wide areas include few people; therefore, the local correlations did not provide clear results. For this reason, we computed the cumulative number of people around each node of the strong motion network. The results indicate that in the range of 0–100 km, the curves generally follow a similar increasing trend, albeit with a certain dispersion (Fig. 11b). The best-located stations (highest population at a distance of 5 km) are widespread, while those surrounded by few people are mainly located in the C region (Fig. 11d and e). For the most-populated cities (Table 1), the curves readily follow the trend that corresponds to a regular distribution of square cells; the mean curve approximates a regular grid distribution of approximately 27×27 km² (black solid line in Fig. 12). Even though it has different densities, the distribution of strong motion stations is quite regular around the main urban centers, and the corresponding density is higher with respect to the overall density, which corresponds to 1/524 km² ($\sim 23\times$23 km²).

Table 2 Classification of seismic zones according to the Italian regulation (PCM 2006) and based on the ground acceleration expected to be reached or exceeded (a_g) with a 10% probability in 50 years (Fig. 10a). The percentage of accelerometric stations, the percentage of the study area, and their ratio are reported for each class

In summary, the Italian Strong Motion Network has a generally uniform distribution over the whole country. The coverage of the network is in agreement with the seismic hazard distribution, although some local inconsistencies occur (parts of the UCW, UCE, C, and S regions). The network is well calibrated with respect to the distribution of the population.

6 Conclusive remarks

In this paper, we use statistical methods and tools for the description and characterization of networks. The proposed approach aims to evaluate the current state of a monitoring network rather than plan its optimal design. The theoretical design of the geometry of a network plays an important role in its future performance; however, the performance should be validated afterwards through recorded data because local-scale factors could play a role (e.g., the hypocentral depth and crustal velocity structure). These factors cannot be included in our analysis, but they do not present limitations on the interpretation of the results at a larger scale. Indeed, this approach could be a very informative tool to assess the degree of coverage of already developed networks as a function of the spatial distributions of other parameters (including complex and irregular distributions).

For both of the networks presented as case studies herein (i.e., the Italian National Seismic and Strong Motion Networks), the proposed approach has been confirmed to be useful for suggesting directions in planning their future optimization.

• The Italian National Seismic Network has a strongly uneven distribution at first sight because of some areas of dense stations and some other areas with scarce coverage. Considering the distributions of seismogenic sources and seismicity, the station coverage appears to be more consistent, even though some local incongruities occur. These areas have been highlighted in this study; therefore, in the future growth of this network, the station sites could be homogenized to enable a uniform and compavrable response of the network over the whole country.

• The Italian Strong Motion Network has a more uniform distribution than the Italian National Seismic Network and is already in good agreement with both the seismic hazard distribution and the population distribution; however, some areas require a re-assessment. Locally, the network could be undersized according to the seismic hazard and, at the same time, oversized according to the exposed population. The network should be strengthened while following the appropriate balance between all the considered factors.