TY - JOUR
AU1 - Chen,, Kejie
AU2 - Liu,, Zhen
AU3 - Liang,, Cunren
AU4 - Song, Y, Tony
AB - SUMMARY Dense strong motion and high-rate Global Navigation Satellite Systems (GNSS) networks have been deployed in central Italy for rapid seismic source determination and corresponding hazard mitigation. Different from previous studies for the consistency between two kinds of sensor at collocated stations, here we focus on the combination of high-rate GNSS displacement waveforms with collocated seismic strong motion accelerators, and investigate its application to image rupture history. Taking the 2016 August 24 Mw 6.1 Central Italy earthquake as a case study, we first generate more accurate and longer period seismogeodetic displacement waveforms by a Kalman filter, then model the rupture behaviour through a joint inversion including seismogeodetic waveforms and InSAR observations. Our results reveal that strong motion data alone can overestimate the magnitude and mismatch the GNSS observations, while 1 Hz sampling rate GNSS is insufficient and the displacement is too noisy to depict rupture process. By contrast, seismogeodetic data enhances temporal resolution and maintains the static offsets that provide vital constraint to the reliable estimation of earthquake magnitude. The obtained model is close to the jointly inverted one. Our work demonstrates the unique usefulness of seismogeodesy for fast seismic hazard response. Space geodetic surveys, Earthquake ground motions, Earthquake hazards, Earthquake source observations 1 INTRODUCTION The central Apennine is a complex tectonic region where the continental collision between Nubia and Eurasia plates leads to the Alpine mountain belt. The region is also affected by the opening of Tyrrhenian basin to the west and the subduction of the Apennines from east to west and the underthrust of Adria microplate beneath Eurasia. Being one of the most seismically active areas in Italy, the region has experienced several medium-sized earthquake sequences over the past decades (Fig. 1) including the recent ones in 1997, 2009 and 2016. Figure 1. Open in new tabDownload slide Recent seismic sequences (Mw > 5.0) with focal mechanisms and distribution of GNSS sites (green triangles), strong motion stations (blue triangles). The inset map is adopted from Lavecchia et al. (2003) showing the regional tectonic setting of Tyrrhenian-Apennine system. Moment tensors are provided by Global Centroid Moment Tensor Catalog. GNSS and strong motion sites information can be found in Avallone et al. (2016) and Michelini et al. (2016). Figure 1. Open in new tabDownload slide Recent seismic sequences (Mw > 5.0) with focal mechanisms and distribution of GNSS sites (green triangles), strong motion stations (blue triangles). The inset map is adopted from Lavecchia et al. (2003) showing the regional tectonic setting of Tyrrhenian-Apennine system. Moment tensors are provided by Global Centroid Moment Tensor Catalog. GNSS and strong motion sites information can be found in Avallone et al. (2016) and Michelini et al. (2016). In response, strong motion (Michelini et al.2016) and continuous Global Navigation Satellite Systems (GNSS) networks have been established over this area (Fig. 1) during the past decades by several local agencies (Avallone et al.2016). The strong motion sensors are mainly designed to capture high sampling rate (100–200 Hz), real-time (or near real-time) coseismic dynamic motions. While originally GNSS stations aim to monitor long term crustal deformation, they have been demonstrated to be able to record coseismic displacement waveform directly as well, facilitating the new concept of ‘GPS seismology’ (e.g. Larson et al.2003). For example, the dynamic ground displacements caused by the 2009 Mw 6.3 L’Aquila and 2016 August 2016 Mw 6.1 events were recorded by a number of high-rate (1–10 Hz) continuous GNSS stations (Avallone et al.2011, 2016). GNSS derived static offsets and displacement waveforms are especially valuable in resolving the rupture plane and determining the slip distribution without bias which can be caused by different rupture velocities (Bock et al.2011). The high-rate GNSS displacement waveforms typically have an accuracy of <1 and 3 cm in the horizontal and vertical components over tens of minutes (Li et al.2013; Chen et al.2015, 2016), millimetre precision can also be expected when few minutes or less displacements are quantified (Geng et al.2017). Thus the deformation caused by medium sized earthquakes like the ones in central Italy is at the detectable limit of high-rate GNSS data. What's more, taking into account the fact that earthquake radiated energy decays very fast over epicentral distance, the coseismic displacements with high signal-to-noise ratio (SNR) should be restricted to small epicentral distance. As showed by Smalley (2009), for earthquakes around Mw 6, 1 Hz GNSS recordings of coseismic displacements will be aliased at very small epicentral distances and a higher sampling rates (e.g. 5 Hz) is needed to provide alias-free recordings. However, currently, the number of high-rate GNSS stations collecting data at higher than 5 Hz sampling rate still remains limited and 1 Hz sampling rate is dominant. As a result, few previous studies attempted to adopt GNSS displacement waveforms for seismic source analysis in this area. For instance, while Avallone et al. (2011) first assessed the high-rate (1–10 Hz) GNSS’s capability to capture P waves for the 2009 Mw 6.3 L’Aquila earthquake, they alternatively chose to use 100 Hz strong motion integrated and 10 Hz GNSS-derived velocity waveforms to retrieve the kinematic rupture history. It is understandable since GNSS velocity waveforms are more accurate than displacement waveforms because related errors (e.g. tropospheric delays) are reduced by epoch-differencing. Nonetheless, this approach is kind of ‘trade-off’ as it does not make direct use of GNSS displacements, and as a result, the fit to the displacement waveforms is often insufficient. As can be seen from Fig. 1, to ensure a good coverage, the two kinds of sensors are mostly distributed separately while a few of them are collocated (distance less than 1.5 km). Previous studies (e.g. Avallone et al.2011, 2016) noted this and emphasized the observation consistency between the collocated strong motion and GNSS sites. Recently, it has been shown that combining the two measurements, often termed as ‘seismogeodetic waveforms’, can generate longer period waveforms with better accuracy and higher sampling rate (e.g. Bock et al.2011; Melgar et al.2013; Tu et al.2013; Geng et al.2013a,b). The combined seismogeodetic waveforms enable the exploitation of the advantages of both sensors. The 2016 August 24 Mw 6.1 event was recorded by a number of well distributed collocated GNSS/strong motion stations (Fig. 2). In this study we focus on assessing the performance and application of seismogeodetic data in earthquake source study in central Italy. Figure 2. Open in new tabDownload slide Distribution of the collocated GNSS and strong motion stations (red squares) and static offsets estimated from high-rate displacement waveforms (blue vectors) and continuous three days’ daily observations (red vectors), both with 68 per cent confidence errors. The beach ball shows the 2016 August 24 main-shock location and black dots are aftershocks in the following three days obtained from eida.rm.ingv.it. The black rectangle outlines the surface projection of the fault plane adopted. Figure 2. Open in new tabDownload slide Distribution of the collocated GNSS and strong motion stations (red squares) and static offsets estimated from high-rate displacement waveforms (blue vectors) and continuous three days’ daily observations (red vectors), both with 68 per cent confidence errors. The beach ball shows the 2016 August 24 main-shock location and black dots are aftershocks in the following three days obtained from eida.rm.ingv.it. The black rectangle outlines the surface projection of the fault plane adopted. We reprocessed the high-rate GNSS data recorded on 23rd and 24th August, and utilized sidereal filtering to reduce the multipath effect. We employed a Kalman filter to combine GNSS displacement waveforms and seismic strong motion data at collocated stations to generate seismogeodetic data. We examined the spatiotemporal rupture details of this earthquake by using high-rate GNSS, strong motion and seismogeodetic data, respectively. To address the advantages of seismogeodetic data for seismic source estimation, we also included the deformation observations from L-band ALOS-2 SAR sensor, which covers the entire fault area, to perform joint inversions with seismogeodetic data and set the preferred model as a benchmark. We note that the InSAR observations were recorded three days after the main shock during which a number of aftershocks occurred. To investigate whether the large aftershocks can cause significant surface deformation and contaminate InSAR observation, we also modelled the largest Mw 5.5 aftershock using strong motion records. Based on these results, we demonstrate the contribution of seismogeodesy and discuss the need for more collocated GNSS and strong motion sites in the study area for fast and reliable earthquake source imaging. 2 DATA AND METHODS 2.1 Integration of high-rate GNSS and strong motion The raw high-rate GNSS observations have been archived by Istituto Nazionale di Geofisica e Vulcanologia (INGV) and are publically available (ftp://gpsfree.gm.ingv.it/amatrice2016/hrgps/data/). Precise point positioning (PPP) was employed to obtain centimetre-level epoch-wise position estimate. We adopted earth rotation parameters, final satellite orbits and 5 s interval precise satellite clocks from Center for Orbit Determination in Europe. Errors such as antenna phase centre variations, phase wind-up, relativistic effect, dry tropospheric delay, ocean tide loading, pole and solid earth tides were corrected according to existing models (Gérard & Luzum 2010). We assumed zenith tropospheric delay as a random walk parameter with a power spectral density of 0.002 mm s−1/2 while the receiver clocks are estimated epoch by epoch as white noise. No constraints between neighbouring epochs were imposed and the positions estimated were aligned to International Terrestrial Reference Frame 2008 (ITRF08). Note that for some stations the observations were interrupted shortly after the main shock. To ensure a high-level confidence of the phase ambiguity resolution, we processed the continuous GNSS data from the preceding day. The PANDA software package (Liu & Ge 2003) was utilized for the GNSS data processing. Figs 3(a) and (b) show an example of PPP position displacement waveforms at station NRCI on 23rd and 24th August. The position time-series display similar low frequency noise (Figs 3a and b) suggesting that they could be caused by local multipath effect. To minimize the effect, we shifted the time-series of 23rd August by 246 s (sidereal period is set to be 23 h 56 m 4 s) and then subtracted it from time-series of the 24th based on the modified sidereal filtering (Choi et al.2004; Geng et al.2017). Figs 3(c) and (d) demonstrate the improvement from this post-processing. Specifically, in Fig. 3(c), while the coseismic signal is clearly seen in all three components, long-period drifts are found in north-south, up-down directions. These are mostly removed in Fig. 3(d). Avallone et al. (2016) processed the high-rate GNSS data using both PPP and double difference approaches (ftp://gpsfree.gm.ingv.it/amatrice2016/hrgps/). But they didn’t mitigate the multipath effect. As a result, the low frequency noise was still present in their time-series, especially in the vertical components. Figure 3. Open in new tabDownload slide Epoch-wise PPP results at station NRCI. (a) 500 s long PPP results on 23rd August; (b) 500 s long PPP results on 24th August, including the co-seismic displacement waveforms; (c) Zoomed view of the co-seismic displacement waveforms shown in (b) over the time of 5800–5880 s.; (d) Filtered co-seismic displacement waveforms after sidereal filtering. Figure 3. Open in new tabDownload slide Epoch-wise PPP results at station NRCI. (a) 500 s long PPP results on 23rd August; (b) 500 s long PPP results on 24th August, including the co-seismic displacement waveforms; (c) Zoomed view of the co-seismic displacement waveforms shown in (b) over the time of 5800–5880 s.; (d) Filtered co-seismic displacement waveforms after sidereal filtering. We corrected the zero baseline offsets from 200 Hz strong motion data and integrate high-rate GNSS displacements with strong motion records to generate seismogeodetic velocity and displacement waveforms at eight collocated sites using a multirate Kalman filter (Bock et al.2011). In the Kalman filtering the variances for GNSS and strong motion data were estimated from previous 60 s GNSS position time-series and 6 s strong motion records, respectively. In Fig. 4, we show the GNSS and seismogeodetic displacement waveforms at station NRCI/NRC. The two displacement waveforms are in good agreement. However, we can still find some differences at some epochs after the arrival of shaking. For example, in the north-south component, at 01:36:43 UTC, GNSS recorded instantaneous displacement is about 5 cm larger than the combined seismogeodetic records. Please note that the horizontal precision of GNSS is better than 1 cm and we should not treat this bias as outlier. This can be explained by the fact that the GNSS station locates at the southeast of the strong motion station and is closer to the epicentre. Coincidentally, the slip is featured by N–NW and S–SE directivity (Tinti et al.2016). Accordingly, we set the strong motion site to be the location of seismogeodetic waveforms in inversion. Figure 4. Open in new tabDownload slide 1 Hz GNSS displacement waveforms (grey dots) at GNSS station NRCI and 200 Hz smoothed seismogeodetic displacement waveforms (blue, green and red lines) at station NRCI/NRC. Figure 4. Open in new tabDownload slide 1 Hz GNSS displacement waveforms (grey dots) at GNSS station NRCI and 200 Hz smoothed seismogeodetic displacement waveforms (blue, green and red lines) at station NRCI/NRC. In addition, three components of coseismic offsets (see in Fig. 2) at each station were extracted from displacement waveforms following the method proposed by Liu et al. (2014), and both displacement and velocity waveforms were cut to be 35 s long with respect to the origin time for inversion. One of the seven stations, INFN, did not have GNSS data on 23rd August, which prevent applying sidereal filtering. As such, this station was given less weight in the inversion. For slip inversion, all of the derived 200 Hz seismogeodetic waveform displacements were low pass filtered with a cut frequency of 0.5 Hz and the sampling rate decimated to 5 Hz. Specially, as noted by Tinti et al. (2016), significant topography exists in this area, to account for the effect of topography, we followed their approach to set the average elevation of the seven stations as 750 m. 2.2 InSAR observations We processed two pairs of InSAR data acquired by JAXA’s L-band ALOS-2 satellite using newly developed software added to the JPL ISCE InSAR processing software (Liang & Fielding 2017a, Liang & Fielding 2017a,b). One pair was acquired on 2015 September 9 and 2016 August 24 on ascending track 197. The other pair was acquired on 2016 May 25 and 2016 August 31 on descending track 92. All the data were acquired in stripmap mode. Two frames of each track were processed. The frames were mosaicked at the single look interferogram level. The 1 arcsec (∼30 m) SRTM DEM (Farr et al.2007) was used to correct the topographic phase. The interferograms were unwrapped by SNAPHU (Chen & Zebker 2002). To reduce computational cost, we downsampled the original InSAR interferogram using a resolution based approach (e.g. Lohman & Simons 2005) and kept all observations within longitude 13.1° E-13.4° E and latitude 42.6° N-42.9° N. As a result, 1465 InSAR samples were included for the joint inversion. 2.3 Fault parameterization and Green's functions calculations We adopted fault geometry with the strike and dip angles (155° and 49°) provided by INGV (http://cnt.rm.ingv.it/event/7073641) based on the moment tensor computation. The hypocentre location is (13.23°, 42.70°) with a depth of 6.8 km (6.05 km below the sea level). The finite fault plane was then discretized into 20 and 10 subfaults along strike and dip, with 2.0 km spacing. We employed a non-negative least-squares inversion in which the rake angle of the slip vector was allowed to vary between 240° and 300° to ensure a normal fault mechanism. In addition, Laplacian regularization was adopted to constrain the second-order gradient of each parameter to be zero and ensure inversion stability. We used a multitime window inversion, in which the source time function of each subfault was parameterized with 10 symmetric triangles with 0.2 s rise time and 0.2 s shift. The frequency-wavenumber integration method (Zhu & Rivera 2002) was adopted to compute Green's functions for static GNSS and InSAR offsets, and seismogeodetic waveforms using a 1-D layered velocity structure model by Herrmann et al. (2011). The rupture speed was set to be 2.8 km/s (the same as adopted by Tinti et al.2016), which is about 0.9 times the maximum shear wave speed of the velocity layer spanned by the fault model. 2.4 Relative weighting between data sets Relative weighting has always been a tricky issue for joint inversion. Up to now, there has not been a well-established rule of thumb to determine an optimized weighting factor. In this study, we applied a relative weighting based on respective variances of each data set when fitting all data sets: coseismic static offsets estimated from seismogeodetic displacement waveforms were assigned to have the largest weight 1, while InSAR observations and seismogeodetic waveforms weights were 0.3. For the station INFN, its static offsets and waveforms had the weights of 0.3 and 0.01, respectively. 2.5 Checkerboard test We conducted a series of checkerboard tests to investigate data resolution and inversion stability. For the input model, the rake angle was 270°, the rupture velocity was 2.8 km/s and the rise time was 1 s. Each patch contains 2 × 2 subfaults (2 km × 2 km) and was prescribed with 1 m of slip for every other patch. The synthetic displacement and velocity waveforms were generated at seven seismogeodetic stations. In the tests, we prescribed all the parameters to be the same as in fault parameterization in section 2.2. Besides, spatial and temporal smoothing was applied to regularize the results, similar to the smoothing factors that we applied to the real data. Results are demonstrated in Fig. 6. It can be seen that the input slip pattern in general is well recovered, especially the shallow portion. 3 RESULTS AND DISCUSSIONS In the following, all kinematic source inversions were conducted by the open source code MudPy (Melgar & Bock 2015), and the preferred model was chosen based on Akaike's Bayesian Information Criterion formalism (Fukahata et al.2003). 3.1 Negligible deformation caused by Mw 5.5 aftershock As mentioned before, in this study, we set the joint inversion result as the benchmark to evaluate the performance of the seismogeodetic data. Note that the InSAR coseismic deformation measurements contains both main shock and aftershock contributions due to the latency of post-earthquake SAR acquisitions (see Fig. 5 and sentinel-1A provided by European Space Agency http://www.esa.int/spaceinimages/Images/2016/08/Ground_displacement_from_Italy_s_earthquake). To avoid potential inconsistency between the displacement measurements of seismogeodetic and InSAR data, first we investigated the impact of Mw 5.5 aftershock using only strong motion waveforms archived by RAN (http://ran.protezionecivile.it/IT/index.php). To reduce the scattering effects due to long propagation distance, we only used stations within 30 km from the epicentre in the inversion. The raw accelerator observations were integrated into velocity waveforms after demeaning, detrending and tapering. The lower frequency response of strong motion sensors is not reliable because of baseline shift. As a result, a bandpass filter with corner frequencies of 0.02 and 0.5 Hz was applied to the obtained velocity waveforms and the sampling rate was then decimated to 5 Hz from original 200 Hz, and the waveforms were cut with a 30 s long time window with respect to the origin time. We adopted the aftershock hypocentre location (longitude 13.15°, latitude 42.79°, depth 8 km), origin time (UTC, 24–08-2016 02:33:28), strike and dip angles (327°, 43°) of fault plane from INGV (http://cnt.rm.ingv.it/en/event/7076161). The preferred model and corresponding data fits are illustrated in Figs 7 and 8, respectively. Figure 5. Open in new tabDownload slide Co-seismic line-of-sight displacement from L-band ALOS-2. The left is from ascending track 192 and the right is from descending track 92. Figure 5. Open in new tabDownload slide Co-seismic line-of-sight displacement from L-band ALOS-2. The left is from ascending track 192 and the right is from descending track 92. Figure 6. Open in new tabDownload slide Checkerboard tests. The left is input model and the right is inverted result. Red star represents the earthquake epicentre. Figure 6. Open in new tabDownload slide Checkerboard tests. The left is input model and the right is inverted result. Red star represents the earthquake epicentre. Figure 7. Open in new tabDownload slide Map-view of slip distribution from the preferred model of Mw 5.5 event. Black triangles represent strong motion stations adopted. Moment release as a function of time for the event is shown in the upper right inset. Figure 7. Open in new tabDownload slide Map-view of slip distribution from the preferred model of Mw 5.5 event. Black triangles represent strong motion stations adopted. Moment release as a function of time for the event is shown in the upper right inset. Figure 8. Open in new tabDownload slide Velocity waveform fits at 10 selected strong motion stations from the preferred model of Mw 5.5 event. The observed and synthetic data are in black and red, respectively. The numbers at right indicate the maximum amplitude values for each waveform. Figure 8. Open in new tabDownload slide Velocity waveform fits at 10 selected strong motion stations from the preferred model of Mw 5.5 event. The observed and synthetic data are in black and red, respectively. The numbers at right indicate the maximum amplitude values for each waveform. Our inverted moment magnitude is ∼5.5. Most of slips concentrate at a quite deep depth (∼9 km), which would cause less than 1.5 cm vertical and 4 mm horizontal deformation (see Supporting Information Fig. S1). Clearly, even as the largest aftershock, the Mw 5.5 event contributes only marginally to the overall deformation. In fact, this event was also modelled by Scognamiglio et al. (2016), and more strong motion satiations were used in their study. We find that our model and their model have general consistencies in the bilateral propagation and two deep slip asperities around the epicentre. In addition, we compared the static offsets estimated from seismogeodetic displacements (caused by main shock) with solutions reported by INGV (INGV Working Group et al.2016) which were obtained based on 3 d’s GNSS observations (including the effects of both main shock and aftershocks) after the event (see in Fig. 2), the estimates are nearly the same. Thus it is reasonable to combine seismogeodetic waveforms and InSAR observations in this study despite the different time span covered by each data set. Similarly, we believe the effects of other smaller aftershocks are negligible. 3.2 Rupture features of the main shock from joint inversion Fig. 9 (bottom right) shows our preferred model of the main shock from individual data set and joint inversion (Fig. 9d). The seismogeodetic displacement and velocity waveform data fits for the joint model are shown in Fig. 10. Fig. 11 shows the InSAR model fit and residuals. The InSAR residuals of both tracks are mostly correlated with the topography in this area, suggesting most residuals are likely related to the tropospheric noise. In particular, there are relatively larger InSAR residuals in descending tracks (>5 cm) to the southeast of the epicentre. These residuals are also correlated with the local topography. Variance reductions (VR) are calculated by the method described in Melgar et al. (2012), and VR for static GPS offsets, InSAR, displacement waveforms, velocity waveforms are 90 per cent, 61 per cent, 57 per cent, 73 per cent. Figure 9. Open in new tabDownload slide Map-view of slip distributions from (a) strong motion, (b) GNSS displacement waveforms, (c) seismogeodetic and (d) joint inversion. Observed and predicted horizontal static displacements used in the inversion are plotted with blue and red arrows, respectively. For subplot (a), the synthetic displacements are forwarded from the model. Moment release as a function of time for the event is shown in the upper right inset. Enlarged map-view of slip distribution is provided in the Supporting Information Fig. S1. Figure 9. Open in new tabDownload slide Map-view of slip distributions from (a) strong motion, (b) GNSS displacement waveforms, (c) seismogeodetic and (d) joint inversion. Observed and predicted horizontal static displacements used in the inversion are plotted with blue and red arrows, respectively. For subplot (a), the synthetic displacements are forwarded from the model. Moment release as a function of time for the event is shown in the upper right inset. Enlarged map-view of slip distribution is provided in the Supporting Information Fig. S1. Figure 10. Open in new tabDownload slide Seismogeodetic displacement (left) and velocity (right) waveform fits at seven collocated GNSS/strong motion stations for the joint model. The observed and synthetic data are in black and red, respectively. The numbers at right indicate the maximum amplitude values for each waveform. Figure 10. Open in new tabDownload slide Seismogeodetic displacement (left) and velocity (right) waveform fits at seven collocated GNSS/strong motion stations for the joint model. The observed and synthetic data are in black and red, respectively. The numbers at right indicate the maximum amplitude values for each waveform. Figure 11. Open in new tabDownload slide Map view of ALOS-2 line-of-sight displacement (after down sampling) (left), model prediction (middle) and residuals (right). The top panels are for ascending track and the bottom ones are for descending track. Figure 11. Open in new tabDownload slide Map view of ALOS-2 line-of-sight displacement (after down sampling) (left), model prediction (middle) and residuals (right). The top panels are for ascending track and the bottom ones are for descending track. As demonstrated in our joint model, the rupture propagated in an asymmetrical pattern with two main asperities, one NW to the hypocentre and another one above the hypocentre. We find most of the slip occurred at depths from 3 to 6 km, with a length of about 20 km along strike, which confirms the results of previous published model. The maximum slip was about 0.9 m, locating just above the hypocentre. The total seismic moment reached 1.8 × 1018 Nm, equivalent to a magnitude of Mw 6.1. The majority of the seismic energy was released within 8 s of the origin time. Fig. 12 depicts the rupture evolution with snapshots at 1 s interval. As can be seen, the rupture initiated near the hypocentre area and propagated bilaterally toward updip very fast in the first 2 s. After that, the main rupture front reached its updip limit, which may be related to shallow clay-rich sedimentary layer (Mattei et al.1997), and extended further along the strike towards southeast and northwest, with the northwest propagation dominating. The northwest rupture accelerated around 6 s. Both slip patches were at the shallow depth with the aftershocks located predominantly beneath coseismic slip (see in Fig. 13). The complementary aftershock distribution was possibly caused by coseismic stress change and early after slip. Figure 12. Open in new tabDownload slide Snapshots of the spatiotemporal history of the slip at 1 s interval, the blue star denote the location of hypocentre, and the grey circles are reference rupture fronts moving out at 2.8, 3.0, and 3.2 km s−1. Figure 12. Open in new tabDownload slide Snapshots of the spatiotemporal history of the slip at 1 s interval, the blue star denote the location of hypocentre, and the grey circles are reference rupture fronts moving out at 2.8, 3.0, and 3.2 km s−1. Figure 13. Open in new tabDownload slide Tile view of slip distribution from joint inversion. The green star denotes the hypocentre. Black circles are the double difference relocated aftershocks from INGV (http://iside.rm.ingv.it/iside/). Figure 13. Open in new tabDownload slide Tile view of slip distribution from joint inversion. The green star denotes the hypocentre. Black circles are the double difference relocated aftershocks from INGV (http://iside.rm.ingv.it/iside/). The joint model parameters are listed in Supporting Information Fig. S2. Furthermore, we compared our model (Fig. 9c) with the one obtained by Tinti et al. (2016) which was inferred from 26 strong motion stations. Overall, the two slip models agree with each other quite well: shallow slip patches with two asperities, bilateral propagation towards to NW and SE. The main difference is the rupture area, our model shows obvious slips extending as far as 15 km to the SE from the epicentre while theirs are limited to only 5 km. Possible cause is the sparse distribution and poor coverage of seismogeodetic stations as there are only seven stations were adopted in our study. Another possible factor is that we did not apply sidereal filtering to the southernmost station INFN due to lack of observations on 23rd August despite that it was assigned the least weight in the inversion. 3.3 Contributions of seismogeodetic data for reliable fast seismic source inversion To evaluate the potential benefits of seismogeodetic data for seismic source imaging, we also present slip inversions using only strong motion, high-rate GNSS waveforms (see in Figs 9a–c), respectively. Compared with our preferred model (Fig. 9d) in Section 3.2, we find that: With respect to strong motion only inversion, the data fits (see in Fig. 14) are much better (VR = 82 per cent). However, the synthetic static offsets from the favoured model slip model cannot fit the GNSS observations quite well (see Fig. 9a). Besides, the maximum slip is larger (around 1.0 m) and a third asperity appears in southeast of the fault plane. The derived magnitude is ∼Mw 6.25, which is also greater. Generally, strong motion data seems to overestimate the event size, which might be caused by the site effect (Bindi et al.2011). Figure 14. Open in new tabDownload slide Strong motion data fits. Black and red represent observations and synthetics, respectively. The strong motion accelerators have been first integrated to velocities and bandpass filtered between 0.02 and 0.5 Hz, then decimated from 200 to 5 Hz. The numbers above each waveform indicate the peak amplitude. Figure 14. Open in new tabDownload slide Strong motion data fits. Black and red represent observations and synthetics, respectively. The strong motion accelerators have been first integrated to velocities and bandpass filtered between 0.02 and 0.5 Hz, then decimated from 200 to 5 Hz. The numbers above each waveform indicate the peak amplitude. In contrast, the data fits (see in Fig. 15) for 1 Hz GNSS displacement waveforms inversion are not that good, especially the vertical components. Besides, only one major compact asperity is retrieved (see in Fig. 9b). The synthetic static offsets, however, agree with the observation nicely, and the magnitude is Mw 6.1. Possible reasons are relatively low sampling rate and SNR of the observations. In this case, 1 Hz sampling rate may not be sufficient to model dynamic rupture when all of the stations have a small epicentral distance (Smalley 2009). The peak ground displacements are only up to several centimetres, just around the detection level of the GNSS. As a result, the pre-event fluctuations can be incorrectly treated as seismic signals. Actually, using high-rate displacement waveforms to constrain rupture process of medium-sized events tends to fit the data not very well. For example, Melgar et al. (2015) use the high-rate GPS to model the 2014 Mw 6.1 Napa earthquake and show a moderate data fit. Figure 15. Open in new tabDownload slide GNSS displacement waveforms model fit. Black and red represent observations and synthetics, respectively. The 1 Hz GNSS displacement waveforms have been low-pass filtered by 0.49 Hz. The numbers above each waveform indicate the maximum amplitude. Figure 15. Open in new tabDownload slide GNSS displacement waveforms model fit. Black and red represent observations and synthetics, respectively. The 1 Hz GNSS displacement waveforms have been low-pass filtered by 0.49 Hz. The numbers above each waveform indicate the maximum amplitude. In comparison, the slip model inferred from seismogeodetic waveforms is close to the joint inversion model. The fits to seismogeodetic data (see in Fig. 16) are quite satisfactory, and the magnitude is also Mw 6.1. We find that there is a clear improvement for displacement waveform fits. As noted, integration of strong motion and high-rate GNSS reduces the noise level of displacement and velocity waveforms, enhances temporal resolution while maintaining the static offsets, which provide vital constrains on the reliable estimation of earthquake magnitude. These advantages as highlighted by kinematic source inversion confirm the unique contribution of seismogeodesy for rapid response and modelling of earthquake rupture. Figure 16. Open in new tabDownload slide Seismogeodetic displacement and velocity waveform fits from the slip model inferred from seismogeodetic data. Black and red represent observations and synthetics, respectively. Both displacement and velocity waveforms have been low-pass filtered by 0.5 Hz. The numbers above each waveform indicate the maximum amplitude. Figure 16. Open in new tabDownload slide Seismogeodetic displacement and velocity waveform fits from the slip model inferred from seismogeodetic data. Black and red represent observations and synthetics, respectively. Both displacement and velocity waveforms have been low-pass filtered by 0.5 Hz. The numbers above each waveform indicate the maximum amplitude. Meanwhile, it should be pointed out that the resolved coseismic slip distribution depends on many other factors such as the data sets adopted, choice a particular inversion scheme, variant of seismological parameterization, geometry of the array, accuracy of earth structure, spatial and temporal smoothing constrain. It has been suggested that one should be cautioned against dogmatic interpretation of inhomogeneous features on inverted slips, except gross characteristics (Beresnev 2003). For example, there are at least four published finite source models for this Mw 6.1 event (Tinti et al.2016; Chiaraluce et al.2017; Huang et al.2017; Liu et al.2017), each model is different but shares some common basic features. For fast finite source modelling, we aim to retrieve gross features of an earthquake. Seismogeodetic data synthesis the advantages of strong motion and GNSS observations, help to reduce slip uncertainties caused by each data set and contribute to refine the main features, which have been demonstrated in our manuscript (see Fig. 9). 3.4 Implications for the deployment of colocated GNSS/strong motion stations in Italy While real-time high-rate GNSS data has been proved to be a valuable complement for fast seismic source inversion and earthquake early warning (see e.g. Allen & Ziv 2011; Colombelli et al.2013), due to model errors such as troposphere delays, multi-path effects, satellite clock bias, particularly in real-time scenario, the accuracy of real-time GNSS is limited in its application to the medium-sized earthquakes (Melgar et al.2015). Italy possesses one of the densest earthquake monitoring network (real-time GNSS, strong motion, etc.) in the world. However, statistically, most of the destructive earthquakes in Italy are not larger than Mw 6.5, which limits the use of the GNSS to reliably resolve real-time coseismic offsets and displacement waveforms especially for the stations at larger epicentral distance. In fact, after reviewing the published models for Italy earthquakes, it is found that very few uses the high-rate GNSS observations. Our own test of rupture inversions using only GNSS displacement waveforms shows the data fit is poor and the results are not trustworthy because of relatively large noise level. Instead, as demonstrated in this study, integration of GNSS and strong motion observations produces more accurate seismogeodetic waveforms and improves the finite source inversion, which has the potential to be used for rapid response and earthquake early warning. The colocation of GNSS and strong motion stations in Italy will be more meaningful than other places given its proneness to magnitude 6 or so earthquakes. For other two events that occurred on 26th and 30th October 2016, the geometry of colocated stations is not as ideal as the 26th August one (see Fig. 17). This reduces the constraint from seismogeodetic records on the rupture details. Our results highlight the importance to have more GNSS and strong motion stations collocated in Italy in the future. Figure 17. Open in new tabDownload slide Distribution of collocated GNSS/strong motion stations for the other two major events of the 2016 central Italy seismic sequence: (a) 26th October Mw 5.9 earthquake; (b) 30th October Mw 6.5 earthquake. Black squares mean seismogeodetic stations and the beach ball locates the epicentre. Figure 17. Open in new tabDownload slide Distribution of collocated GNSS/strong motion stations for the other two major events of the 2016 central Italy seismic sequence: (a) 26th October Mw 5.9 earthquake; (b) 30th October Mw 6.5 earthquake. Black squares mean seismogeodetic stations and the beach ball locates the epicentre. 4 CONCLUSIONS In this study, we demonstrated the application of seismogeodesy for fast and reliable seismic source analysis through the example of the 2016 Mw 6.1 Central Italy earthquake. We reanalysed the GNSS data, and employed sidereal filter to reduce low frequency noise in the displacement time-series. Through Kalman filter based combination, we generated more accurate seismogeodetic waveforms and used them in the inversion. InSAR provided 2-D high resolution ground deformation data (while the most of the errors come from troposphere) for constraining and validating the seismogeodetic inversion. Using strong motion data, we find that the largest Mw 5.5 aftershock contributed marginally to the whole deformation, which ensures the feasibility of including InSAR observations for joint inversion. The joint inversion model indicates that the total moment for the main event was 1.8 × 1018 Nm (∼Mw 6.1) with a peak slip of nearly 90 cm locating southeast to the epicentre. The rupture spread bilaterally with two slip asperities being revealed, one NW and the other SE to the epicentre. Taking the joint inversion model as a reference, we find that the slip model inferred from individual strong motion data tends to overestimate the magnitude, and cannot fit the coseismic static offsets well. By contrast, the GNSS displacement waveforms are too noisy but provide some valuable constraints on magnitude estimation. Seismogeodetic data enhances temporal resolution and retains the long period information. The inferred slip model from seismogeodetic data is close to the jointly inverted one. Our results suggest more GNSS/strong motion stations should be colocated in Italy where most of the earthquakes are medium-sized to fully exploit the high-rate GNSS for earthquake study. ACKNOWLEDGEMENTS The research described here was conducted at the Jet Propulsion Laboratory (JPL), California Institute of Technology, under contracts with the National Aeronautics and Space Administration (NASA). The authors also thank Dr Han Yue from California Institute of Technology, Dr Shengji Wei from Nanyang Technological University and Dr Wanpeng Feng from Natural Resources Canada for beneficial discussions on seismic inversions. The inversion codes were originally developed by Dr Diego Melgar at University of California, Berkeley. JPL high-performance computing provided the computational resources. The work also made use of GMT software. We also thank two anonymous reviewers, whose comments and suggestions are very helpful for improving the manuscript. REFERENCES Allen R.M. , Ziv A. , 2011 . Application of real-time GPS to earthquake early warning , Geophys. Res. Lett. , 38 , doi:10.1029/2011GL047947. https://doi.org/10.1029/2011GL047947 WorldCat Avallone A. et al. , 2011 . Very high rate (10 Hz) GPS seismology for moderate-magnitude earthquakes: the case of the Mw 6.3 L’Aquila (central Italy) event . J. geophys. Res. , 116 , doi:10.1029/2010JB007834 . https://doi.org/10.1029/2010JB007834 WorldCat Avallone A. et al. , 2016 . Coseismic displacement waveforms for the 2016 August 24 Mw 6.0 Amatrice earthquake (central Italy) carried out from High-Rate GPS data , Ann. Geophys. , 59 , doi:10.4401/ag-7275 . WorldCat Beresnev I.A. , 2003 . Uncertainties in Finite-Fault Slip Inversions: to what extent to believe? (A critical review) , Bull. seism. Soc. Am. , 93 , 2445 – 2458 . https://doi.org/10.1785/0120020225 Google Scholar Crossref Search ADS WorldCat Bindi D. , Luzi L. , Parolai S. , Di Giacomo D. , Monachesi G. , 2011 . Site effects observed in alluvial basins: the case of Norcia (Central Italy) , Bull. Earthq. Eng. , 9 , 1941 – 1959 . https://doi.org/10.1007/s10518-011-9273-3 Google Scholar Crossref Search ADS WorldCat Bock Y. , Melgar D. , Crowell B.W. , 2011 . Real-time strong-motion broadband displacements from collocated GPS and accelerometers , Bull. seism. Soc. Am. , 101 , 2904 – 2925 . https://doi.org/10.1785/0120110007 Google Scholar Crossref Search ADS WorldCat Chen C.W. , Zebker H.A. , 2002 . Phase unwrapping for large SAR interferograms: statistical segmentation and generalized network models , IEEE Trans. Geosci. Remote Sens. , 40 , 1709 – 1719 . https://doi.org/10.1109/TGRS.2002.802453 Google Scholar Crossref Search ADS WorldCat Chen K. , Ge M. , Li X. , Babeyko A. , Ramatschi M. , Bradke M. , 2015 . Retrieving real-time precise co-seismic displacements with a standalone single-frequency GPS receiver , Adv. Space. Res. , 56 , 634 – 647 . https://doi.org/10.1016/j.asr.2015.04.029 Google Scholar Crossref Search ADS WorldCat Chen K. , Ge M. , Babeyko A. , Li X. , Diao F. , Tu R. , 2016 . Retrieving real-time co-seismic displacements using GPS/GLONASS: a preliminary report from the September 2015 Mw 8.3 Illapel earthquake in Chile , Geophys. J. Int. , 206 , 941 – 953 . https://doi.org/10.1093/gji/ggw190 Google Scholar Crossref Search ADS WorldCat Chiaraluce L. et al. , 2017 . The 2016 central Italy seismic sequence: a first look at the mainshocks, aftershocks, and source models . Seismol. Res. Lett. , 88 , 757 – 771 . https://doi.org/10.1785/0220160221 Google Scholar Crossref Search ADS WorldCat Choi K. , Bilich A. , Larson K.M. , Axelrad P. , 2004 . Modified sidereal filtering: Implications for high-rate GPS positioning , Geophys. Res. Lett. , 31 , 1 – 4 . https://doi.org/10.1029/2004GL021621 Google Scholar Crossref Search ADS WorldCat Colombelli S. , Allen R.M. , Zollo A. , 2013 . Application of real-time GPS to earthquake early warning in subduction and strike-slip environments , J. geophys. Res. , 118 , 3448 – 3461 . https://doi.org/10.1002/jgrb.50242 Google Scholar Crossref Search ADS WorldCat Farr T.G. et al. , 2007 . The shuttle radar topography mission , Rev. Geophys. , doi:10.1029/2005RG000183. https://doi.org/10.1029/2005RG000183 WorldCat Fukahata Y. , Yagi Y. , Matsu’ura M. , 2003 . Waveform inversion for seismic source processes using ABIC with two sorts of prior constraints: comparison between proper and improper formulations , Geophys. Res. Lett. , doi:10.1029/2002GL016293 . https://doi.org/10.1029/2002GL016293 WorldCat Geng J. , Bock Y. , Melgar D. , Crowell B.W. , Haase J.S. , 2013a . A new seismogeodetic approach applied to GPS and accelerometer observations of the 2012 Brawley seismic swarm: Implications for earthquake early warning , Geochem. Geophys. Geosyst. , 14 , 2124 – 2142 . https://doi.org/10.1002/ggge.20144 Google Scholar Crossref Search ADS WorldCat Geng J. , Melgar D. , Bock Y. , Pantoli E. , Restrepo J. , 2013b . Recovering coseismic point ground tilts from collocated high-rate GPS and accelerometers , Geophys. Res. Lett. , 40 , 5095 – 5100 . https://doi.org/10.1002/grl.51001 Google Scholar Crossref Search ADS WorldCat Geng J. , Jiang P. , Liu J. , 2017 . Integrating GPS with GLONASS for high-rate seismogeodesy , Geophys. Res. Lett. , 44 , 3139 – 3146 . Google Scholar Crossref Search ADS WorldCat Gérard P. , Luzum B. , 2010 . IERS Conventions (2010) , IERS Technical Note No. 36 . WorldCat Group I.W. , Geodesy G.P.S. , Luca D. , 2016 . Preliminary co-seismic displacements for the August 24, 2016 ML 6, Amatrice (central Italy) earthquake from the analysis of continuous GPS stations , Ann. Geophys. , doi:10.5281/zenodo.61355 . WorldCat Herrmann R.B. , Malagnini L. , Munafo I. , 2011 . Regional moment tensors of the 2009 L’Aquila earthquake sequence , Bull. seism. Soc. Am. , 101 , 975 – 993 . https://doi.org/10.1785/0120100184 Google Scholar Crossref Search ADS WorldCat Huang M.-H. , Fielding E.J. , Liang C. , Milillo P. , Bekaert D. , Dreger D. , Salzer J. , 2017 . Coseismic deformation and triggered landslides of the 2016 Mw 6.2 Amatrice earthquake in Italy , Geophys. Res. Lett. , 44 , 1266 – 1274 . https://doi.org/10.1002/2016GL071687 Google Scholar Crossref Search ADS WorldCat Larson K.M. , Bodin P. , Gomberg J. , 2003 . Using 1-Hz GPS data to measure deformations caused by the Denali fault earthquake , Science 300 , 1421 – 1424 . https://doi.org/10.1126/science.1084531 Google Scholar Crossref Search ADS PubMed WorldCat Lavecchia G. , Boncio P. , Creati N. , Brozzetti F. , 2003 . Some aspects of the Italian geology not fitting with a subduction scenario , J. Virtual Explor. , 10 , 1 – 14 . https://doi.org/10.3809/jvirtex.2003.00064 Google Scholar Crossref Search ADS WorldCat Li X. , Ge M. , Guo B , Wickert J. , Schuh H. , 2013 . Temporal point positioning approach for real-time GNSS seismology using a single receiver , Geophys. Res. Lett. , 40 , 5677 – 5682 . https://doi.org/10.1002/2013GL057818 Google Scholar Crossref Search ADS WorldCat Liang C. , Fielding E.J. , 2017a . Measuring Azimuth Deformation With L-Band ALOS-2 ScanSAR Interferometry . IEEE Trans. Geosci. Remote Sens. , 55 , 2725 – 2738 . https://doi.org/10.1109/TGRS.2017.2653186 Google Scholar Crossref Search ADS WorldCat Liang C. , Fielding E.J. , 2017b . Interferometry With ALOS-2 Full-Aperture ScanSAR Data , IEEE Trans. Geosci. Remote Sens. , 55 , 2739 – 2750 . https://doi.org/10.1109/TGRS.2017.2653190 Google Scholar Crossref Search ADS WorldCat Liu C. , Zheng Y. , Xie Z. , Xiong X. , 2017 . Rupture features of the 2016 Mw6.2 Norcia earthquake and its possible relationship with strong seismic hazards , Geophys. Res. Lett. , doi:10.1002/2016GL071958 . WorldCat Liu J. , Ge M. , 2003 . PANDA software and its preliminary result of positioning and orbit determination , Wuhan Univ. J. Nat. Sci. , 8 , 603 – 609 . https://doi.org/10.1007/BF02899825 Google Scholar Crossref Search ADS WorldCat Liu Z. , Owen S. , Moore A. , 2014 . Rapid estimate and modeling of permanent coseismic displacements for large earthquakes using high-rate global positioning system data , Seismol. Res. Lett. , 85 , 284 – 294 . https://doi.org/10.1785/0220130174 Google Scholar Crossref Search ADS WorldCat Lohman R.B. , Simons M. , 2005 . Some thoughts on the use of InSAR data to constrain models of surface deformation: Noise structure and data downsampling , Geochem. Geophys. Geosyst. , 6 , doi:10.1029/2004GC000841 . https://doi.org/10.1029/2004GC000841 WorldCat Mattei M. , Sagnotti L. , Faccenna C. , Funiciello R. , 1997 . Magnetic fabric of weakly deformed clay-rich sediments in the Italian peninsula: Relationship with compressional and extensional tectonics , Tectonophysics 271 , 107 – 122 . https://doi.org/10.1016/S0040-1951(96)00244-2 Google Scholar Crossref Search ADS WorldCat Melgar D. , Bock Y. , 2015 . Kinematic earthquake source inversion and tsunami runup prediction with regional geophysical data , J. geophys. Res. , 120 , 3324 – 3349 . https://doi.org/10.1002/2014JB011832 Google Scholar Crossref Search ADS WorldCat Melgar D. , Bock Y. , Crowell B.W. , 2012 . Real-time centroid moment tensor determination for large earthquakes from local and regional displacement records , Geophys. J. Int. , 188 , 703 – 718 . https://doi.org/10.1111/j.1365-246X.2011.05297.x Google Scholar Crossref Search ADS WorldCat Melgar D. , Bock Y. , Sanchez D. , Crowell B.W. , 2013 . On robust and reliable automated baseline corrections for strong motion seismology , J. geophys. Res. , 118 , 1177 – 1187 . https://doi.org/10.1002/jgrb.50135 Google Scholar Crossref Search ADS WorldCat Melgar D. , Geng J. , Crowell B.W. , Haase J.S. , Bock Y. , Hammond W.C. , Allen R.M. , 2015 . Seismogeodesy of the 2014 Mw 6.1 Napa earthquake, California: rapid response and modeling of fast rupture on a dipping strike-slip fault , J. geophys. Res. , 120 , 5013 – 5033 . https://doi.org/10.1002/2015JB011921 Google Scholar Crossref Search ADS WorldCat Michelini A. et al. , 2016 . The Italian National Seismic Network and the earthquake and tsunami monitoring and surveillance systems , Adv. Geosci. , 43 , 31 – 38 . https://doi.org/10.5194/adgeo-43-31-2016 Google Scholar Crossref Search ADS WorldCat Scognamiglio L. , Tinti E. , Quintiliani M. , 2016 . The first month of the 2016 central Italy seismic sequence: Fast determination of time domain moment tensors and finite fault model analysis of the ML 5.4 aftershock , Ann. Geophys. , doi:10.4401/ag-7246 . WorldCat Smalley R. , 2009 . High-rate GPS: how high do we need to go? Seismol. Res. Lett. , 80 , 1054 – 1061 . https://doi.org/10.1785/gssrl.80.6.1054 Google Scholar Crossref Search ADS WorldCat Tinti E. , Scognamiglio L. , Michelini A. , Cocco M. , 2016 . Slip heterogeneity and directivity of the ML 6.0, 2016, Amatrice earthquake estimated with rapid finite-fault inversion , Geophys. Res. Lett. , 43 , 10 745–10 752 . Google Scholar Crossref Search ADS WorldCat Tu R. et al. , 2013 . Cost-effective monitoring of ground motion related to earthquakes, landslides, or volcanic activity by joint use of a single-frequency GPS and a MEMS accelerometer , Geophys. Res. Lett. , 40 , 3825 – 3829 . https://doi.org/10.1002/grl.50653 Google Scholar Crossref Search ADS WorldCat Zhu L. , Rivera L.A. , 2002 . A note on the dynamic and static displacements from a point source in multilayered media , Geophys. J. Int. , 148 , 619 – 627 . https://doi.org/10.1046/j.1365-246X.2002.01610.x Google Scholar Crossref Search ADS WorldCat SUPPORTING INFORMATION Supplementary data are available at GJI online. Figure S1. Map-view of slip distributions from strong motion inversion. Observed and predicted horizontal static displacements used in the inversion are plotted with blue and red arrows, respectively. Figure S2. Map-view of slip distributions from 1 Hz GNSS displacement waveforms. Observed and predicted horizontal static displacements used in the inversion are plotted with blue and red arrows, respectively. Figure S3. Map-view of slip distributions from seismogeodetic inversion. Observed and predicted horizontal static displacements used in the inversion are plotted with blue and red arrows, respectively. Figure S4. Map-view of slip distributions from joint inversion. Observed and predicted horizontal static displacements used in the inversion are plotted with blue and red arrows, respectively. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.
TI - Towards the application of seismogeodesy in central Italy: a case study for the 2016 August 24 Mw 6.1 Italy earthquake modelling
JF - Geophysical Journal International
DO - 10.1093/gji/ggy089
DA - 2018-06-01
UR - https://www.deepdyve.com/lp/oxford-university-press/towards-the-application-of-seismogeodesy-in-central-italy-a-case-study-xrcjwC6fmf
SP - 1647
VL - 213
IS - 3
DP - DeepDyve
ER -