1 Introduction

A country-scale earthquake monitoring network usually extends over areas in the range of 105–106 km2 (D’Alessandro 2019) and, from its initial planning to its final configuration may be necessary for several years to decades. Throughout this process, the final geometrical arrangement of the nodes could turn remarkably different from the initial, planned configuration. For example, the Italian National Seismic Network (D’Alessandro et al. 2011a), the Red Sísmica Nacional in Spain (D’Alessandro et al. 2013), the Romanian network (D’Alessandro et al. 2012), and also the Global Seismographic Network (GSN) took on average 20–30 years (or even more) to reach their current configurations, and they will keep developing further in the future (Mezcua 1995; Gee and Leith 2011; Popa et al. 2015; Michelini et al. 2016). In some other cases, a network could be the result of the merge between several local sub-networks, as for instance happened in Greece (D’Alessandro et al. 2011b). As a counterpart, in some cases, the long time span necessary to build up the networks turned into a benefit because of the concurrent technological improvements in sensors, transmission systems, etc. In all the cases, it should always be opportune to verify periodically whether the distribution of nodes within a given network copes with its final objectives, mainly because the knowledge about the seismicity of a given area evolves over time.

In this work, we propose the application of the approach proposed by Siino et al. (2020) for the evaluation of the earthquake monitoring network in Taiwan. This approach considers the spatial distribution of the nodes of the network together with the ancillary information related to the aims of the network itself. Taiwan island is characterized by a relevant seismicity (Theunissen et al. 2010; Shin and Teng 2001; Wu et al. 2013), accountable with the active geodynamics of this region which involves active subduction and collision (Barrier and Angelier 1986; Sibuet and Hsu 1997; Rau and Wu 1998; Shyu et al. 2005; Chin et al. 2019). The seismic release is very high: more than one hundred M>6.0 earthquakes are recorded in the last 25 years; the largest earthquake that occurred in recent times is the Mw 7.6, 1999 Chi-Chi earthquake (Shin and Teng 2001; Kim et al. 2010).

The beginning of the instrumental earthquake monitoring in Taiwan dates back to 1897 when the very first seismograph was installed (Shin et al. 2013). Since then, a network has been established and expanded over the whole country, mainly boosted later by the technological advancements in the 1970s and 1990s (i.e., digital signals, real-time data transmission, and telemetry), and it is still expanding and improving at present (Wang 1989; Shin et al. 2013). The map in Fig. 1 shows the main seismotectonic features of the region together with the locations of the earthquake monitoring network in Taiwan. For the abovementioned reasons, the island of Taiwan represents an interesting case study to evaluate the earthquake monitoring networks. Moreover, the last paper dealing with this topic in Taiwan dates back to 1997 (Tsai and Wu 1997). In detail, we considered separately the part of the network which includes the broad-band velocimeters (hereafter Real-time Seismic Monitoring Network ), from the part of the network which includes the accelerometers (hereafter Strong-Motion Earthquake Observation Network). Each of the two specific networks is then analyzed in comparison with various data sets (e.g., seismicity, seismogenic sources, the magnitude of completeness, seismic hazard, and population distribution) according to their mutual relationships as suggested by Siino et al. (2020). Contributing to the assessment of the degree of coverage of the network, this work can possibly address also its future optimizations.

Fig. 1
figure 1

Seismicity distribution in the Taiwan region; both historical and instrumental earthquakes are shown together with the active faults. The colored triangles indicate the locations of the earthquake monitoring networks. The arrow indicates the direction of the plate motion of the Philippine Sea Plate relative to the Eurasian Plate. Elevation data from GEBCO Bathymetric Compilation Group (2020).

Full size image

2 Data

In this section, we introduce the data sets and their sources. The main information concerns the spatial distribution of the earthquake networks which are processed together with the ancillary information related to the objectives.

All the information on the earthquake monitoring network in Taiwan (i.e., location, type of sensor, international code) and about seismicity is available at the website of the Seismological Centre of the Central Weather Bureau (see Data and Resources).

