Spatial forecasting of seismicity provided from Earth observation by space satellite technology

Understanding the controls on the distribution and magnitude of earthquakes is required for effective earthquake forecasting. We present a study that demonstrates that the distribution and size of earthquakes in Italy correlates with the steady state rate at which the Earth’s crust moves. We use a new high-resolution horizontal strain rate (S) field determined from a very dense velocity field derived from the combination of Global Navigation Satellite System (GNSS) and satellite radar interferometry from two decades of observations. Through a statistical approach we study the correlation between the S and the magnitude of M ≥ 2.5 earthquakes that occurred in the same period of satellite observations. We found that the probability of earthquakes occurring is linked to S by a linear correlation, and more specifically the probability that a strong seismic event occurs doubles with the doubling of S. It also means that lower horizontal strain rate zone can have as large earthquakes as high horizontal strain rate zones, just with a reduced probability. The work demonstrates an independent and quantitative tool to spatially forecast seismicity.

At seismically active plate boundaries, forecasting the distribution and sizes of earthquakes is fundamental for seismic hazard assessment. Approaches for seismic hazard estimation were traditionally based on forecasting the probability of future earthquakes from the statistical analysis of historical and instrumental earthquake catalogues 1 . However, this approach is flawed in most regions of continental deformation since the recurrence interval for earthquakes exceeds the length of the historical seismic record 2 . Instead, an alternative and independent approach in forecasting earthquake activity in regions that are deforming relatively fast involves using geodetic measurements of interseismic strain-rate as a proxy for future seismicity 3 . The theory of elastic rebound forms the basis of this method by considering that the elastic potential energy budget in the seismogenic crust is released by seismic events. With this assumption, the rate of elastic strain accumulation should correlate with the rate of strain release by earthquakes. Such an approach has successfully been implemented qualitatively for the San Andreas Fault system of California 4 and North Anatolian system of Turkey 3 . Since regions of high strain rate should theoretically rupture more frequently, a quantitative correlation between increased frequency of earthquakes of any magnitude and strain rate should exist even for an earthquake catalogue that does not sample the full seismic cycle.
We test whether the rate of seismicity can be forecast from interseismic strain-rate by using the distribution and magnitude of crustal earthquakes reported in the Istituto Nazionale Geofisica e Vulcanologia (INGV) earthquake catalogue 5 during 1990 to 2017 complete above magnitude (M) 2.5 ( Fig. 1). We also compute a new and high-resolution strain rate map of Italy using geodetic data from the same time period. To do this we compute ground velocities using Global Navigation satellite Systems (GNSS) permanent stations and campaigns [6][7][8] , and integrate with ground velocities from Persistent Scatterers (PS) Interferometric Synthetic Aperture Radar (InSAR) 9 . The merging of GNSS and PSInSAR 10-13 is used to produce an improved map of the horizontal strain rate (S), that quantifies the change in horizontal strain (deformation) of the Earth's crust. The excellent seismicity catalogue combined with the well constrained and relatively fast strain rates of Italy make our data ideal to quantitatively test the correlation between position and size of earthquakes and S.
Tectonically, Italy is on the plate boundaries between the colliding African and Eurasian plates, that has been forming the Alpine and Apennine belts since around 90 Myrs ago 14 . Seismological data, recent geodetic

Results
Seismicity catalogue and strain rate computation. We used seismicity in Italy from January 1990 to December 2017 reported by the INGV catalogue, which is complete above M 2.5. We first isolated 24255 earthquake of magnitude M ≥ 2.5 and those that occurred in the crust using the Moho depth map reported by Labrousse 22 . Earthquakes typically have horizontal uncertainties in their position of less than 5 km. The size, low magnitude of completeness, and well constrained earthquake locations makes the INGV it one of the most comprehensive earthquake catalogues for a tectonically active region on Earth.
The map of strain rates of Italy was derived from a geodynamic velocity field of Italy calculated over the period 1994 to 2014. To compute the velocity field, we integrated the velocities of 621 GNSS sites with 25 million PS data points 13 . The resultant single and very dense velocity field is then gridded at 20 × 20 km and used to compute the two-dimensional velocity gradient tensor (Fig. 1B), from which the second invariant of the horizontal strain rate is derived (see methods section for the full description of how S maps were derived).
With our earthquake catalogue and S maps we use a Bayesian, or conditional approach, to determine the probability function for the occurrence of an earthquake of a given magnitude within a given S. From the analysis that we describe below we are able to determine that the seismicity rate increase with S by a linear trend. We show that the probability of an earthquake of any magnitude is highest at high strain rates following the Gutenberg-Richter law 23,24 . Probability density functions of earthquake magnitude and horizontal strain rate. The first step in our analysis is to compute the probability density function of magnitudes in the earthquake catalog from the Gutenberg-Richter distribution; where N is the cumulative number of earthquakes ( Fig. 2A) with magnitude (M) higher or equal to the magnitude of completeness (m c ) which is the lowest magnitude at which the seismic-event catalogue is complete. The parameters a and b are respectively the earthquake productivity of a volume and the relative size coefficient distribution 25 The second stage of the analysis was to determine the distribution of S. The distribution of horizontal strain rates of the study area d S ( ) was determined from the normalized histogram of S (Fig. 2D,E). We hypothesise that S follows the Weibull distribution (Fig. 2F): www.nature.com/scientificreports www.nature.com/scientificreports/ Correlation between seismicity and horizontal strain rate. To correlate seismicity and horizontal strain rate (S) we associated to each earthquake the strain rate of the 20 × 20 km grid cell in which the earthquakes occurred. We therefore obtained a dataset in which for each earthquake we know the earthquake magnitude and strain rate. The magnitude of earthquakes and the S are plotted against each other in Fig. 3A. Since we had to determine the mutual condition probability f S M ( ), which is the probability function of the strain rate S given an earthquake of a particular magnitude, we calculated the distributions of f S M ( ) for different classes of magnitudes and display the results in Fig. 3B. The results show that the distribution of strain rates is stochastically independent of the magnitude of events. This means that earthquakes grouped into classes of magnitude have the same probability density as a function of the strain rate S. In addition, the probability of earthquakes occurring increases only with the strain rate S and we can therefore We interpret that f S eqk ( ) is a linear trend until the maximum value of the strain rate s max of the study area which is ~40 nstrain/a (Fig. 3B): Finally, we determined the probability that a seismic event of magnitude M occurs in an area with a given S as a combination of independent events: www.nature.com/scientificreports www.nature.com/scientificreports/ Replacing equations Eq. 2 and Eq. 4 in Eq. 6: The resultant Eq. 7 preserves the Gutenberg-Richter law for the magnitude distribution of earthquakes. The results from applying Eq. 7 shows that the probability of earthquakes of any magnitude is linearly dependent on the strain rate (Fig. 3C).

