In the following we describe and analyze, using synthetic examples, the two main types of errors affecting the localization process, namely (i) the uncertainty of arrivaltime estimation (picking error) and (ii) the imperfect knowledge of the compressional and shear velocity models. We assume that the nonlinear nature of the localization problem is properly accounted for by the inversion algorithm (localization software) thus implying a negligible contribution to the whole error. The geometry of the monitoring network, obviously, influences the error propagation in the localization process.
2.1 Uncertainty of arrivaltime estimation
There is a quite vast body of literature discussing the factors influencing the identification accuracy of the P and S phases on a seismogram (e.g., Husen and Hardebeck, 2010; Abakumov et al. 2020). The consensus is that there are subjective factors affecting the picking, such as the experience of the operator performing it, as well objective ones such as the strength of the seismic event and dominant frequency.
Assuming that the subjective factors can be at least limited, the signaltonoise ratio (SNR) becomes the main factor controlling the quality and the accuracy of the picking. Hence, the improvement of the quality of the picking is directly connected to the preprocessing applied to the seismograms with the aim of improving their SNR.
Picking errors define the precision of the hypocenter location since they lead to a relative location scattered around the true earthquake point location (Husen and Hardebeck, 2010). The pick is commonly assumed to be a random variable symmetrically and normally distributed around the true arrival time, even though some authors have demonstrated that the distribution can be highly asymmetrical (Diehl et al. 2012, Abakumov et al. 2020).
Following these considerations, we assume that the estimation error follows a Gaussian distribution with zero mean and variance \({\sigma }^{2}\). This latter can be quantified using the information theory and the theorem of ShannonHartley (Aki, 1976), and described by the equation
\({\sigma }^{2} = {\left({\text{log}}_{2}\left(1+\frac{SNR}{K}\right)\cdot {2\cdot f}_{max}\right)}^{2}\)  (1) 
where \({f}_{max}\) is the maximum frequency of the recorded signal, SNR the signaltonoise ratio, and \(K\) represents an empirical factor equal to \(K=20\) (Aki, 1976). The derived picking errors, for different values of the parameter \({f}_{max}\), are shown in Fig. 1 as a function of the SNR. It is worth noting that for good quality events (i.e., SNR > 10) the picking error curves become progressively flat with increasing frequency and ultimately nearly independent of \({f}_{max}\). In practical applications, fmax can be taken as the maximum of the simulated spectra of the seismic event (using, for example, Brune’s 1976 model) whereas SNR is a design parameter. It is important to note that both the seismic spectrum and the SNR are also used in the computation of network sensitivity, hence a link between picking uncertainty and event detectability can be established. In other words, the picking error distribution computed from Eq. (1) is directly associated to the minimum magnitude that the seismic network will be able to detect for a given SNR.
2.2 Estimation of the traveltime modelling errors
The errors in modelling of the traveltimes of P and S phases arise from the inexact knowledge of the compressional and shear waves velocities. This uncertainty is even more pronounced at locations where a new monitoring network must be designed and installed, because the knowledge of the subsurface conditions is generally rather limited at that stage of development.
Seismic velocity can be derived from a number of information sources, ranging from very localized, high frequency, seismic logs in boreholes to estimates derived from distant seismic events and recorded by the local network. Active seismic acquisition methods (based on refraction, reflection, or surface waves) can also be utilized and cover the intermediate scale range. The variety of scale, frequency range, acquisition geometry and nature of the propagating waves (P, S, Rayleigh) pose a severe challenge in integrating the data into a single, potentially threedimensional, model of P and S wave seismic velocities that is also, potentially, anisotropic. This is particularly true for the depth range between the scope of nearsurface active acquisitions (~ 100 m) and deep investigations relevant to the reservoir zone (~ 1000 m). Whatever the intricacies of this data integration process, a reliable P and S velocity model is essential to provide reliable localizations of seismic events, and possibly not enough attention is traditionally paid to this aspect that proves absolutely critical in the case of smallscale, highresolution monitoring networks focused on a specific area of interest.
Although it is feasible to use threedimensional velocity models for localizing earthquakes, the use of layered models (i.e., onedimensional velocity models as a function of depth) is a very common and a widely accepted practice in both global and local seismology (e.g. Husen et al., 2011; Theunissen et al. 2018), in particular during the network design phase. For local networks, this approach is valid when lateral variations of elastic properties are negligible. This may depend on the local geological situation. In addition, it is also useful to recall that during the design of a new seismic monitoring network a 3D velocity model at the scale of interest is seldom or never available, therefore this situation has clear practical implications.
Considering that the velocities of each layer are affected by a certain degree of uncertainty, to exemplify how the uncertainty in the velocity model affects traveltimes in the following we will assume that both P and Swave velocities are distributed symmetrically and uniformly in the range of ± 10% about the estimated layer velocities.
Based on the assumption above of a uniformly distributed velocity error of ± 10%, \({v}_{i}\) (the velocity of the \(i\)th layer) is a random variable uniformly distributed between \(0.9{v}_{i}\) and \(1.1{v}_{i}\), while \({p}_{i}\), the slowness of each layer, is a random variable with inverse uniform distribution. Consequently, traveltime, defined as the integral of slowness along the raypath, is distributed as the weighted sum of inverse uniform random variables. The weights involved in the computation of the traveltime depend on the lengths of the raypaths in each layer which in turns depend on the position of the hypocenter as schematically shown in Fig. 2a for two different hypocentral locations. The position of the hypocenter with respect to the station (the black triangle in Fig. 2a) also controls the number of random variables to be summed.
As the velocity, and consequently the slowness, changes from layer to layer, the corresponding random variables span different ranges and thus the probabilistic distribution of the traveltime errors cannot, in general, be derived analytically. The effect of the velocity errors on traveltime can be evaluated numerically by a stochastic approach. For the sake of an example, we generated an ensemble of 1000 1D velocity models (500 P and 500 Swave velocity models) by perturbating the velocity value of each layer in a range of \(\pm 10\%\) using a uniform distribution. Figure 3b shows the average onedimensional P and S wave velocity models, red and blue curve respectively, with their corresponding velocity ranges. The modelling is performed for 24 hypocenter locations per depth slice (at depths of 2, 3, 4 and 5 Km, hence a total of 96 hypocenter locations) on a regular coarse grid with 2 km spacing in both \(x\) (Easting) and \(y\) (Northing) directions, ranging from 2 to 12 km Easting and from 2 to 8 km Northing (see Fig. 5a as reference). Exploiting the onedimensional nature of the velocity models, thus in a cylindrical coordinate system, the epicenters locate at distances ranging from 1.2 to 12 Km for a reference station (see Fig. 4d). A schematic view of the hypocenter positions with respect to the reference station is shown in Fig. 2b.
Figure 3a and 3c show examples of the perturbation distribution for a single hypocenter position and for Svelocity and Pvelocity, respectively. The histograms of the velocity values in each layer are juxtaposed vertically, with the shallowest layer at the top of the display. It can be noted that at deeper layers the velocity distributions are no longer uniform and are no longer symmetrical about the true velocity values. This is due to a limitation of the used raytracing method that only accepts velocity models that increase monotonically with depth. At deeper layers (e.g., below 5 km in Fig. 3b) the differential velocity between two adjacent layers decreases. This leads to a narrower and asymmetric interval about the mean value where random values are drawn.
For each hypocenter and for each perturbed model in the ensemble, we computed the traveltime using the forward model of Hypoinverse (Klein, 2002). Each run is then compared against the reference traveltime provided by the unperturbed model. This yields an ensemble of traveltime realizations that shows the uncertainty associated with traveltimes for the reference station and for the different hypocentral distances and depths considered.
Figure 4a and Fig. 4b show the histograms of the traveltime errors for a set of hypocenters at 2km and 5km depth, respectively. The Gaussian distribution overlaid to the histograms has mean and standard deviation computed from the population. Figure 4c displays the histogram of the traveltime errors, overlaid by the Gaussian curve estimated from the entire model ensemble, for all the hypocenters. Figure 4d shows, for each modeled hypocenter, the relative error associated to the traveltime (ratio between the absolute traveltime error (\({\Delta }t\)) and the reference time \({t}_{ref}\)) expressed as a percentage. The points are plotted in the distancedepth plane together with a map obtained via interpolation of the points. Interestingly the percentage ratio is between 2% and 3% and shows a clear decreasing trend with depth. Figure 4 considers only the effect on Ptraveltime errors, but a similar behavior is observed using S velocities.
2.3 Localization uncertainty assessment
Based on the considerations and the results presented in the previous sections, here we evaluate the performances of a monitoring network by computing the spatial uncertainty of the localization. To this end we consider a real monitoring network installed by STOGIT in Italy. The network is composed of 7 surface stations installed with the spatial configuration shown in Fig. 5a. Figure 5b shows the corresponding azimuthal coverage (AG, Azimuthal Gap) for this network.
We estimate the earthquake hypocentral locations on a fine grid in both \(x\) (Easting) and \(y\) (Northing) directions (with 150 m spacing in both directions) at 4 different depths: 2, 3, 4 and 5 Km. The resulting grid has more than 24000 hypocenter locations.
The traveltime for each of these hypocenters is perturbed, for each realization \(r\), according to the equation
$${\left({T}_{P,S\left(pert\right)}^{i,j,k, s}\right)}_{r} ={\left({T}_{P,S \left(true\right)}^{i,j,k,s}+{\epsilon }_{pick}+{\epsilon }_{tt}\right)}_{r}$$
2
where \({T}_{P,S \left(true\right)}^{i,j,k,s}\) is the true traveltime (for both P and S waves) from the hypocenter at the coordinates \(i,j,k\) along the axes \(x, y, z\), respectively, to the station \(s,\) \({\epsilon }_{pick}\) is the picking error, and \({\epsilon }_{tt}\)is the error caused by the uncertainty associated to the velocity models. In this work the total number or realizations is \(r=500\).
The picking error \({\epsilon }_{pick}\) is calculated from Eq. (1) with parameters \({f}_{max}=25\)Hz, and SNR=5 and SNR=2 for shear and compressional wave traveltimes, respectively. These SNR are realistic for a situation when the picking is manually performed by an analyst.
The traveltime modelling errors \({\epsilon }_{tt}\) are estimated, according to the results obtained in the previous section, assuming a relative error \({\Delta }\text{t}/{t}_{ref}\) constant and equal to 3%. For each hypocenter location, each traveltime realization is inverted using the software package Hypoinverse (Klein, 2002).
To prove how the location accuracy is sensitive to the Azimuthal Gap (AG), we show the analysis of the inverted hypocenters for all model realizations in two localization cases. The first case comprises the localization of seismic events with a low AG (approximately \({115}^{\circ }\)), whereas the second scenario shows the localization of events with a high AG (approximately \({175}^{\circ }\)).
The first case, displayed in Fig. 6, shows the inverted hypocenters at 4 different depths (2 Km, 3 Km, 4 Km, ad 5 Km) but with the same epicentral location and with an AG of approximately \({115}^{\circ }\). It can be observed that the inverted hypocenters have a low dispersion (a low variance) about the true hypocenter locations both in the horizontal and in the vertical directions and picture with a high accuracy the true hypocenter location. Also note that the variance of the scattered inverted hypocentral locations increases with depth.
The second case (Fig. 7), although the events’ epicenters are near one seismic station (ST6), has an AG approximately equal to \({175}^{\circ }\), that is close to the 180° threshold as indicated in the literature for a bad network performance (Bondár, et al, 2004; Tiira et. al, 2016). In fact, the localization of the hypocenters has a wider distribution both in the vertical and horizontal directions, i.e., the variance of the inverted hypothetical increases as function of the AG.
Extending the localization analysis to all epicentral locations on the fine grid and for all depths, we have an ensemble of 500 locations for each hypocenter location. Each of these ensembles is then used to derive statistical measures such as:

