Full three-dimensional relocation and tomographic inversion of the 1977-2008 earthquakes in north-eastern Italy: a feasibility study

In a regional seismological network, the estimation of the epicenter is usually robust, especially for events inside or close to the network boundaries. In contrast, the hypocentral depth is very sensitive to the assumed velocity field. In this study, we compare the hypocenter estimates obtained by a classical algorithm in a simple one-dimensional (1D) model with a recently developed full 3D model that is based on shrinking grids. This study is preliminary, as the 3D Earth model is based on limited data from the literature; however, it demonstrates that different patterns show up when a more representative geological model is adopted. This encourages further studies, based on fully integrated 3D models from active surface seismic, well data and other geophysical measurements. Such an integrated approach has been successfully adopted by the oil and gas industries for decades, which has increased the exploration success rate and the production of hydrocarbon reservoirs.


Introduction
Seismic tomography has been used by seismologists during the past century over extended regions with very sparse networks.Their main applications, i.e., crustal studies and civil protection, do not require high resolution, as errors of 1 km in an epicenter location are acceptable.A different scenario has become apparent over the last few decades, when seismic tomography was adopted first by exploration seismologists [Bishop et al. 1985, Zhang andToksöz 1998, among others] and later by reservoir geophysicists [Laurent et al. 1993, Zhou et al. 1993, Laurent et al. 1995, Zhou et al. 1997, Wittrisch and Maisons 1998, Maxwell and Urbancic 2001, Rossi et al. 2001a, Rossi et al. 2001b, Vesnaver et al. 2003, Fortier et al. 2005, among others] to improve the threedimensional (3D) velocity models of hydrocarbon reservoirs.The resolution needed for reservoir monitoring is higher than for seismological cases: the bin size for a 3D surface seismic survey is about 30 m × 30 m in the x and y directions.As a survey can reach a size of 20 km × 20 km, the resulting spatial resolution is about 0.15%.This value approximates the seismological case if a 1-km error occurs for a network sized 600 km × 600 km, or so; thus, the relative resolution can be comparable.The key difference is the reliability of the 3D velocity model: it can approach the bin size for a surveyed reservoir, but it is often poor at a seismological scale.Only a few countries (e.g., The Netherlands) are fully covered by 3D seismic surveys.In the USA, the related data acquisition is still ongoing within the USArray Project (www.usarray.org).
In this study, we investigate the full 3D relocation accuracy for earthquakes recorded from 1977 to 2008 in northeastern Italy (Friuli) by the network of the Italian National Institute of Oceanography and Applied Geophysics (OGS).Using the actual coordinates of the recording stations, we first simulate the travel-times in a 3D synthetic model, which approximates the major geological structures in the Friuli region.As a second step, we locate the hypocenters and compare these with the modeled ones.In this way, we can measure the relocation error as a function of the adopted velocity model and the algorithm parameters.Another goal was to compare the accuracy between the classical Hypo71 code, used by the OGS for cursory preliminary analysis, and the recent Hypo3D code, which was developed for reservoir monitoring of hydrocarbon production [Vesnaver et al. 2008[Vesnaver et al. , 2009[Vesnaver et al. , 2010]].The Hypo3D code can jointly invert the traveltimes from active and passive seismic surveys, with receivers located both at the Earth surface and in boreholes, using direct, reflected, refracted and converted waves.

Algorithm overview
The adopted algorithm Hypo3D was described in detail by Vesnaver et al. [2008Vesnaver et al. [ , 2009Vesnaver et al. [ , 2010]].Here, we briefly review a few key features that explain part of the relocation differ-ences we can get when adopting other popular tools, such as VELEST [Kissling et al. 1994], HypoEllipse [Lahr 1999], Hy-poInverse [Klein 2002], NonLinLoc [Lomax et al. 2000] or the classical Hypo71 [Lee and Lahr 1972].This last was used for the standard relocations reported in the OGS catalogs.Other methods have been presented more recently, using cross-correlation and spectral properties [Shearer 1997, Du et al. 2004], flexible gridding [Thurber and Eberhart-Phillips 1999], and the double-difference algorithms [Waldhauser and Ellsworth 2002], among others, which are not discussed further here.The NonLinLoc algorithm provides a spatial probability for hypocenters, and so it is not directly comparable with the other methods mentioned above.These last methods estimate the hypocenter coordinates and the origin time by minimizing the misfit between the measured arrival times with those simulated for a given Earth model at all of the stations.
The philosophy of the Hypo3D code we used is to split as much as possible the time from the space unknowns; i.e., the time origin from the hypocentral coordinates.The velocity field and the origin time are kept constant during the relocation of one or several events, and it can be updated before or after a relocation phase is completed.The relocation is carried out with the 'shrinking grids' method, which comprises the following steps: (i) A first guess for the hypocentral coordinates (x 0 , y 0 , z 0 ) is made, as described above.
(ii) Centered at the current hypocenter guess, a 3D regular grid is defined with points at (x 0 ±i Dx, y 0 ±i Dy, z 0 ±i Dz), where the space intervals Dx, Dy and Dz are chosen arbitrarily (and are not necessarily equal), as are the number of points for each dimension.
(iii) The difference between the measured and modeled travel-times is computed at each grid point, and summed into the object function we want to minimize (the norm of this difference can be linear or squared, but the results obtained are very similar in most cases).
(iv) The grid point where this function is a minimum becomes our new guess for the hypocenter.
(v) If the object function is lower than a user-defined threshold, if its decrease-rate approximates to zero, or if we reach a maximum iteration number, we stop the procedure and the current hypocenter guess becomes our final estimate.
(vi) Otherwise, we reduce the grid intervals Dx, Dy and Dz by a user-defined shrinking factor, and then return to step 2.
A first major difference between the Hypo3D algorithm [Vesnaver et al. 2010] and Hypo71 (or its evolutions, as Hy-poInverse and HypoEllipse, or VELEST) is the minimization procedure.The Hypo71 algorithm simultaneously estimates the spatial hypocentral coordinates (x H , y H , z H ) and the time origin t 0 based on the method of Geiger [1912].The object function to be minimized depends on station coordinates, the picked arrival times, and the quadri-vector (x, y, z, t), which is composed of the hypocentral coordinates and the time origin.These last variables are handled as mutually independent.This hypothesis appears reasonable from a geometrical point of view: why should we process the spatial coordinates in a different way, or even the time origin?
There are indeed good reasons for this.Seismological stations are mostly located at the Earth surface, so they correctly span the x and y coordinates, and much less the vertical direction (the opposite happens in microseismic surveys for hydrocarbon reservoir monitoring, where receivers are mostly located in vertical wells, so spanning the z coordinate much better than the other two dimensions).Another reason is the limited thickness of the Earth crust, and so of the hypocenter depth, which are normally smaller than the horizontal extension of regional seismological networks.Thus, both sources and receivers are sampling the x and y direction better than the z direction.For this reason, the Hypo3D method estimates the epicenter coordinates (x H , y H ) first, because these unknowns are much better constrained by data than the depth and time origin.The initial position for the epicenter are the x and y coordinates of the station with the minimum P travel-time; the depth is obtained by tracing a vertical ray from there, according to the assumed Earth model, such that the resulting travel-time is equal to the observed travel-time.
The method of Wadati [1933a, b] can estimate the time origin t 0 by a cross-plot of P versus (S-P) arrival times.This estimate is reliable when its basic assumption holds; i.e., that the ratio between the average P and S velocities is constant along the ray paths of the picked events.In this case, the estimate is independent of the interface shape and velocity changes in three dimensions.Therefore, mixing such a wellconstrained parameter with others that are not so well constrained, such as the hypocentral depth, is a weakness [e.g., see Bishop et al. 1985].In the Hypo3D algorithm, the Wadati [1933a, b] method is used as a first pass, and later marginal corrections are allowed with respect to the initial value according to a user-defined threshold, to allow for cases when the Wadati [1933a, b] assumption is violated.Normally our knowledge of the P and S velocity field in three-dimensions is not precise, especially at a crustal scale; if so, we can only guess at the cases when the Wadati [1933a, b] method can fail; i.e., when the regional structure of the Earth crust cannot be approximated by plane horizontal layers with a constant ratio of P and S velocities.Instead, in reservoir monitoring for hydrocarbon production, where detailed information about these velocities is available, such an assessment can be accurate, and can be quantified by seismic modeling.
In the quadri-vector (x H , y H , z H , t 0 ) that identifies an event, the hypocentral depth z H is heavily dependent on our 3D velocity model for the P and S waves.When the time origin t 0 is known, so are the travel-times from the hypocenter to the recording stations.The hypocentral depth depends almost linearly on the average velocities along the ray paths; thus, it is interpreter-biased or data-dependent as much as our velocity model is.
Another key difference between the Hypo71 and Hypo3D algorithms is how they account for elevation.Hypo71 assigns a station delay, redatuming the arrival times to a plane, in a similar way that static corrections are applied to seismic records in exploration surveys.Instead, Hypo3D computes the travel-times in a full 3D model, where the elevation of the recording station is used directly, without any redatuming.As a result, even shallow ray-bending effects or azimuthal changes of velocities are allowed for.
Other algorithms have been presented over the last decades that are not considered further here, such as the 'simulPS' program [Thurber 1993] and the 'tomoDD' program [Waldhauser and Ellsworth 2000, Waldhauser 2001, Zhang and Thurber 2003].The VELEST code [Kissling et al. 1994] has similar features to Hypo71, as it simultaneously estimates the hypocenter and time origin in a 1D model.Such a simplification allows simultaneous inversion of the model velocity to be tried too, which can introduce relevant ambiguities when it comes to a more realistic Earth model.A possible comparison is left for future studies.However, we can get a good idea about the reflection tomography from the classical study of Bishop et al. [1985].They showed that there is severe cross-talk between velocity and depth errors, even in the particularly favorable case of active seismic surveys, where we know not only the receiver positions (as in earthquake seismology), but also the location, timing and radiation pattern of the sources.

Two models
Representing the Earth by a 1D model, where the velocities depend on depth alone, is a crude simplification of the Hypo71 algorithm, which was dictated by the modest computing resources available when this code was written [Lee and Lahr 1972].The errors so introduced depend on the difference between the 1D and the 'true' 3D model, and thus they changing from case to case.Eberhart-Phillips and Michael [1993], among others, showed the value of a full 3D approach with pioneering applications to a Californian region.To quantify the possible improvement in a case of practical interest, we built a simple 3D model for north-eastern Italy, based on the lithological and geological information available in the literature.We understand that such a model is certainly wrong in terms of the details, but hopefully it is not so wrong in terms of the regional trends; e.g., having higher velocities in the Alps and lower velocities in the Friuli Plain.Such an obvious difference, however, is not allowed in the Hypo71 algorithm, and is partly mitigated by its successors, HypoInverse and HypoEllipse.
The Friuli region is one of the most seismically active regions in Italy, with relevant seismic hazard posed by infrequent earthquakes with magnitudes >5.5.The seismic ac-tivity is driven by the anti-clockwise rotation of the Adriatic microplate within the Africa-Eurasia collision zone, with under-thrusting and wedge indentation of the Adriatic microplate with respect to the Africa-Eurasia collision zone.Bressan et al. [1998] reported the regional seismicity as shallow and confined to the upper 30 km of crustal depth.For details of the tectonic evolution, we recommend Galadini et al. [2005] to the reader.In our feasibility study, we assumed simple point sources, and focused our attention on the kinematic aspects.Point-source approximation acts as a lower limit for the solution precision: the smaller is the rupture area, the higher the achievable resolution of the relocation grid.
Figure 1 shows a simple layered model for the area, which is often used for a preliminary cursory relocation of earthquakes.However, it does not include the subduction zone, with its strong heterogeneities.In the 3D layered model we built (Figure 2), the curved dipping interfaces mimic the subduction trend.The model area is roughly square, extending from ca. 45°30' N to 46°40' N, and 12° E to 14° E, spanning a distance of 127 km in the northerly direction and 230 km in the easterly direction.These boundaries were chosen to include the whole of the OGS network, although without embracing too large an area.The vertical boundaries are the Earth surface for the top, and the Moho discontinuity for the bottom, which is located at a variable depth of between 35 km in the southern part, down to 60 km  below the Alpine zone.
The layers are homogeneous, with the P and S velocities obtained from petrophysics data [Faccenda et al. 2007] and four well logs (in Cargnacco, Bernadia, Lavriano and Terenziano).The two uppermost layers define the Alpine area, and their thickness vanishes in the plain zone.The third and fourth layers represent the formations that fill the Friuli Plain and subduct under the Alps in the north, below the bottom of the seismogenetic zone.The files of these models can be downloaded from the website www.earth-prints.org,together with other numerical results presented in this study.
Figures 1 and 2 allow the comparison of the 1D and 3D models we used to build the comparisons.The simpler model is used for the standard relocation, and is composed of two horizontal layers over a half-space; however, as its only interface is quite deep (at 22 km), most of the events are located in the upper layer.The 3D model is composed of five layers with smooth interfaces.In both cases, the layers were modeled as homogeneous, with the P and S velocities re-  1. One-dimensional model parameters.The P and S velocities, and their ratios, as illustrated in Figure 1.
Table 2. Three-dimensional model parameters.P and S velocities, and their ratios, as illustrated in Figure 2. ported in Tables 1 and 2. We note that the Vp/Vs ratio is not constant among the different formations in the 3D model.
Figure 3 shows the geological model along the Transalp Profile, as proposed by Cassinis [2006], and based on previous work by Scarascia and Cassinis [1997] and on new data from Lueschen et al. [2005]; the dotted lines in Figure 3 denote the limits of the area covered by our model.This model is relatively similar to a N-S section of our model close to the Transalp Profile (Figure 4), at a distance of 29.9 km from our model origin.

Relocation comparisons
The relocation of hypocenters depends greatly both on the velocity model and on the inversion algorithm.To single out these factors, we carried out a few tests, where we changed one of these factors at a time.As the hypocenter location of real earthquakes cannot be checked experimentally (e.g., by drilling, as in seismic exploration), we modeled synthetic events as point sources, in such a way that the hypocenter coordinates are perfectly known, and compared these with the estimates provided by the different methods.
To mimic a realistic distribution of hypocenters, those estimated from the real data and published in the OGS bulletins from the years 1977 to 2008 were assumed to be exact, and were used for our modeling exercise.Published locations are mainly calculated on a 1D model with the Hypo71 algorithm.Only events with small magnitude (M L <3.5) are considered, to maintain the validity of the point-source approximation.Synthetic travel-times for all of the available stations in the existing network (Figure 5) were computed, keeping the hypocentral coordinates fixed, but with the two different models described in the previous section.Table 3 reports the absolute average root mean square error in the time residuals, hypocenter coordinates, and time origin ob-tained by the relocation of these events.The iteration number is 10 for all cases, and the other relocation parameters are the same too: the shrinking grid resolution is 5 × 5 × 9, and its initial size is 10 × 10 × 8 km, which implies an initial grid sampling of 2 × 2 × 0.9 km.The minimum final grid interval depends on the estimation convergence for each of the combinations of models and travel-times [Vesnaver et al. 2010], although in most cases this is close to 25 × 25 × 25 m, a value that approaches the expected source size for the magnitude range considered.
In the 10 top rows of Table 3, the travel-times and models have a consistent dimensionality: 1D with one dimension, and 3D with three dimensions.The input data are the traveltimes in the 1D or 3D models in Tables 1 and 2, respectively, although the estimations are carried out also by models where the velocities are 10% and 20% lower or higher with respect to these values.We note first that the time origin is always exact when dealing with 1D models and 1D-modeled travel-times.In these cases, the hypothesis of the Wadati [1933a, b] method holds, and the precision obtained is amazing.The average error is very low in the full 3D case too, as 2 ms approximates the typical sampling rate of the traveltimes.As expected, the precision is superb in the simplest case; i.e., for the 1D model with its correct travel-times there are average errors of 1 m in the epicenters and 5 m in the hypocenter depth, which might appear ridiculously good to seismologists, but will appear appropriate to petroleum engineers for the monitoring of the microseismicity in a producing reservoir.
Both the P and S velocities are perturbed by the same factor, so their ratio does not change, nor does their compliance with the Wadati [1933a, b] hypothesis change.Indeed, the error in the estimated origin time depends only on the model dimensionality, while it is constant for different ve-   Table 3. Absolute average root mean square time and coordinate misfit for the different combinations of 1D and 3D travel-times and relocation models.The 'Time misfit' column shows the average squared differences between the 'known' (synthetic) and the estimated travel-times, the 'XY misfit' and 'Z misfit' are the average squared differences between the actual and estimated hypocentral coordinates, and the 'T 0 misfit' is the different between actual and estimated time origin.locity fields.The two cases have different complexities of the misfit function to be minimized, while using the same number of data.The 1D case allows a global minimum to be reached quickly, while for the 3D case, there are probably several local minima that can prevent the convergence of the iterative minimum search.
Figure 6 shows the tight correlation of the errors in the hypocentral coordinates with the velocity model perturbation.In the 1D case (Figure 6, left), a quasi-linear dependence shows up: a 10% error causes a 4-km error in both the epicentral and hypocentral coordinates.The curve of time misfit is very flat, but its minimum is correctly located at the unperturbed model position.The 3D case (Figure 6, right) is quite different.First, the minimum time misfit is located at the wrong model, where the velocities are 10% higher than those used to model the travel-times.The epicentral estimate is instead robust: even following the criterion of the minimum time misfit, it is lower than the 1D case, except at the plot origin.The major problem shows up with the hypocenter depth, the trend of which is relatively flat, and all values exceed 5 km.This implies that the available data for this model complexity and parameter choice cannot accurately estimate the hypocentral depth.At the same time, we note that this is not due to the chosen algorithm, as in the simpler model, the accuracy achieved is outstanding.
The mismatch between the adopted model and the real Earth structure is another factor that limits the relocation accuracy.The last two rows in the Table 3 provide us with a flavor of the types of errors that we can expect in this area.We exchanged models and travel-times, relocating in 1D model events with travel-times computed in the 3D model, and vice versa.The errors are significant, but not huge, and again, they are about double for the depth with respect to the epicentral error.This means that the 1D model is a good approximation for the more complex 3D model, which mimics the real geology to some extent.Therefore, this approach is viable for a real-time alert system, when the computational speed is critical and when such precision is acceptable.It is not viable, instead, for engineering applications and in oil and gas production.
The results in Table 3 are directly comparable, as they are obtained using the same parameters.One of these parameters is the number of iterations for the relocation process.Very marginal improvements can be obtained by increasing this beyond the chosen value of 10. Figure 7 shows the time misfit as a function of the iteration number for the 3D model.We see that after six iterations, the convergence becomes extremely slow.Therefore, the 10 iterations adopted for the computing of the values in Table 3 are appropriate for practical applications.
The performances in Table 3 only give an idea of the  upper limits that can be approximated under ideal conditions, which are never met in the real world: point sources, high frequency signals, no noise or picking errors, readability of all S arrivals, and so on.Therefore, these quantify the errors due only to the model choice, the origin time estimation, and the relocation algorithm.The actual errors might be larger.On the other hand, these 'numerical errors' can be reduced by increasing the resolution of the shrinking grid, thus approaching the extreme case of an exhaustive search.The higher this resolution, the lower is the probability of being trapped in local minima, but the higher is the computational cost.Local minima in a 3D model are harder to determine, as they do not depend simply on the relative locations of the receivers, as in the 1D velocity model, but also on the velocity model in use.Nonetheless, this problem can be solved, particularly in reservoir monitoring, by recording data according to the correct design of the receiver net-work.The chosen values for the shrinking grid resolution are a reasonable compromise, which reach (and surpass) the physical resolution limits in jobs running for about 10 h on a laptop computer and 7 h on a Linux workstation.As this problem is highly parallel, we expect that this estimation would run on a small PC cluster within 1 h.

Relocation and tomographic inversion of real events in north-eastern Italy
Figures 8 to 10 allow the comparison of the epicenter distributions obtained with the same data, but using Hypo71 in the 1D model (Figure 1, Table 1) and Hypo3D in the 3D model (Figure 2, Table 2).In Figure 8, we highlight a linear trend in the upper left corner and a group of events close to S. Donà di Piave, which were obtained by the 1D model.In Figure 9, these patterns change into an alignment with a different dip in the upper part, and to a similar alignment in the     Differences are even more apparent when the hypocentral depths estimated by Hypo71 (Figure 11) and Hypo3D (Figure 12) are compared over the years.In both cases, during the early years of the network operations (from 1977 to 1990), most of the events are estimated as shallow.In that period, fewer stations were available and the receivers were analog, so the data quality was lower than in the following years.However, we note that in Figure 9 there is some horizontal alignment that might be due to processing artifacts, as Hypo71 estimates can depend on the initial guess for the hypocentral depth.Such alignments are barely visible in Figure 12, where instead the general trend suggests a decreasing depth over the years for most of the hypocenters.Changes in receiver quality and in their spatial distribution can have a strong influence, although they are not easy to evaluate in the temporal pattern, due to the creation of local minima in the vertically homogeneous velocity model.For this reason, the two dashed orange lines in Figure 12 provide qualitative trends only.Any more rigorous analysis will require assimilation of the hypocenter estimates, while compensating for different network geometries and recording instruments, which is beyond the scope of the present study.
Keeping the hypocenter coordinates fixed, we inverted the travel-times for the P and S arrivals using the 3D relocation model (Figure 2, Table 2) as our initial guess for the tomography.Figures 13 and 14 show the changes obtained in this way for the P velocities.In the first layer (Figure 13), a low-velocity anomaly shows up across the borders between Italy and Slovenia.In the second layer (Figure 14), instead, a higher velocity is obtained in the Alps.While this second result is quite obvious, the low-velocity anomaly de-   Figure 15 highlights once more the sensitivity of the hypocentral depths to the chosen velocity model.The three sequences show the evolution over the years of the average hypocentral depths using Hypo71 and a 1D model (Figure 15, blue dots), and using Hypo3D with full 3D models; i.e., the initial (Figure 15, yellow dots) and the final (Figure 15, orange dots) model of the tomographic inversion.Their trends over time show relatively close correlation, but their average levels are very different, and these depend on the different average velocities in each model, as a function of depth.We note that the differences among the three trends vary across the years, and that they do not represent a simple time shift or scale factor.This will be due to the non-linearity of the relocation process, especially when a 3D model and a tomographic inversion are involved.
Figure 16 is a histogram plot with the hypocentral depths estimated with the initial tomographic model and the final tomographic model.This reveals an interesting detail: not only are most of the events of the pre-tomography model (Figure 16, cyan bars) shallower than those of the post-tomography model (Figure 16, grey bars), but a kind of bi-modal distribution appears: a major peak with shallow events at 5 km, and a second deeper event at 12 km in depth.If this can be confirmed in further studies using more accurate velocity models for the Earth, this phenomenon might open new discussions about the geodynamics of this area.

Computational aspects
The inversion algorithm we used exploits the sparsity of the tomographic matrices and the algebraic reconstruction technique or the simultaneous iterative reconstruction technique (see van der Sluis and van der Vorst [1987], among others).In the figures presented here, the simultaneous iterative reconstruction technique was used, because its estimate does not depend on the order of the travel-times in the input.The approach adopted has two main computational advantages.First, it is very parsimonious in terms of the required computer memory, which is linearly proportional to the model parameters.For other approaches that require the inversion of the corresponding fully populated matrix, such a function might be quadratic.Secondly, it provides a very finegrain parallelism: the travel-time computation along each ray can be carried out independently by a different CPU, and the  number of possible travel-times is almost without limit.The computational speed-up in a parallel implementation is nearly linear as a function of the available CPU numbers.Vesnaver et al. [2010] presented an example where active and passive seismic data were jointly inverted: 3.5 million traveltimes took ca.12 hours to be inverted on a small cluster of eight CPUs.We are aware that jobs with over 260 million travel-times required a few days on a larger and faster cluster of 64 CPUs [Vesnaver 2008].Much larger jobs are needed nowadays, as seismic acquisition systems with 200,000 channels are used in the field, and new ones with 1 million channels are being introduced [Poggiagliolmi et al. 2012].As a full 3D seismic survey might include a hundred thousand shots, the number of travel-times to be inverted might exceed 1 billion.Without a massively parallel implementation, the computational challenge appears to be impossible.Furthermore, these methods that require the inversion of large, non-sparse matrices, the singular-value decomposition, or the minimization of object functions by the conjugate gradient method might become prohibitive despite the significant computational power mostly available nowadays.
In the case of the Friuli data, the tomographic imaging is very rapid, taking less than a minute on a single CPU, as it only involves a few thousand events.The grid resolution was chosen to be coarse because the sparsity of the available data does not guarantee a reliable solution without significant ambiguities, unless arbitrary damping factors are introduced.In this way, instead, the number of events per domain is large enough to average out, to some extent, miss-picked or misslocated events.
The relocation code Hypo3D uses the same ray-tracing routines as the tomographic part, and requires a similar iteration number to converge; i.e., <10 in most cases.Only a scalar version has been developed so far, but a performance similar to the tomographic part can be expected, because the computational granularity is the same.

Conclusions
The new algorithm for full 3D relocation and traveltime inversion of seismological events provides estimates for hypocentral coordinates that are different from those provided by the Hypo71 code.These differences are not large for the epicenters, so this can be used for an initial cursory location and alert for civil protection purposes.However, the depth of the seismogenic layers allows the prediction of the lateral extension and the local intensity of possible damage at the surface, and therefore it is still provides relevant information that we might need to be accurate and available in a short time.When carrying out studies relating to seismogenesis and tectonic movements, the resolution improvements that can be achieved are significant.In a simple 1D exercise, the resolution achieved is superior to the possible physical resolution of real experiments; in three dimensions, it is com-parable laterally, although still poor vertically, with average errors of almost 6 km.
The 1D model we analyzed has been used to relocate the hypocenters published in the OGS bulletin.This provides a good equivalent model to the more complex 3D model for the Friuli area, as it achieves lower errors in both the time misfit and the hypocentral coordinates.Its simplicity improves the convergence of the relocation algorithm towards the local minima that are closer to the optimal ones.The complexity of the 3D model can probably be better controlled by a larger number of stations.
This feasibility study shows that the hypocentral estimates in this area mainly depend on the velocity model: the average depths range from 5 km to 15 km, according to the adopted model.These models highlight different features, both in the spatial patterns and in the time evolution of the hypocenter coordinates, with consequently different geodynamic implications.The key way to remove these ambiguities is by adding new geophysical and geological data, and building a solid and more detailed 3D Earth model, since a model with horizontal homogeneous layers hardly provides information about local structures and lateral heterogeneities.The newly developed Hypo3D code has relevant potential applications in terms of the resolution and accuracy of the Earth model, as it can jointly invert any existing and future active and passive seismic data.Such an integrated approach has been used in the last decades by the oil and gas industries for their so-called reservoir characterization, which allows rock properties, such as density and elastic parameters (which can be estimated from seismic data), to be linked to other properties, such as azimuthal anisotropy (which is mostly related to fault patterns), intrinsic anelastic absorption (which is mostly related to fractures), and apparent anelastic absorption (which is mostly related to fine layering and scattering) [Liner 2012, among others].
Figure 1.A simple 1D model composed of two homogeneous layers overlying a homogeneous half-space.White dots, hypocenters estimated using the Hypo71 algorithm.

Figure 2 .
Figure 2. A smoothed 3D model for the Friuli area composed of five homogeneous layers with curved interfaces.White dots, hypocenters estimated using the Hypo3D algorithm.

Figure 4 .
Figure 4. Vertical and horizontal slice across the model in Figure 2 (left), and vertical section in the N-S direction (right) at a distance of 29.9 km from the model origin.

Figure 5 .
Figure 5. Map of the OGS seismological network in north-eastern Italy.Red rectangle, modeled area in Figures 1 and 2.

Figure 6 .
Figure 6.Time and hypocentral misfits for a 1D model (left) and a 3D model (right).The vertical axis indicates meters for the misfits in space (green and red lines), and milliseconds for the misfit in time (blue line).

Figure 7 .
Figure 7. Time misfit for a 3D relocation as a function of the iteration number.

Figure 8 .
Figure 8. Epicenters estimated using the Hypo71 algorithm.The orange lines highlight patterns to be compared with the corresponding ones in Figure 9.

Figure 9 .
Figure 9. Epicenters estimated using the Hypo3D algorithm.The orange lines highlight patterns to be compared with the corresponding ones in Figure 8.

Figure 11 .
Figure 11.Time evolution of the hypocentral depths according to the Hypo71 algorithm.Orange dashed line shows a flat, constant trend.

Figure 12 .
Figure 12.Time evolution of the hypocentral depths according to the Hypo3D algorithm.Orange dashed line shows a dipping, increasing trend.

Figure 13 .
Figure 13.Tomographic inversion of the first layer in the model in Figure 2. Homogeneous initial model (left) and estimated P velocities (right).

Figure 14 .
Figure 14.Tomographic inversion of the second layer in the model in Figure 2. Homogeneous initial model (left) and estimated P velocities (right).
serves further attention from a geological and geodynamic perspective.
Figure 15.Time evolution of the average hypocenter depth using different algorithms and models.Blue dots, Hypo71 in a 1D model; yellow dots, Hypo3D in the 3D model in Figure 2; orange dots, Hypo3D in the 3D model obtained by the tomographic inversion displayed in Figures 9 and 10.

Figure 16 .
Figure16.Histograms for the estimated hypocenter depths using the model in Figure2(blue) and the model with tomographic velocities (gray).