Discussion
The linear dependency that links the occurrence of earthquakes of a particular magnitude with the increase in strain rates shows that regions of high inter-seismic strain rate are more likely to be seismically active across all earthquake magnitudes (Fig. 3A). This is consistent with the hypothesis that earthquakes are the sudden release of accumulated elastic strain. Our study is the first to statistically show a positive link between inter-seismic strain rate and the occurrence of earthquakes across a range of earthquake magnitudes. The relationship between strain rates and seismicity is also visually clear in map view. In Fig. 4 we plot the map of ∩f M S s ( ) coefficient  Our work advances efforts to link seismic activity to interseismic strain rates. Previous studies from the San Andreas Fault in California have demonstrated that regions of high strain rate correlated with the occurrence of relatively large (M 6.5) earthquakes but with lower magnitude events showing a more complicated pattern 4 . Our results also differ from a previous analysis of the link between strain rates and seismicity in Italy 27 , which showed a statistical link between the occurrence of the largest magnitude earthquakes in Italy with low to moderate strain rates. Their study suggest that regions of locked faults have relatively lower inter-seismic strain rates but rupture in relatively larger events. In contrast to both these studies, we demonstrate that the probability of an earthquake occurring of any magnitude is highest in regions of elevated strain rate. Therefore, the results support the simple hypothesis that fault zones follow simple elastic rebound theory and exhibit Gutenberg-Richter behavior.

conclusions
We use a very dense velocity field derived by the integration of GNSS and InSAR over 20 years, and the INGV earthquake catalogue complete above magnitude 2.5, and demonstrate the relationship linking earthquake magnitude with horizontal geodetic strain rates. We found the probability of seismic events at all magnitudes linearly increases with the strain rate. Our approach quantitatively demonstrates that measuring the rate of inter-seismic strain accumulation is an effective method for forecasting the distribution and magnitude of earthquakes.

Methods
Geodetic velocity maps and strain rate computation. The map of the strain rate of Italy was derived from a geodynamic velocity field of Italy calculated over the period 1994 to 2014 by the integration of the velocities of 621 GNSS sites (Fig. 5A) with 25 million of PS 6 (Fig. 5B). PSInSAR is a technique applied to satellite SAR images that produce persistent scatterers (PS), which are sparse ground point-wise radar benchmarks characterized by long-term stability of the electromagnetic backscattered signal and high reflectivity [28][29][30] . This accurate and very dense velocity field is the base for the determination of a correct geodetic strain field. Full details of the geodetic dataset are reported in Farolfi et al. 13 .
The two-dimensional velocity gradient tensor is calculated by applying the infinitesimal strain approach 31,32 . Due to the fine scale detail of horizontal velocity, and to reduce the size of the large dataset, we have mapped the velocity field and strain tensor with a grid of 20 km × 20 km. The known horizontal incremental velocity vector V i of the i vertex polygon is defined as: where A i is the unknown velocity at the origin of the coordinate system, x j is the position of the station, and t ij is the displacement gradient tensor. Tensor theory states that any second-rank tensor can be separated into a symmetric and an anti-symmetric tensor then t ij can be additively decomposed asfollows: For infinitesimal strain rates, the symmetric and anti-symmetric parts can be associated with the infinitesimal strain e ij and rotation ω ij tensors. Principal strains (Fig. 6) were computed as: