A first look at the Gargano ( southern Italy ) seismicity as seen by the local scale OTRIONS seismic network

On April 2013, a local scale seismic network, named OTRIONS, composed of twelve short period (1 Hz) three component seismometers, has been located in the northern part of the Apulia (southern Italy). In the first two months of data acquisition, the network recorded about one hundred very small (ML<2) magnitude earthquakes. A three-layer 1D VP velocity model was preliminarily computed, using the recordings of earthquakes occurred in the area in the period 2006-2012 and recorded by the national seismic network of INGV (Istituto Nazionale di Geofisica e Vulcanologia). This model was calibrated by means of a multi-scale approach, based on a global search of the minimum misfit between observed and theoretical travel times. At each step of the inversion, a grid-search technique was implemented to infer the elastic properties of the layers, by using HYPO71 to compute the forward models. In a further step, we used P and S travel times of both INGV and OTRIONS events to infer a minimum 1D VP velocity model, using a classical linearized inversion approach. Owing to the relatively small number of data and poor coverage of the area, in the inversion procedure, the VP/VS ratio was fixed to 1.82, as inferred from a modified Wadati diagram. The final 1D velocity model was obtained by averaging the inversion results arising from nine different initial velocity models. The inferred VP velocity model shows a gradual increase of P wave velocity with increasing the depth. The model is well constrained by data until to a depth of about 25-30 km.


Introduction
Despite the Gargano promontory is a part of the Adria foreland (Figure 1), it is characterized by an unusual seismicity rate, comparable with that of seismically active part of the Italian peninsula [e.g., Di Bucci and Angeloni 2013].However, this part of Italy has been monitored less than other areas, probably owing to the smaller magnitude of the events.Instrumental observations are in fact available only for low to moderate size events of moment magnitude M W not greater than 5.7 [Del Gaudio et al. 2007], although historical documentation reports cases of catastrophic events which killed people in the order of thousands [e.g., Patacca and Scandone 2004].Moreover, historical catalogues report at least eleven events having an estimated M w >5.5 in the last millennium [Gruppo di Lavoro "Mappa della Pericolosità Sismica" 2004].The present day seismicity seems to be related to tectonic activity along the approximately E-W Mattinata fault and adjacent faults [e.g., Del Gaudio et al. 2007] (Figure 1).
Based on these evidences, in the frame of an European Territorial Cooperation Programme Greece-Italy 2007-2013 (acronym OTRIONS, INTEREG III), on April 2013 a seismic network was installed on the Gargano promontory (Figure 2a).The geometry and instrumental properties of the OTRIONS network will be presented in a next section.The positions of the recording sites were chosen to cover the part of the Gargano promontory that is affected by the higher rate of seismicity, according to the actual knowledge of the Italian seismicity as reported in the INGV seismic catalogue.The sites were selected among those made available by local institutions controlled by the Puglia Region government.The installation of this regional seismic network was aimed at improving the knowledge of the seismogenic potential of the area, through either the geometrical and dynamical characterization of the active faults and the elastic and inelastic properties of the crustal rocks.
In this study we focus on the inference of a velocity model for the upper crust lying below the Gargano promontory.The inference of a reliable velocity model is the first step in the assessment of the seismic hazard of a region, allowing to locate the recorded seismic events, to compute their magnitude and to image the geometry of the active faults.
We carried out two different types of inversion of travel times using two different data sets.A first threelayer model was obtained before the installation of the OTRIONS network.To this aim, we used the dataset of the events recorded in the area during the period 2006-2012.We re-picked all P and S phases of 220 events (2<M L <3.5) recorded in the area by the national seismic network managed by INGV.Owing to the peculiar configuration of the INGV network (Figure 2b), that has only two stations steadily operating in the Gargano area, the ray sampling of the crust is not optimal to infer a 1D velocity model from a standard linearized inversion of travel times.For this reason, we used a multi-scale approach, based on the calculation of several thousand forward layered models, by progressively increasing the number of layers of the crust.The obtained 1D velocity model was then implemented in the SEISCOMP3 software [Olivieri and Clinton 2012] to allow the automatic and preliminary localization of the events recorded by the OTRI-ONS network.
In a second inversion, we combined P and S travel times of both INGV and OTRIONS events to infer a 1D velocity model of P-waves from the inversion of P and S travel times.The problem can be formulated as a linearized inverse problem around initial values of unknown model parameters [Kissling et al. 1994].The data are represented by the travel times of first arrival P and S phases of a set of earthquakes recorded at an array of seismic stations.The unknown model parameters are represented by the hypocenter location and the origin time of the earthquakes, the V P and/or V S velocities of a one-dimensional layered medium and the station delays.
Many papers [e.g., Husen et al. 2011, Matrullo et al. 2013] have focused on the difficulties of finding a stable solution to this nonlinear inverse problem.In fact, this problem does not admit an unique solution and, as an effect of the linearization, the final model may depend on the starting model.The most used approach to overcome this problem is due to Kissling et al. [1994] that proposed to compute the so-called "average minimum 1D model" as the average of a suite of inverted models arising from different initial velocity models that account for the information content coming from geology.The VELEST code [Kissling 1995] was used to compute the minimum 1D V P model.