the horizontal localization error, given by the \({L}_{2}\)norm of the true and median of the estimated epicenters, which will be indicated as \({\delta }_{h}\).

the vertical localization error, given by the \({L}_{1}\)norm of the true and median values of the estimated hypocenter depths, indicated as\({\delta }_{z}\)

the standard deviation of the horizontal localization error, indicated as \({\sigma }_{h}\) and

the standard deviation of the vertical localization error, indicated as \({\sigma }_{z}\).
Since these quantities are evaluated on a fine grid, it is convenient to display them in form of maps as done in Fig. 8–11, where each panel represents one of the aforementioned quantities for each given hypocentral depth.
The horizontal localization error at all depth slices, as shown in Fig. 8, consistently follows the AG values for this network configuration. Therefore, the minimum error is found inside the seismic station network where the AG has lower values (see Fig. 5b). Although the error shape map remains the same at different depths, on average it increases with depth.
Inside the network the error that we make in the localization of epicenters is smaller than 200 m and increases up to values greater than 600 m at the edges of the model. This is observed for all hypocentral depths.
Differently from the horizontal localization error, the vertical error (Fig. 9) shows a different pattern with respect to the hypocentral depths.
To a hypocentral depth of 2 km the error shows lower values inside the seismic network with absolute values smaller than 150 m. Conversely, for deeper hypocenters the minimum values of the error progressively move towards the peripheral regions of the model. At 3 km depth the localization uncertainty shows a narrow region of minima that surrounds the seismic network, while at 4 km and 5 km the region of minimum error moves progressively away from the seismic network, showing consequently an increase of the localization uncertainties inside it. In fact, the uncertainty associated with a seismic event at 5 km depth, with an epicenter falling inside the seismic network, reaches values greater than 400 m (Fig. 9d).
The maps showing the statistical dispersion of the inverted hypocenters at all evaluated locations, quantified with the measure of the standard deviation, for the horizontal and vertical directions are displayed in Fig. 10 and Fig. 11, respectively. In all cases, for both vertical and horizontal errors and for all evaluated depths, the standard deviation increases from the center of the network to the peripheral regions of the model. This evidence agrees with the AG map where it reaches values greater than 180° outside the seismic network.