Concerning the observation network, we focused on two different typologies, namely the 24-bit Real-time Seismic Monitoring Network (also known as Central Weather Bureau Seismographic Network, CWBSN) which counts 119 ground stations and 9 OBS at present, and it is in charge of the seismic surveillance of the country and the Strong-Motion Earthquake Observation Network (also known as Taiwan Strong Motion Instrumentation Program network, TSMIP) which counts 420 monitoring sites at present. This last network collects full records of strong earthquakes to understand the seismic activities under different geological conditions. It also serves as the basis for the earthquake intensity assessment and earthquake-resistant design norms which are tools to effectively reduce and prevent earthquake disasters. Technical details about these networks are provided by Tsai and Lee (2005), Hsiao et al. (2009), Chang et al. (2012), Shin et al. (2013), Guan et al. (2020), and Central Weather Bureau (CWB, Taiwan) (2012).

The selected catalogue of instrumental seismicity covers the period from January 1, 1991, to December 31, 2022, and counts 259,290 events, with earthquake magnitude ranging in the interval 2.2–7.3. The minimum magnitude corresponds to the magnitude of completeness (Mc) calculated for the entire catalogue of 787,674 events with a magnitude within the range of 0–7.3. The values of (Mc) are coherent with the ones for inland Taiwan provided by Mignan et al. (2011).

The catalogue for the historical seismicity results from the combination of two different lists, namely the one from Cheng and Yeh (1989) and the one from Chen and Tsai (2008). In the case of multiple magnitude estimation for a single event, we took into account the worst scenario considering the highest magnitude for those events with different values. The merged list counts 33 earthquakes with M> 5 that occurred in the period from 1604 to 1988. Besides, the events with M >= 5.5 from both instrumental and historical earthquakes have been merged into a unique list of 382 events including all those earthquakes usually considered beyond the damage threshold. For these events, we estimated the seismic moment (M0) by converting the provided values of local magnitude by means of the relation proposed by Hanks and Kanamori (1979):

$${M}_0={10}^{\left(M+10.7\right)}\ast 1.5$$

The mapped active faults within the study area are retrieved from the global active faults database (Styron and Pagani 2020). They have been collected from several geological and geophysical studies at various scales. Most of the seismic sources are from Shyu et al. (2016) which are also the same that have been used as the base for the evaluation of the seismic hazard model. The map of the active faults contains 44 elements with lengths ranging from ∼10 to ∼200 km. The main fault systems can be directly framed within the broad-scale geodynamic setting of region 1.

This study also considers the earthquake hazard of Taiwan. Seismic hazard maps provide information about the likelihood of ground shaking in a given area and represent fundamental tools for estimating the likelihood of damage and losses. The more recent Taiwan probabilistic seismic hazard model is the one from Chan et al. (2020) named ”TEM PSHA2020” and hereinafter briefly ”TEM.” The TEM includes the updated earthquake database and seismogenic structure database together with a new set of ground motion prediction equations and site amplification factors. The probabilistic model assesses the likelihood that a given ground acceleration threshold is overcome in a given time span. The TEM model is included in the data analysis considering the probability of exceedance equal to 10% in 50 years, which is the most common reference to display a seismic hazard model.

Finally, we considered the distribution of the population over the study area. Data, in the form of population density (i.e., number of people per map unit), come from official United Nations population estimates and are updated to 2015 (see Data and Resources).

3 Method

In this paper, data are analyzed by means of descriptive spatial statistics and point process methods following the approach proposed by Siino et al. (2020). The methodologies are briefly recalled in the following paragraphs. The analysis is performed with R statistical software (R Development Core Team 2005), and the functions employed are in the packages raster (Hijmans 2018), spatstat (Baddeley and Turner 2005), and spatstat.local (Baddeley 2018).

3.1 Point process methods

A point process model X is a random finite subset of W ⊂ Rd, where |W| < ∞. A spatial point pattern υ = {u1, . . . , un}, which is a realization of 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 1du goes to 0:

$$\lambda \left(\boldsymbol{u}\right)=\underset{\left|d\boldsymbol{u}\right|\to 0}{\lim}\frac{E\left(N\left(d\boldsymbol{u}\right)\right)}{\left|d\boldsymbol{u}\right|}$$
(1)

where E is the expected number of points. The intensity (λ(u)) is assumed to be inhomogeneous, so it is not constant over the considered 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}\left(\textbf{u}\right)=\frac{1}{e\left(\textbf{u}\right)}{}_{i=1}{}^nk\left(\textbf{u}-{\textbf{u}}_i,\boldsymbol{h}\right)$$
(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 = (hx, hy). The bandwidth controls the degree of smoothing, and a larger set of values gives more smoothing. The estimation of the bandwidth vector h = (hx, hy) is computed with Scott’s rule (Scott 1992):

$$\begin{array}{cccc}{{\mathrm h}_{\mathrm x}=\mathrm n^{-1/6}}^\checkmark&\overline{\mathrm{var}\left(\mathrm X\right)}&{{\mathrm h}_{\mathrm y}=\mathrm n^{-1/6}}^\checkmark&\overline{\mathrm{var}\left(\mathrm Y\right)}\end{array}$$
(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 linear form is assumed following Baddeley et al. (2015):

$$\lambda \left(\textbf{u}\right)={\beta}_0+{\beta}_1D\left(\textbf{u}\right)$$
(4)

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

Moreover, we locally infer the spatially varying estimates of these two parameters of the inhomogeneous Poisson process model in Eq. (4). This technique is called “Geographically Weighted Regression” (GWR). GWR is a data exploration technique that allows understanding changes in the importance of different variables over space which may indicate that the model used is mis-specified and can be improved (Fotheringham et al. 1998). Regression models are typically “global”: that is, all data are used simultaneously to fit a single model. In some cases, it can make sense to fit more flexible “local” models. Such models exist in a general regression framework (e.g., generalized additive models), where “local” refers to the values of the predictor values. In a spatial context, local refers to a location. Rather than fitting a single regression model, it is possible to fit several models, one for each (out of possibly very many) location.

3.2 Descriptive spatial statistics

In the second step of the analysis, the intensities of the network are compared with the pertinent information to assess their coherence. The estimated intensities of the point patterns, the hazard maps, and the population distribution are considered raster data sets. Pairs of raster 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 resolution of 3 km and a neighborhood of 5×5 cells (15×15 km2) around each cell. This resolution has been calibrated according to the case study. In addition, as proposed by Siino et al. (2020), we compute some descriptive statistics to infer further insights about the distribution of the earthquake networks with respect to the related data sets.

4 Results

The results are presented considering separately the Real-time Seismic Monitoring Network and the Strong-Motion Earthquake Observation Network. Each of the two networks is then analyzed in comparison with the various data set described in the previous paragraph. All the maps resulting from the analysis are projected according to the Mercator projection UTM51 in km, and all the gridded maps have pixel size resolution equal to 3×3 km. Results are displayed over an area of 300 km by 450 km. Of course, being the stations mainly located in the mainland area, here results are more consistent. In the eastern off-shore where seismicity is consistent and where OBSs are placed, results are robust in any case. In the NW corner of the study area, results should be considered not reliable. First, we calculated the inhomogeneous intensities with a kernel estimator as shown in Eq. (2), given the spatial point patterns of the two networks. Having fixed the grid resolution to 3 km, the values reported in the intensity plots indicate the number of stations over a 9 km2 area.

4.1 Real-time seismic monitoring network

The intensity maps of the real-time seismic monitoring network, briefly the “seismic” network is shown in Fig. 2. The intensity appears to be almost constant over the mainland area: relatively higher intensity values mark the northern portion of the country and part of the western. This seismicity intensity map shows very high values of intensity along a portion of the eastern coast (Fig. 2). The relationship (i.e., local correlation) between intensity maps of seismic stations and instrumental seismicity is evaluated considering the pairs of values around each cell for both maps within a square area of 5×5 cells; the correlation coefficient is then recorded at the central cell of the focal area. Both zones with positive and negative correlation emerge (Fig. 2). The positive correlation indicates coherence between the two considered variables; conversely, the negative correlation indicates discordance between the distributions of the stations and the instrumental seismicity. At this stage, it is not possible to discriminate whether the negative correlations areas are due to the lack of stations in a seismic area or poor seismicity in a well-covered area; however, some considerations will be discussed later.

Fig. 2
figure 2

Left: Kernel intensity estimation for the seismic monitoring network. Center: Kernel intensity estimation for the instrumental seismicity. Right: Local correlation coefficients of the two previous intensity maps. Bandwidth selected with Scott’s rule; pixel size resolution is 3×3 km

Full size image

Furthermore, the cumulative number of events with increasing distance is calculated around each station (Fig. 3). The euclidean distance between each pair of station-epicenter has been considered. We selected the 10% of stations with the highest number of events at a distance of 5 km (green triangles) and the 10% stations with the lowest number of events at a distance of 100 km (red triangles). These threshold values of distance and percentage are taken from Siino et al. (2020) and have the purpose to distinguish which stations are best and worst located 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. 3.

Fig. 3
figure 3

Top: cumulative number of instrumental earthquakes located at increasing distance from each seismic station and corresponding intensity. Bottom: Top 10% of all stations with the highest cumulative number of earthquakes at a distance of 5 km (green triangles) and the top 10% of the stations with the lowest cumulative number of earthquakes at a distance of 100 km (red triangles); seismic stations classified according to the cumulative number of earthquakes at a distance of 50 km

Full size image

To further investigate the relationship between the seismic network and the seismicity, we also mapped the magnitude of completeness (Mc) over the study area. Mc is defined as the minimum earthquake magnitude above which the earthquakes are reliably detected. The completeness magnitude has been mapped for the time interval from 1991/01/01 to 2022/12/31 which includes approximately 790,000 events with magnitude ranging from 0 to 7.3. The estimation of Mc usually resorts to the assessment of the for Guternber-Richter (GR) law (Gutenberg and Richter 1944) for the earthquake catalogue. In particular, assuming that the GR law establishes a linear relationship between the magnitude M and the logarithm of the number of earthquakes that have a magnitude greater than M in a given area, it is possible to estimate Mc as the lower threshold of the linear behavior. Based on this assumption, there are many techniques for the calculation of the completeness magnitude (Wiemer and Wyss 2000; Cao and Gao 2002; Wiemer and Wyss 2000; Woessner and Wiemer 2005; Amorese 2007; Taroni 2021; Godano et al. 2022). In this study, Mc has been assessed by computing the maximum curvature of frequency-magnitude distribution (Wiemer and Wyss 2000) using the ZMAP software (Wiemer 2001). To map the distribution of Mc over the area, we fixed a minimum number of 100 events within a constant radius equal to 50 km around each map unit. The Mc map has been finally smoothed with a moving average window (Fig. 4, center). Mc ranges from ∼1.8–2.0 in the mainland, decreasing outwards, reaching values greater than ∼3.5 in the peripheral areas. This pattern is consistent with the Mc map provided by Mignan et al. (2011). The spatial relationship between the nodes of the seismic network and the active faults has been also 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 Fig. 5.

Fig. 4
figure 4

Left: Kernel intensity estimation for the seismic monitoring network. Center: Smoothed map of the completeness magnitude (Mc). Right: Local correlation coefficients of the two previous intensity maps. Bandwidth selected with Scott’s rule; pixel size resolution is 3×3 km

Full size image
Fig. 5
figure 5

Left: cumulative density function of the distance to the nearest fault of each seismic station. Right: cumulative density function of the distance to the nearest seismic station of each fault

Full size image

Among all the seismic stations, approximately 55% are located within 10 km from a fault, and approximately 90% are less than 25 km away from the nearest fault (Fig. 5). In total, about 70% of the faults are closely monitored by at least one seismic station, and 90% of the faults are monitored from a distance shorter than 10 km. (Fig. 5).

As explained in the previous section, the seismic station intensity is also computed as a function of an external covariate, namely the distance to the nearest fault. The main goal is to assess whether the stations are spatially arranged to occur more frequently near mapped faults. The geological information (i.e., red lines in Fig. 1) is transformed into a spatial variable defined at all locations u ∈ W, namely, D(u), which is the distance to the nearest active fault (Fig. 6). We expected that, increasing the distance to the nearest fault, the station intensity decreases (i.e., negative coefficient of the geographically weighted regression). However, positive values of the slope coefficient mark large portions of the study area, especially in the mainland, indicating, at first glance, an unexpected increase in the number of stations at increasing distance from the faults’ positions.

Fig. 6
figure 6

Left: distances (in km) to the nearest fault for each grid point. Center: spatial distribution of β0 (i.e., intercept) of the geographically weighed regression model. Right: spatial distribution of β1 (i.e., slope) of the geographically weighed regression model; contour line (dashed gray line) corresponds to β1 = 0

Full size image

4.2 Strong-Motion Earthquake Observation Network

The intensity maps of the Strong-Motion Earthquake Observation Network, briefly “strong-motion” network is shown in Fig. 7. The intensity is clearly higher in the western and northern parts of the country, while it is sensibly lower in the inner and southern parts. The intensity map of the M > 5.5 earthquakes shows a gradient moving from west to east (Fig. 7). The relationship (i.e., local correlation) between intensity maps of strong-motion stations and potentially damaging seismicity has been calculated considering the pairs of values around each cell in both maps. This focal square area is equal to 5×5 cells, and we recorded the correlation value at the central cell of the focal area (Fig. 2, right). The positive correlation indicates coherence between the two considered variables; conversely, the negative correlation indicates discordance between the distributions of the stations and the historical seismicity. At this stage, it is not possible to discriminate whether the negative correlations areas are due to the lack of stations in a historically seismic area or, on the contrary, to the poor seismicity in a well-covered area; however, some considerations will be discussed later.

Fig. 7
figure 7

Left: Kernel intensity estimation for the strong motion network. Center: Kernel intensity estimation for the events with M >= 5.5. Right: Local correlation coefficients of the two previous intensity maps. Bandwidth selected with Scott’s rule; pixel size resolution is 3×3 km

Full size image

Similarly to the instrumental seismicity, the cumulative seismic moment M0 with increasing distance is computed around each node of the strong-motion network. The obtained curves provide information about the coherence between the location of each station and the historical seismic release (Fig. 8). Among all the stations, we selected 5% having the highest values of M0 at 5 km (green triangles) and those having the lowest values at 100 km (red triangles). Again, such thresholds are taken from Siino et al. (2020) and have the purpose to recognize the best and worst–placed strong-motion stations with respect to the historical seismic release. A map of all the stations ranked according to the cumulative M0 at 50 km is also presented (Fig. 8), providing a comprehensive explanatory framework for the effectiveness of the strong-motion network. The relationships between the spatial distribution of strong motion stations and the hazard maps and population have been also evaluated. In fact, because one of the goals of a strong motion network is to assess the earthquake intensity (i.e., the effects at the surface), this network should be also planned according to the distribution of people to the expected ground motion. The local correlation between the strong-motion station intensity and the population density (Fig. 9a, b) shows a high spatial variability of the area with positive and negative correlation with short wavelength. Moreover, we have also calculated the cumulative number of people living around each station with increasing distance (Fig. 10). By selecting the 5% of stations with the highest number of people at a distance of 5 km (green triangles) and the 5% stations with the lowest number of people at a distance of 100 km (red triangles), it is possible to distinguish which stations are placed 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 presented (Fig. 10).

Fig. 8
figure 8

Top: cumulative M0 of events M >= 5.5 located at increasing distance from each strong motion station and corresponding intensity. Bottom: Top 5% of all the strong-motion stations with the highest cumulative M0 at a distance of 5 km (green triangles) and top 5% of the stations with the lowest cumulative M0 at a distance of 100 km (red triangles); strong-motion stations classified according to the cumulative M0 at a distance of 50 km

Full size image
Fig. 9
figure 9

Top row: population density (displayed as log10 of the number of people per map unit) and local correlation with the strong motion network. Bottom row: expected PGA from the TEM model considering 10% exceedance probability in 50 years (Chan et al. 2020) and local correlation with the strong motion network. The pixel size resolution is 3×3 km

Full size image
Fig. 10
figure 10

Top: cumulative population at increasing distance from each strong-motion station and corresponding intensity. Bottom: Top 5% of all the strong-motion stations with the highest cumulative population at a distance of 5 km (green triangles) and top 5% of the stations with the lowest cumulative population at a distance of 100 km (red triangles); strong-motion stations classified according to the cumulative population at a distance of 50 km

Full size image

Finally, the local correlation between the strong-motion station intensity and the expected peak ground acceleration in 50 years has been computed (Fig. 9). The results indicate a general positive correlation, whereas small-scale regions of negative correlation occur especially in the northern half of the island.

5 Discussion

The assessment of the quality of earthquake monitoring networks is necessary to recognize both weak and strong points and to address the continuous development. We performed this task for the seismic and strong motion networks in Taiwan following the approach proposed by Siino et al. (2020) which apply point-process analysis together with descriptive spatial statistics.

The characteristics of the seismic network have been assessed taking into account the instrumental seismicity, the Mc, and the active faults. Results of the local correlation between the station intensity and the instrumental seismicity highlight zones with negative correlation (red areas in Fig. 2). Considering the almost uniform station intensity (Fig. 2), most of these negative areas are likely the result of a significant network coverage in zones with limited instrumental seismicity. The positions of the stations compared with the distributions of the earthquakes indicate that the best-placed stations (i.e., the highest number of events at a distance of 5 km) have a good match with the intensity of instrumental seismicity (Fig. 3). The less coherent stations (i.e., the lowest number of events at a distance of 100 km) are located in the peripheral areas and some of them are the OBS stations. In particular, the latter sensors are essential to better locate and constrain the high seismic area located in the eastern offshore of the main island 3. Considering the cumulative number of earthquake at increasing distance, only a few stations have a low gradient whereas most of them have a notable, even though variable, gradient (Fig. 3). For some of them, more than 150,000 events have been located at a distance shorter than 100 km.

The correlation between the station intensity and the Mc indicates a general positive correlation in the inner and northernmost parts of the island, and a general negative correlation in the outer ones (Fig. 4). The density of the station and Mc is expected to be uncorrelated because the distance from stations is among the controlling factors in reducing the Mc (c.f. Mignan et al. (2011)). For this reason, the areas with positive correlation would represent anomalies (Fig. 4, right).

We found a good coherence also checking the distribution of the seismic network nodes against the distribution of active faults. This is first confirmed by the cumulative density functions (Fig. 5), while the results of the geographically weighted regression would indicate an apparent disagreement. In fact, in some areas, the number of stations increases at increasing distance from faults (Fig. 6). This can be framed within the almost uniform station distribution (Fig. 2) so that many stations are located in areas where no faults are mapped. In a future perspective, this configuration can be considered a strong point because the fault database cannot be considered definitively complete, and other faults will be mapped in the future also in the areas where they are not mapped at present.

The characteristics of the strong motion network have been assessed taking into account the distribution of the strongest events (M >= 5.5) and their associated seismic release (M0), the population, and the seismic hazard model. The results of the local correlation between the intensities of stations and of the M >= 5.5 events highlight large portions with negative correlation (red zones in Fig. 7, right). In particular, the large red area in the inner part of the island is likely the result of a relative lack of strong motion stations in an area experiencing a relatively abundant strong seismicity; as a counterpart, this area corresponds to a low-populated area (Fig. 8). Conversely, the smaller area along the upper-eastern coast is likely the result of the high density of M >= 5.5 earthquakes.

The stations with the higher gradient and greater cumulative M0 at a 50 km distance (Fig. 8) are placed coherently with the expected values of ground motion (Fig. 9). The worst-placed stations according to the lowest M0 at a distance of 100 km are exactly located where the expected ground motion is lower and where few people live 9.

The distribution of the population is highly uneven, being characterized by small and highly populated zones, and large areas with few inhabitants (Fig. 9 top-left). The resulting local correlation shows small-scale variations whose interpretation is not straightforward. The gradient of the cumulative population at increasing distance highlights three well-defined groups of stations (Fig. 10 top) clearly recognizable in the map (Fig. 10 bottom). The best-located strong motion stations with respect to the population are clustered in a very restricted area in correspondence with the capital city (green triangles in Fig. 10), while the worst-placed stations are located along the lower eastern coast and in some islands (red triangles in Fig. 10). The local correlation between the strong motion stations and the seismic hazard model (TEM) results in general positive values (Fig. 9), with the negative values grouped in small portions clustered mainly in the northern part of the study area. This is also the area in which the seismic hazard is generally lower (Fig. 9) and uncorrelated areas are the effect of a relatively greater number of stations.

6 Conclusive remarks

The geometry of a network plays a role in its performances. In particular, earthquake monitoring networks should be validated according to their final objectives in consideration of associated information like the earthquakes’ and faults’ distributions, as well as seismic hazard model and populations.

The planned design of the geometry of a given network is crucial in its future performance; however, its assessment should be validated afterwards through recorded data because local-scale factors could play a role (e.g., the hypocentral depth and crustal velocity structure) (Scudero et al. 2021). Although these factors have not been included in our analysis, they do not present limitations on the interpretation of the results at a larger scale (Siino et al. 2020). Indeed, this approach represents a very useful tool to assess the degree of coverage of monitoring networks as a function of the spatial distributions of correlated parameters (including complex and irregular distributions), and it has been applied for the characterization of the earthquake monitoring networks in Taiwan.

For both of the Taiwan networks (i.e., seismic and strong motion), the adopted approach revealed useful to assess the present-day quality of the earthquake monitoring and to suggest directions in planning their future optimization. The main considerations achieved with this work can be summarized as follows:

  • The shape of the island inevitably influences the geometry of the earthquake networks and this, in turn, clearly influences many results showing NNE-SSE elongated features. However, at a large scale, does not emerge any evidence of critical weak areas for both the seismic and strong-motion network in Taiwan. They both are coherent with their own objectives of seismic surveillance and intensity/hazard assessment, respectively.

  • The seismic network extends even over zones apparently not justifiable at present because of the lack of mapped faults. This might turn out in a strong point in a future perspective if further studies will reveal previously unmapped faults, as revealed in the experiment by Wu et al. (2016). The OBSs deployed in the eastern offshore area reveal essential to locate the relevant seismicity of this sector. However, all the stations are essential to better locate and constrain all those earthquakes occurring in the more external parts of the network. The hazardous area along the eastern coast is well covered, and the weaker parts of the networks in the inner part of the island, all things considered, are low seismic zones. Another advantage of such a high-density network is enabling the scientific research like seismic tomographies and other seismological studies with great coverage and resolution (Wu et al. 2007, 2008; Sun et al. 2015; Kuo et al. 2016). Similarly, the strong-motion network well covers zones where the current seismic hazard model indicates moderate hazard, without overlooking the most hazardous zones. Moreover, especially for the eastern coast, there is a relevant simultaneous coherence among all the information: high expected ground motion, low population, and appropriately placed strong motion stations.

  • For both networks, only restricted portions of the country show a limited coherence between the networks with the associated information, but not with multiple parameters simultaneously.