Geological, geodynamical and structural setting of the Gargano promontory
The Gargano promontory is a part of the Adriatic plate and forms an isolated massif, distinct from the Apennines by both structural and morphological setting, with elevations of more than 1000 m above sea level, still not involved in the accretion of central-southern Apennines [e.g., Mostardini and Merlini 1986].The Adriatic plate is formed by continental lithosphere and is subducting towards west below the Apennine chain; it represents a promontory formed by the collision of Africa and Eurasia plates [Channel et al. 1979].The Adriatic plate principally extends beneath the Adriatic sea [Anderson and Jackson 1987], albeit it is exposed in southern Italy in the Apulia region.In the complex geodynamic context of the area, the Adriatic plate is considered as the foreland of both the Apennines and the southern Alps, at west, and of both the Dinarides and Albanides thrust belts, at east (Figure 1) [Di Bucci and Angeloni 2013].
The Apulia region mainly consists of three deformed carbonatic plateaus (Gargano, Murge, Salento), separated by transversal morpho-structural depressions, which allowed the accommodation of differential vertical movements.These plateaus are cutted by normal fault systems of different orientation and age (from Mesozoic to Pleistocene) [Del Gaudio et al. 2007].The Apulia foreland shows a rather uniform structure, with a Variscan crystalline basement and an approximately 4-6 km thick Mesozoic sedimentary cover and it is discordantly overlain by thin, discontinuous Late Pliocene-Pleistocene deposits [Funiciello et al. 1991, Bosellini et al. 1993].The Moho is located at a depth of about 30-35 km [Piana Agostinetti and Amato 2009].
Geological and geophysical data indicate that the Gargano is a region of local crustal uplift and anomalous contractional deformation within the relatively less deformed Apulian foreland [Brankman and Aydin 2004].In particular, gravity surveys revealed a positive Bouguer gravity anomaly of 110 mGal, coincident with the location of the Gargano uplift [Finetti and Morelli 1973].Detailed information on the properties of the carbonate platform have been inferred by four deep drilling performed by the Italian Company for Oil Exploration (AGIP) [Mostardini and Merlini 1986, Bosel-lini et al. 1993, 2000, Improta et al. 2000].These data have been recently reanalyzed by Festa et al. [2013] that proposed a model consisting of four seismostratigraphic units, with V P ranging from 1400 m/s to 6400 m/s, whose thicknesses varies from a well to another one and may indicate a strong heterogeneity of the crust at a local scale.
From a structural point of view, the Gargano promontory is characterized by a widespread brittle deformation with normal and strike-slip faults [Brankman and Aydin 2004].A major active E-W striking shear zone, known as Molise-Gondola shear zone, cuts a foreland zone characterized by an unexpectedly high level of seismicity [Di Bucci and Angeloni 2013].This shear zone is formed by well-known faults, as the Mattinata fault that cross-cuts the southern part of the Gargano promontory.It is possible to identify three main fault trends: NW-SE, ENE-WSW, E-W (Figure 2a).One of the most prominent NW-SE fault is the Apricena normal fault, that extends about 30 km and was recognized by Patacca and Scandone [2004] (Figure 2a).The most important E-W faults are represented by the Mattinata and the Tremiti faults (the last recognized principally offshore), respectively towards south and towards north of Gargano.These two faults are characterized by a strike-slip kinematic.The Tremiti line is an approximately 50 km long fault and is described as a dextral strike-slip fault [e.g., Argnani et al. 1993].The 60 km Mattinata fault (on shore) is composed of two main active fault segments, dipping at high angle toward the north: the San Marco in Lamis fault to the west and the Monte Sant'Angelo fault to the east, connected through a right-step.

