Supplemental Material for COSEISMIC DISPLACEMENTS ON ISCHIA ISLAND , REAL-TIME GPS POSITIONING CONSTRAINTS ON EARTHQUAKE SOURCE LOCATION

We use GNSS data to simulate the early seismic source location of a Mw 3.9 event occurred in the island of Ischia, Italy. The study suggests that real-time GNSS data can support the seismic location system in the early stage of the emergency phase. We demonstrate that this very shallow earthquake, triggered significant displacements at a few stations in less than half an hour. Using exclusively GPS data, the first location of the hypocenter was possible with a latency of only 20 minutes. Early upgrades of the offset field in the first two hours confirm the source location confined within 1-2 km in the horizontal plane and less than 1 km depth.


INTRODUCTION
GPS data have proven the potential of early tsunami warning [Blewitt et al., 2007;Sobolev et al., 2007;Singh et al., 2012;Hoechner et al., 2013;Savastano et al., 2017].In literature, so far, the focus is mainly on the advantages of the GPS data in recovering large earthquake parameters, but recently there is a growing interest in developing automated algorithms for real-time determinations of seismological parameters and deformation analysis for smaller seismic events [Melgar et al., 2015;Crowell et al., 2016].The real-time monitoring of ground deformations is not a technically demanding effort because of the broad spreading of GPS data, but may be extremely useful for assessing seismic location estimates, or for cross-checking purposes with other geophysical data in the early emergency phase.In this study we investigate the capability for singleepoch GPS positioning to characterize small to medium earthquake sources.We also assess the latency for the rapid seismic source location using only geodetic data for a recent event occurred on Ischia island.
Ischia is a volcanic island located in the northern sector of the Gulf of Naples, Italy.Thermal activity characterizes diffusively the island, fumarolic and thermal fluids emissions are widespread [Chiodini et al., 2004], mostly concentrated along the faults boarding Mt.Epomeo, a resurgent dome in the center of the island.The seismic activity is rather low and generally with low magnitude, however the historical record shows that Ischia island was struck by several earthquakes, a few of them causing severe damages and fatalities.The most outstanding in 1883 in which more than 2,000 people were killed and causing the total destruction of the village of Casamicciola [Carlino et al., 2010].On August 21th, 2017, an earthquake (Md=4; Mw 3.9) occurred,damaging a limited portion of the island in the Casamicciola area, reaching a maximum intensity of 8 on the European Macroseismic Scale (EMS) and causing two fatalities [Azzaro et al., 2017].
The Neapolitan volcanic area is constantly monitored by the INGV -Osservatorio Vesuviano (INGV-OV) institution.It operates a geodetic monitoring system, which includes a continuous GPS network of 37 permanent stations [Tammaro et al. 2013;De Martino et al. 2014].At present, 6 of these GPS stations are located on Ischia Island (Figure 1).An additional station on the island, part of the Regione Campania GNSS network, is also considered in this work (ISCH in Figure 1).The GPS observations were used to compute the coseismic displacements caused by the earthquake.We demonstrate that the coseismic offsets could be available in less than half an hour after the earthquake and could be useful in constraining the hypocenter position even for such small events.The proposed method is especially effective in volcanic or geothermal areas or for those events that are poorly constrained by the observing seismic network (e.g. on islands).

GNSS DATA ANALYSIS
The GPS data are sampled at 30 seconds and transmitted to the processing center at INGV-OV at Naples every 24 hours.Station positions are routinely computed and updated at daily rates using the Bernese v.5.0 software [Beutler et al., 2007].The routine precise processing demands standard products delivered by the International GNSS Service (IGS) with latency of about 2 weeks, the most relevant being, precise orbits, satellite clock corrections and Earth rotation parameters.Less The purple square box displays the source location obtained from the inversion of the GNSS displacement.
precise and predicted products are also released four times per day for real time and near real time use, these are the so called rapid and ultra-rapid IGS products.In particular the ultra-rapid predicted products are available in real time and allow adequate position estimates, virtually at every observation epoch, relative to a control point located hundreds kilometers away.In this work we estimate GPS time series using both ultra-rapid and precise products assessing the real-time capability to detect coseismic displacements.