The OTRIONS Seismic Network
On April 24, 2013, the OTRIONS Seismic Network was installed on the Gargano promontory.The OTRI-ONS network is composed of 12 three component seismic stations whose position is shown in Figure 2a.Each station consists of a 24 bit SL06/SARA data-logger (dynamic range equal to 124dB at 100 sps) equipped with a short-period Lennartz 3D-V seismometer (flat response above 1 Hz).The acquisition system allows the recording of data on an external USB device and their real time transfer, through a modem MOXA, to a seismic laboratory, located at the Dipartimento di Scienze della Terra e Geoambientali, Università di Bari "Aldo Moro".The realtime data transfer is realized by using a GPRS/UMTS connection (Figure 3).Data are transferred also to INGV and Regione Puglia centers.The transfer of data is managed by SEED link protocol.A server collects data from the stations by using the software "OnCell Central Manager", that allows to archive data in SEED format.Moreover, these data are sent to a PC where they are managed by the SeisComp3 software.
SeisComp3 is a widespread software aimed at acquiring and exchanging seismic data on the web.It allows to visualize the seismic records in real time [Hanka et al. 2010].SeisComp3 performs also the automatic picking of P waves on the traces, allowing for the automatic location of the seismic events and the computation of their magnitude.Moreover, it allows both the manual reprocessing of P and S phases, the relocation of the events and their storage in a seismic bulletin.The automatic location of seismic events is based on the automatic recognition of P phases and the consequential inversion of travel times.The automatic picking of traces is performed through a STA/LTA algorithm [Withers et al. 1998].A minimum of 6 phases is necessary to detect a seismic event.The automatic localization of the events is performed by using the LoCSAT method [Bratt and Bache 1988] and the IASPEI 91 velocity model [Kennett 1991].The software allows also to locate the events by using the NonLinLoc global search algorithm [Lomax et al. 2000].Moreover it allows to allocate in input an arbitrary Vp velocity model.
In the first month of data acquisition, the detection of the events was supported by five INGV seismic stations.During the first two months of data acquisition, 67 seismic events have been recorded by the OTRIONS seismic network (Table 1).After the re picking of P waves and the picking of S waves with Seis-Comp3, the events were relocated using the LOCSAT method and the IASPEI 91 velocity model.The spatial and temporal location of these events is reported in Table 1. Figure 4 shows the location of all the events.Figure 5 shows the location of the 27 events that were recorded by the OTRIONS network but were not recorded by INGV or EMSC networks.

Data
Before the installation of the OTRIONS network, we analyzed the waveforms of 220 seismic events (2≤M≤ 3.5) recorded by the seismic network of INGV, localized in the Gargano promontory in the period 2006-2012.These data were integrated with the seismic recordings of more than 100 small magnitude events (0.3<M L <1.9) recorded by the OTRIONS network in the period ranging from April 24 to May 31, 2013.The detection of the events was visually performed using the SAC [Goldstein and Snoke 2005] seismic software.The visual analysis of traces allowed us to detect a number of earthquakes greater than those automatically inferred by SeisComp3.We performed the manual picking of P and S travel times on each trace of INGV and OTRIONS seismic events.A weighting factor inversely proportional to the uncertainty associated to the travel time was assigned to each data.Table 2 reports the correspondence between the weighting factor and the corresponding range of error on data.Figure 6 shows the picking of P and S phases on the seismograms of a small magnitude event (M L =0.9) recorded by the OTRIONS network.The histogram representing the quality of the overall dataset is shown in Figure 7.   Table 1 (continued from previous page).Spatial and temporal location of seismic events recorded by the OTRIONS seismic network in the period March 24 -April 30, 2013.
Since a magnitude scale for the Gargano region is not yet available, the magnitude of these events was computed, after the relocation of the events in the 1D V P model, using the equation recommended by the IASPEI: In Equation (1), R is the hypocentral distance in km, typically less than 1000 km and A is maximum trace amplitude (in nm) that is measured on the horizontal com-ponents of waveforms, after the deconvolution for the instrumental response of the available seismometers and the convolution with the response of a Wood-Anderson standard seismograph, but with a static magnification of 1. Equation ( 1) is equivalent to that used by INGV for the calculation of the event magnitude and is an expansion of the Hutton and Boore [1987] formulation.

Layered V P velocity models
A layered V P and V P /V S velocity model for the area, was calibrated by the inversion of the travel times of P and S phases of the seismic events localized in the Gargano region by INGV in the period 2006-2012.This model was then used in the preliminary location of the events recorded by the OTRIONS network.Figure 2b shows the position of the INGV seismic stations that recorded at least one available datum.The length scale of the INGV network (Figure 2),that is of the order of several hundreds kilometers, indicates that the inferred model will only represent the averaged elastic properties of an area that covers several regions of southern Italy.
In our analysis we used a multi-scale approach that consists of progressively increasing the degree of complexity of the crust.First, we inferred the best fit half-space LOOKING AT THE GARGANO SEISMICITY      reduction with respect to the homogeneous model of about 15% and a RMS reduction of 9% was inferred.
The three obtained models are shown in Figure 10.

The 1D velocity model
We used the Velest code [e.g., Kissling et al. 1994] to determine the "minimum" 1D velocity model.This code is based on a damped least square approach and several iterative inversion steps and allows to infer the model parameters through a linearized approach.As an effect of linearization, the inferred parameters (seismic velocities and hypocenter locations) are generally dependent on the starting values.
As concerns the actual dataset of P and S travel times, the main problems are represented by both the relatively small number of available events (320) and the greater source to receiver distance range of the events recorded by INGV network, that has a typical length scale (of the order of hundreds kilometers) greater than the local scale (of the order of ten kilometers) of the OTRIONS network (Figure 2).For this reason, we decided to reduce the number of degrees of freedom of the problem by fixing the V P /V S value and, therefore, the V S 1D profile.For the same reason, we did not consider the effect of the station delays, even because the number of INGV stations that recorded at least one event is higher than 100.
Following Matrullo et al. [2013], a first selection of data was performed by removing from the dataset all P and S travel times having a residual higher than a fixed threshold (1.5 s in this study), after their localization with a simple homogeneous velocity model (V P =5.5 km/s; V P /V S =1.8).Moreover, we included in the dataset only those events that have at least four P travel times and 2 S travel times.The total number of events available for the study reduced to 280 (200 recorded by INGV and 80 recorded by OTRIONS).After the removal of outliers, a total number of 3580 P wave travel times and 1800 S wave travel times was selected.
Before of carrying out the inversion of data, the overall dataset of P and S waves of the seismic events recorded by both the INGV and the OTRIONS networks was used to compute the V P /V S ratio.To this aim, we used the method proposed by Chatelain [1978], that consists of determining the slope of the straight line that best fits the difference between couples of S wave travel times t S,i − t S,j vs. the difference between couples of P wave travel times t P,i − t P,j for each couple (i,j) of stations and for each event.Data are plotted in Figure 11 and are well interpolated by a straight line with V P /V S =1.82, with a linear correlation coefficient R 2 =0.98.
As concerns the inversion of V P velocity profile, based also on several previous studies [Kissling et al. 1994, Husen et al. 2011, Matrullo et al. 2013], the final minimum 1D V P model was obtained as the average of several inverted 1D models, arising from different starting velocity models.The following nine starting LOOKING AT THE GARGANO SEISMICITY   d) The model of Costa et al. [1993], that was inferred after a zoning of the Italian territory.e) A regional velocity model of the Apulian plate [Venisti et al. 2005].
f ) The three-layer velocity model described in the previous section (named hypo71).
Both the two homogeneous velocity models and the three gradient velocity models were chosen on the basis of the available values of the seismic velocities of the upper crust in the Gargano promontory, by taking into account the well known heterogeneity of the crust at a local scale in the area [e.g., Festa et al. 2013].
In each inversion process (i.e., for each starting velocity model) we followed the guidelines prescribed by Kissling et al. [1994], by using different damping coefficients for hypocenter parameters and velocity model.The inversion steps can be briefly summarized as it follows.In a first step, we used a damping coefficient 0.01 for hypocenter parameters and a damping coefficient 0.1 for the velocity model.In this step, we jointly inferred all the parameters several times, by updating the starting velocity model with the computed velocity model.After the relocation of the events, we performed a further inversion step, using a damping of 0.01 for the hypocenters and 1.0 for the velocity model, with the aim of finding the velocity model that minimizes the total estimated location error [e.g., Kissling et al. 1994].We do not allowed the presence of low velocity layers in the inversion.In fact, after several preliminary trials, we deduced that the use of low velocity layers gives rise to unstable solutions, as often described in literature [e.g., Kissling et al. 1994, Matrullo et al. 2013].
The nine inferred velocity model are summarized in Figure 13a.With the exception of the velocity model of Costa et al. [1993], which gave rise to a higher final  RMS (0.43 s), the adherence of model to data is comparable for the remaining eight velocity models (RMS of the order of 0.38 s).Therefore the averaged 1D velocity model was computed using only these eight models (Figure 13b).