PRECISE POSITION SOLUTION
Two different precise station position solutions are available from the INGV analysis centers, hereafter labeled OV and RM, both processing phase double differences.These two solutions are used to estimate a combined offset field for the seismic event occurred on the 21 st of August, and used here as the reference offset field for kinematic solution assessing.
The standard processing scheme for dual frequency GPS receivers uses the Ionosphere-free linear combination (L3) of phase observables while the phase ambiguities are solved using the QIF (Quasi Ionosphere Free) strategy [Mervart, 1995].
The troposphere delay is modeled using the dry-Niell model for the hydrostatic troposphere and the wet-Niell mapping function is used for the estimation of a hourly zenith delay correction [Niell, 1996].To model tropospheric azimuthal asymmetries, an additional daily horizontal gradient at each station is also estimated.Tidal ocean loading is modelled using the FES2004 coefficients [Lyard et al., 2006] provided by the Ocean Tide Loading web service (http://holt.oso.chalmers.se/loading/).At each station, the IGS08 absolute antenna phase corrections are applied.GPS satellite orbits and Earth Rotation Parameters are fixed to the combined official IGS precise products.
The reference frame realization differs substantially for the two solutions.The RM solution is obtained as loosely constraint estimates and subsequently projected on the ITRF2008 reference frame by a 4-parameter (3 translations and a scale factor) Helmert transformation.The OV solution instead, is obtained imposing three nonet translation conditions on a set of ITRF2008 reference stations (Minimum Constraint Solution).
The two precise positioning solutions are used to assess the relative quality of the individual solutions and in particular, to build a merged (in a least squares sense) solution for the coseismic offsets.The combined "consensus" solution is a rather common procedure for the geodetic product generation (see e.g.www.iers.org).
A rationale for this approach is that the consensus solution minimizes the chance of systematic errors and permits an assessment of the repeatability of the estimates.

KINEMATIC POSITION SOLUTION
We use the same processing software (Bernese v. 5.0) to build the kinematic solution.We select 7 stations at Ischia Island as rover stations using 3 additional stations as reference sites, located far from the rover stations in the Neapolitan area.To simulate a real-time application we use the IGS ultra-rapid products, processing data at different sampling rates (30s, 180s and 300s).The processing scheme may be summarized by the following steps: cycle-slip screening and outlier removal using the ionosphere-free linear combination (L3); ambiguity resolution using a sigma dependent strategy [Beutler, 2007] and finally the epoch wise position estimation based on L3 single differences whereas the reference sites are fixed to ITRF2008 coordinates.The dry part of the troposphere is modeled using the dry-Niell a-priori model and estimating two troposphere zenith delay parameters per day, using the wet-Niell mapping function.

PRECISE COSEISMIC OFFSET COMBINATION AND SOURCE LOCATION ESTIMATION
Two coseismic displacement fields (OV and RM) are estimated using the daily precise time series of the GPS stations.The coseismic offsets are computed as differences between the average positions, respectively 15 days before and 4 days after the event.Thus, the two displacement fields are combined in a least squares sense [Devoti, 2012], obtaining a "consensus" solution that minimizes possible sources of systematic errors.Figure 1 shows the displacements observed at the GPS stations for the horizontal and vertical spatial components.Among all stations, only two of them (OSCM and MEPO, Casamicciola Observatory and Mt.Epomeo, respectively) show significant horizontal offsets, characterized by a slightly divergent northward displacement, approaching each other in the S-N direction, whereas the vertical offsets are generally negligible.The average errors associated to the offsets in the vertical and horizontal components are respectively 8 mm and 2.5 mm. Figure 2 displays the residuals of the two input solutions with respect to the combined solution, respectively labeled RM and OV.It turns out that the residuals are well below their uncertainties and are unbiased, thus testifying an adequate repeatability and reproducibility of the displacement field.Acknowledg-ing this offset field and assuming a normal fault source in an elastic half-space with area proportional to the given magnitude, the estimated source location is provided in Figure 1.Unfortunately this conclusions can be drawn a few days after the seismic event, depending on how many daily positions are averaged out and when the rapid products are available.We will now push the real-time potentiality of the method, exploring the feasibility of estimating the offsets as soon as they become significant, given that real time data connection is available.

REAL TIME DISPLACEMENT DETECTION
We have to face the problem of detecting a transient signal in a growing time series, when an instantaneous change of the position coordinates (i.e.offsets) occurs.In this validation case we exploit the advantage of knowing the epoch at which the offset may occur, i.e. we recognize a trigger event after which we should be able to detect positions changes as soon as they will become reliable.Due to the noise process affecting the position estimates, the displacement can be measured only when an adequate number of observations are available.To quantify this time span, we implement a Student's t-test to check the null hypothesis that the population mean of the offset estimated after the trigger event, is equal to zero.In this context, the offsets are computed at each subsequent epoch providing the sample of offsets, that will be tested against the hypothesis of zero offset.The test will reject the null hypothesis as soon as the offset sample becomes significantly distinct from zero.The amount of time needed to detect an offset (i.e. the time at which the null hypothesis will be rejected) is the latency for the GPS observations to provide a significant displacement at the chosen significance level.Obviously the latency depends strongly on the signal-to-noise ratio (i.e. the displacement magnitude over the noise), the larger the signal (e.g. higher magnitude of the seismic event or shorter distance from the epicenter), the shorter will be the latency.Figures 3 and 4 show the updated offsets for MEPO and OSCM stations respectively, the two stations showing the most stable offsets.Similar offset plots for the other stations are reported in supplementary figures: a key point is in fact not only the matching of the significant offsets in the off-line and real-time approach, but also the lacking of false positive values in the real-time approach.The solid red line shows the offset estimated every 180 seconds after the trigger event in the three spatial components (up, east and north) and the inset on top of each graph shows the ratio (t/t 1% ) of the test statistics over the critical value (1% significance level), so that the region of rejection is t/t 1% >1.We recall that rejection of the null hypothesis implies the detection of a significant offset, thus when t/t 1% passes the unit threshold the displacement is significantly different from zero.Due to their impact on the latency of the estimations, we investigate three different time sampling intervals for the position estimates, namely the 30 sec (the GPS receiver sampling rate), 180 sec and 300 sec.The three time sampling do not differ much in the low frequency band (supplementary Figure S1), evidenced by the nearly equal efficiency to detect the same offset, but differ in high frequency noise which implies different latency time for the offset detection.
The MEPO station shows a clear northward offset of 1.0±1.2cm, recognized after 47 minutes, whereas the vertical and east components do not show any significant displacement, albeit a persistent negative easting offset of a few mm (i.e.towards West) is noticeable and the vertical is trending towards negative values.The only other station that shows any significant displacements is OSCM located about 2km north of MEPO.All three components show a positive offset, of 5.7±2.4,0.5±0.5 and 1.4±0.7cm,respectively in the vertical, east and north directions.The vertical uplift, initially of about 6 cm (significant after 17 minutes) reduces gradually in the following hours to zero (0.3±2.3cm) and the test statistics, after first overtaking the unit threshold, rebounds back to the region of acceptance (no detectable offset).The east and the north offsets, instead, show more stable estimates detected at about 30 and 20 minutes respectively after the trigger event.The only other station of the Ischia network that shows a non zero offset is FORI, where the east component is dis- '8.%/9:#;<=9#>=<?@∆N=1.0±1.1 (47min) placed by a small amount (0.2±0.6cm) that becomes significant about 3 hours after the trigger event (see supplementary Figure S2).

ASSESSMENT OF THE REAL TIME GNSS SEISMIC SOURCE FOR THE TEST CASE
The coseismic displacement obtained from the daily GPS time series is noticeable only in two stations of the Ischia GNSS network, at the Mt.Epomeo station (MEPO) and at the Casamicciola Observatory station (OSCM).The offsets are directed towards N (15.6 mm) and NNE (10.3 mm), respectively, whereas the vertical displacements are negligible apart from a slight subsidence observed at MEPO station.The observed displacements are compatible with the coseismic subsidence measured by the InSAR technique that shows a localized (~1 km 2 ) subsidence up to 4 cm, SW of Casamicciola (Gruppo di Lavoro INGV sul terremoto dell'isola di Ischia, 2017, De Novellis et al., 2018].Both results can be achieved only after a time delay that depends on data collection issues and products availabil-ity for the respective techniques. To prove the usefulness of the GPS data at the early stage of the earthquake detection, we perform a fast computation of the hypocenter by inverting the significant GPS displacements assuming a rectangular fault source [Okada, 1985].In the Apennines and Tyrrhenian area normal faulting is quite a common tectonic style [e.g.Malinverno and Ryan, 1986] but in volcanic domains also a point-like Mogi inflating/deflating source may be appropriate [Mogi, 1958].To study the more general problem we use a simple square fault plane and try to minimize the number of source parameters to be estimated.We fix the source dimensions according to the given magnitude using the relationship of the seismic moment M 0 =μAd [Aki, 1966], where μ is the rigidity modulus, A is the area of the fault surface and d is the slip dislocation.The surface area A is chosen according to the empirical relationship of Wells and Coppersmith, 1994 and fixed to 1 km 2 .The rake is fixed to normal faulting and the location, depth, dip and strike are chosen as the unknown parameters.The source parameters are solved using a nonlinear minimum search function fmincon [Matlab, 2012a].Figure  5 shows the projections of the source surfaces estimated at increasing time spans after the earthquake (20, 40, 60 and 120 minutes) using the pertinent offset field estimated at that epochs.The weighted root-mean-squared (wrms) of the misfit residuals, in all three spatial components, are on the order of 3 mm.The source locations are mainly scattered along the N-S direction, being roughly confined in a 1x3 km 2 area.The early source depth never exceeds 0.5 km, a result which is constrained by the few stations that show a non-zero offset.The estimated dip is very low, indicating a nearly horizontal source plane.However the geometry of the source (strike and dip) is weakly constrained, since the actual displacement field is quite under-sampled by the sparse GPS stations.

CONCLUSIONS
The estimated GNSS positions using 24h of satellite observations are highly accurate, stable and replicable under different processing approaches and currently used for calibration purposes in conjunction with InSAR data; however the potentiality of continuous and real-time nature of GNSS data is not commonly fully exploited.To assess the potentiality of a real-time positioning system, we set up a simulated kinematic solution and a procedure suited for offset detection that checks for a significant offset at any epoch step.The algorithm, based on classical statistical hypothesis testing, demonstrates the possibility of revealing displacements on the order of 1 cm in less than 20 minutes (50 and 90 minutes for different sampling rates) after the earthquake.The results concerning the intermediate sampling rate (180s) for the position time series are shown in this paper, enlightening also that different sampling rate does not affect substantially the estimates (e.g.see supplementary Figure S1), but affects the latency at which the estimate becomes reliable.This is an expected behavior since the noise does not decrease much as the sampling lengthens, providing less data points per unit of time and therefore allowing a longer latency.The real-time displacements are very close to the precisely estimated offsets but they are not identical.The differences may arise mainly for two reasons, first because of the high noise affecting the kinematic solution, possibly biased by orbital and tropospheric modeling errors; secondary because the kinematic solution provides the instantaneous timevarying displacement, whereas the static solution describes the long lasting deformations.
In conclusion, GPS kinematic positioning if linked with real-time data transmission, is suitable for detecting sudden displacements in the time series down to 1 cm or even below with sufficient reliability.We demonstrate that even a small magnitude earthquake, because very shallow, can trigger detectable GPS displacements useful for characterizing a preliminary source.In this case, the location of the hypocenter using exclusively GPS data, is possible as early as 20 minutes after the earthquake and is persistent after that time, thus being a valuable piece of information to be cross-validated with classical seismological inferences.This result testifies the capability of GPS data to provide quick evidence of earthquake source location, which can be crucial in assisting decision makers during the emergency management.
In the Italian area the hypocentral depth is on average deeper (8-12 km), so that higher magnitudes (>M5) will cause the same amount of displacements and the proposed method could be applied as well to monitor the low to moderate seismicity in other areas of the Italian peninsula.We evaluate that earthquake frequency in Italy reaches a few events per year for magnitude higher than M5 (and about 1 event per year for Mw > 5.4) [Gasperini et al., 2013].Thus it could be worthwhile to pursue real-time GNSS monitoring at national scale, or in special areas (e.g.Alpine and Apennines orogens) where major seismogenic structures have been identified [DISS Working Group, 2015) with the aim to possibly detect geodetic-derived seismic sources in the very first few minutes of the emergency response phase.

FIGURE 1 .
FIGURE 1. Map of the horizontal (blue arrows) and vertical (red arrows) coseismic displacement field of the Ischia GNSS network.The purple square box displays the source location obtained from the inversion of the GNSS displacement.

FIGURE 2 .
FIGURE 2. Residuals of the two displacement fields (OV and RM) with respect to the combined field at each GNSS station.

FIGURE 3 .
FIGURE 3.Real-time progression of the offset detection algorithm at MEPO station (Mt.Epomeo).The thinner blue line shows the significance ratio (t/t 1% ) whereas the heavy red line shows the offset estimated at each epoch in millimeters.The grey bar indicates the epoch at which the offsets become significant (t/t 1% >1), next to it, the value of the offset is given along with its 1-sigma error.

FIGURE 4 .
FIGURE 4. Real-time progression of the offset detection algorithm at OSCM station (Casamicciola).The thinner blue line shows the significance ratio (t/t 1% ) whereas the heavy red line shows the offset estimated at each epoch in millimeters.The grey bars indicate the epochs at which the offsets become significant (t/t 1% >1), next to it, the value of the offset is given along with its 1-sigma error.

FIGURE 5 .
FIGURE 5.Estimated source locations at subsequent epochs (in black) after the earthquake, as obtained from the real-time detected offsets.The purple square box displays the estimated source location shown in Figure 1.

FIGURE S2 .
FIGURE S2.Real-time offset detection algorithm worked out at five remaining stations of Ischia island.The thinner blue line shows the significance ratio (t/t 1% ) whereas the heavy red line shows the offset estimated at each epoch in millimeters.The grey bars indicate the epochs at which the offsets become significant (t/t 1% >1), next to it, the value of the offset is given along with its 1-sigma error.The sampling rate of the kinematic solution is 180 sec.