Discussion and conclusion
The events were relocated using the average minimum 1D model (Figure 14).The most part of the events occurs in the Gargano area.In this area, the event distribution is rather spread, even if a higher number of events is located between San Giovanni Rotondo, Monte Sant'Angelo and Manfredonia.Therefore, the position of the epicenters confirms that the seismic activity is mainly related to the tectonic activity in the shear zone that comprises the Mattinata Fault and the Apricena Fault and some minor lineaments, as found in previous studies [Di Bucci and Angeloni 2013].
The events tend to concentrate until to a depth of about 30 km (Figure 14).Moreover, in the three-layer velocity model (Figure 10) V P abruptly increases to 7.3 km/s at a depth of about 27-30 km.These two results seem indicate that the Moho is located at a depth of about 27-30 km, as previously inferred in a teleseismic receiver function analysis [Piana Agostinetti and Amato 2009].
Moreover, we note the similarity of V P /V S value inferred from the grid search technique (V P /V S =1.81 for the half-space model in Figure 8) with the value inferred from the Chatelain [1978] method for the whole dataset (V P /V S =1.82 in Figure 11).These values are in close agreement with the results obtained by Piana Agostinetti and Amato [2009] and may indicate that the crust, in the Gargano area, is characterized by a moderate fluid content.If we compare this value with the average V P /V S =1.89 [Chiarabba and Amato 2003] of the near Umbria-Marche Apennine, we conclude that the Gargano promontory is characterized by a minor fluid content, that could be indicative of a minor degree of fracturing of the crust.
The errors on source parameters were computed using two different approaches.First, we estimated the formal errors on horizontal and vertical coordinates of the event foci using Hypo71 (Figure 15a).A further calculation (Figure 15b) was made by con-LOOKING AT THE GARGANO SEISMICITY  ferred three-layer model (30 % of RMS reduction).An average total residual of 0.36 s is inferred in the 1D V P model, that reduces to 0.24 s for the data recorded by the OTRIONS network (Table 3).
The small number of events considered in this study does not allow us to image the geometry of the active faults.This objective will require the further analysis of the about one thousand events further recorded by the OTRIONS network in the period from May 2013 to March 2014.The analysis of these events is still in progress and could help us to better constrain the elastic properties of the crust in a future study.
As final consideration we note that the use of a local scale array in a region of apparently moderate seismicity allows the detection of very small magnitude (minimum M L =0.3) events, that are extremely important to better extend the range of completeness of seismic catalogues and therefore in the evaluation of the seismic hazard of an area.Even if the magnitude of the events has been computed using Equation (1), a further study has to be carried out to calibrate a magnitude relationship for the Gargano area, as it has been done for the southern Italy [Bobbio et al. 2009]

Figure 2 .
Figure 2. (a) Location of seismic stations of the Otrions seismic network.The surface projections of the most important faults are also represented.A.F: Apricena fault; M.F: Mattinata fault; T.F:Tremiti fault; C.F.F: Cerignola-Foggia fault; S.F: Sannicandro Garganico-Apricena fault (redrawn from Del Gaudio et al. [2007]).(b) Location of INGV seismic stations considered in this study.Only stations recording at least one available datum are shown.The size of triangles is proportional to number of T P and T S readings at each station.

Figure 3 .
Figure 3. Architecture of the real time system of transfer, visualization and archiving of data of the OTRIONS seismic network.

Figure 5 .Figure 6 .
Figure 5. Automatic location of 27 seismic events recorded at the OTRIONS seismic event that were not detected by INGV and EMSC.

Figure 7 .
Figure 7. Histogram representing the overall quality of the selected dataset of P-and S-phase arrival times.

Figure 8 .
Figure8.Initial RMS plot in the two-dimensional (V P and V P /V S ) parameter space.

Figure 9 .
Figure 9. Final RMS plot for the search of the best-fit three-layer Vp velocity model.

Figure 10 .
Figure 10.The layered Vp velocity models obtained with the gridsearch method.

Figure 11 .
Figure 11.Modified Wadati diagram (on the top) and residuals on S wave travel times (on the bottom).
velocity model were chosen to account for the present day knowledge of the crust in the area (Figure 12): a) two homogeneous velocity models (V P = 4.5 km/s and V P = 5.5 km/s).b) three velocity models characterized by three different constant velocity gradients.c) The velocity model used by the Nation Institute of Geophysics and Volcanology (INGV) to compile the instrumental catalogue of Italian earthquakes [Gruppo di lavoro CSTI 2001].

Figure 12 .
Figure 12.The nine starting Vp velocity models used in the linearized inversions.

Figure 13 .
Figure 13.(a) The nine inverted 1D velocity models; (b) the minimum (blue line) 1D velocity model.Red lines represent the error bounds.

Figure 14 .Figure 15 .
Figure 14.Position of the earthquake foci in the minimum 1D velocity model.

Figure 16 .
Figure 16.Plot of residuals between observed and theoretical travel times.Orange points refer to the three-layer velocity model; blue points refer to the minimum 1D velocity model.

Table 2 .
Weights associated to errors on P and S phase reading.