Uses of SAR-derived information to support various components and actions of seismic risk management

## 1. Introduction

One of the scientific disciplines in which SAR data have generated a real breakthrough is Solid Earth Geophysics. As soon as the first differential SAR interferograms appeared on scientific journals, describing the surface crustal changes due to earthquakes or inflation of volcanic edifices, geophysicists realized that their science was to greatly benefit from these new data.

The reasons were clear: before the InSAR era crustal deformation was measured using very time-consuming, ground based methods, as leveling, triangulation, trilateration; all requiring costly networks and measurement campaigns. The same was true for the GPS technique, appeared only few years before. SAR interferometry made things simpler and cheaper, although with some limitations; probably the most appealing features of InSAR deformation measurements were their spatial imaging capacities and the high ground resolution. These allowed an unprecedented way to look at all the natural phenomena causing, or derived from, surface deformation.

In this chapter we start describing technical aspects related to the way InSAR data are exploited to derive the parameters of the sources of an observed phenomenon, focusing the attention on the seismic and volcanic activity. We present the state-of-the-art techniques to this inference problem, describing the most common inversion algorithms and analytical models, especially those implemented for the description of the seismic and volcanic cycles. Other technical aspects, as the role played by the data uncertainty, the existing strategies to reduce the data redundancy and the way sources of different nature interact with each other, are also presented to provide useful basics for the reader interested InSAR data modeling.

The second part of this chapter focuses on the practical use of InSAR data and derived models, describing their assimilation in the seismic risk assessment and prevention or during the disaster response. Examples of already existing operational procedures, research and development pilot projects, current and future SAR missions, including Sentinel-1, are also presented to complete this detailed overview on the role played by SAR and SAR-derived information in Solid Earth Geophysics.

## 2. Analytical source modeling

InSAR modeling is an inference process used to define the properties of a realistic geophysical source, or combination of sources, from a set of geodetic measurements. First aspect to clarify is the definition of “realistic”, strongly related to *how* resulting models are exploited. Few hours after an earthquake, the availability of an approximated fault model can be more useful to a decision maker than a complex finite element solution obtained in a week. Conversely, the study of a tectonic setting in a seismogenic area might need a detailed source description, including mechanical discontinuities, the effects of topography, non-elastic rheology and a stratified crust, characteristics not considered by analytical models.

Accepting a simple approximations, as a free surface without topography overlying an homogeneous and isotropic half-space, ordinary geophysical models can be expressed as a set of equations in explicit form, where the modeled surface displacement **dmod** for a given point located in (x,y) is obtained as a function

being *m*_{1}, *m*_{2},… *m*_{n} the model parameters.

Though the above approximations might appear excessive, in most cases they provide very good description of the sources, a fact that is indeed reflected in their massive use in Geophysics. Commonly used sources in the literature are:

**elastic dislocation in a finite rectangular source** [1]: it is definitely the most used model to predict the surface displacement due to an earthquake, represented as a shear dislocation over a finite rectangular fault. The Okada model, however, can be also used to describe magma intrusion like sills or dykes [2],[3],[4], interseismic and post-seismic deformations (see Section 4), landslides [5] and ground subsidence induced by fluid extraction [6]. Source parameters are: East and North position, depth, length, width, strike angle, dip angle, dislocation (or slip), dislocation angle (rake), opening (Figure 1);

**point pressure source** [7]: it is one of the simplest and effective source used in volcanology, as its description requires only 4 parameters: depth, east and north position, volume variation or pressure variation[1] - (Figure 2)

Several other sources have been proposed in the literature, with the aim of providing more realistic solutions to describe geophysical phenomena: dislocation over a finite triangular source [8]; volume variation of a dipping finite prolate spheroid [9]; inflation of an arbitrarily oriented triaxial ellipsoidal cavity [10]; pressure change in a penny-crack source [11]; closed vertical pipe [12]; stress induced by a finite spherical source [13]. A description of the differences among all these sources is beyond the scope of this article, and we refer the reader to the cited literature.

Lastly, we remark the existence of semi-analytical solutions for a layered crust [14] [15], even with a visco-elastic rheology [16]; however they often require a considerable computational time, preventing their use in non-linear inversion schemes, where thousands of forward calculations are often needed. The same limitation affects finite-element models, unless used to build the Green’s function matrix as explained in the next section.

For almost all the listed models, the geometric parameters (position, depth, dimension, orientation, etc…) are non-linearly related to the surface displacement. On the contrary, “intensity” parameters, as the dislocation for the Okada model or the pressure change for a Mogi model, have a linear dependency with the surface displacement. The retrieval of non-linear and linear parameters from geodetic data follows different inversion strategies, explained in the next section.

A further Okada-derived model, widely used to model earthquakes, consists of an array of single Okada sources, or *patches*, for which all the non-linear parameters are already constrained (Figure 3). This compound source is used to reproduce a fault rupture with a variable shear dislocation (or tensile opening), thus improving the modeling performances compared with the uniform slip Okada source.

InSAR measurements are relative, therefore the displacement maps are often related to a reference point considered stable or with a known deformation. This point sometimes turns out to be inappropriate: it could fall in an area affected by atmospheric artifacts; the non-zero displacement field could be larger than the image frame; the GPS value used to tie the reference point might contain a very long wavelength tectonic signal. Furthermore, orbital inaccuracies might introduce artificial linear ramps or even quadratic surfaces. Thus a further non-geophysical source must be often considered, to account for possible improper reference points and/or possible orbital artifacts.

This model can be implemented as a second order surface of the type

so that the apparent displacement **dmod** for a given point (x,y) is a linear combination of the coefficients A, B… F. Second order terms can be neglected when orbital artifacts are well reproduced only with a linear ramp.

We remark that all the geophysical and orbital sources[1] - must always be solved simultaneously in the inversion, otherwise the first source used to model the data will tend to improperly reproduce the signal we might fit with other sources.

## 3. Non-linear and linear inversions

Current approaches to the modeling of geodetic data are the evolution of pioneer studies where the amount of geodetic measurements where certainly lower than nowadays [17][18][19] and the availability of geophysical models were still limited: solutions for dislocation in an elastic medium were available only for some simple fault configurations and for point sources [20][21]. However, non-linear and linear inversion strategies were essentially the same presently adopted. What deeply changed through the decades is the amount of observed data, with InSAR playing a lead role.

An inversion algorithm is a procedure to infer some source parameters so that the modeled, or predicted, surface displacement best reproduces the observed one. We must firstly provide a definition for “best fit”; almost every model shown in literature assumes that the best parameter estimate (**mest**) occurs when the square of the residuals between observed (**dobs**) and modeled (**dmod**) data, i.e. a cost function based on the L_{2} norm, is minimized:

This choice implies that **dobs** data obey to a Gaussian statistics, so that the presence of large outliers is strongly unlikely [22]. However, uncertainties in InSAR data may arise at different stages until the unwrapped and geocoded displacement map is generated, and large spurious artifacts are not so unlikely (see [23] and references therein). Thus the assumption of an L_{2} norm criterion should be taken with care, since it may not be the best choice to build the cost function to minimize.

The inversion strategy must be chosen according to the source parameters** m** linearity with respect to the surface displacement** dmod**. In this section, we assume that the surface displacement is a non-linear combination of the source parameters we want to infer from InSAR data.

Any non-linear inversion scheme is based on the realization of a sequence of forward calculations as in (1) until the condition (3) is satisfied or, in other terms, a set of rules to iteratively change the **m** parameters until the best estimate **mest** is obtained. A peculiar aspect of InSAR data is that displacement measurements are provided in the line-of-sight direction, while every modeled dataset **dmod** is calculated in a Cartesian reference system; therefore the calculation of the L_{2} norm must be preceded by projection of modeled data into the line-of-sight.

A typical issue to face in a non-linear inversion is the presence of an unpredictable number of local minima, corresponding to unsatisfactory solutions; a robust algorithm should be able to take a cost function out from local minima. Unfortunately, only few strategies guarantee the achievement of this goal, for instance Simulated Annealing [24]; however, they are incompatible with a reasonable processing time[1] -. We also remark that the level of non-linearity depends on the parameters to invert: for instance, when periodic parameters like rake or strike angles in the Okada model are fixed or forced to vary in a narrow range, the impact of local minima decreases.

In our experience, the Gauss-Newton or gradient descent methods, as the Levemberg-Marquardt algorithm which is a mix of both [25], provide excellent results with a short processing time; although this algorithm does not consider any stratagem to escape from local minima, its implementation with multiple restarts provides an efficient strategy to identify the global minimum.

We incidentally remark that the repeated calculation of forward models makes unfeasible the use of finite element models in any inversion scheme, because of the long time needed to setup and calculate even a single forward model.

A different approach is adopted when all the parameters to invert are linearly related to the surface displacement, as occurs with the strike-slip, dip-slip and opening components for an Okada model, with the volume variation in a Mogi source or with the coefficients of equation (2). In this case, the inversion can be set up in a matrix form of the type:

where the vector collecting all the parameters **m** is related to the known displacement values **dobs** through the Green’s function matrix** G**.

It is worth noting that the number of unknowns **m** is almost always lower than the number of equations in the linear system (4), given the abundance of InSAR measurements** dobs**. Therefore the problem is generally over-determined and the solution, in the simplest case, is calculated in the least square sense on the data, as follows

being **G**^{-g} the generalized inverse of **G** [22].

However, if we consider the practical problem of finding the coseismic slip distribution as in Figure 6e, solutions found via (4) and (5) are generally not acceptable because of the weak control that InSAR surface measurements have on deep parameters, leading to highly scattered **mest** values. This problem can therefore be treated as partially underdetermined, adding constraints on the model parameters to get a more reliable result. In this case, since the problem is not mathematically underdetermined, the model reliability is obtained at the expense of the data fit. The system of (4) can be arbitrarily extended to introduce almost any type of *a priori* constraint; in earthquake modeling, for instance, we can force the slip to gently vary across the fault and this can be achieved introducing a damping parameter in the form:

where the Green’s functions** G** are extended with a Laplacian operator *ε* [19][26]. In this *damped least square* solution, the *ε* coefficient is problem dependent and must be chosen via trial and error or through the construction of an empirical fit vs. smoothness curve, to find the desired compromise.

Other inequality constraints, such as the positivity (Non-Negative Least Square), can be introduced as well, to further increase the solution reliability; in the modeling of a coseismic displacement field, the constraint **m**>0 is used to avoid unrealistic “back-slip” conditions. For a complete review of the strategies used to add conditions to the linear system (4), refer to [22].

Unlike non-linear inversion, in linear inversions the use of finite elements models (FEM) is not prohibitive in terms of computational time. In this stage, the heavy forward calculation must be done only once to build the kernel matrix **G** [27]. After the Green’s function matrix is available, the problem is reduced to the solution of equations (4), (6) or equivalent.

## 4. Seismic cycle imaging and modeling

Based on the elastic rebound theory formulated by Harry Fielding Reid in 1910, after the 1906 San Francisco earthquake, a seismic cycle is formed by a slow accumulation of stress and deformation, as consequence of the forces acting from the underlying Earth mantle, followed by an impulsive release of stress and energy when the internal strength is exceeded such that the brittle crust breaks (Figure 4). These two phases, which we refer to as inter- and co-seismic, are completed by a post-seismic phase, where different phenomena may induce further deformations during a short- to mid-term period after the earthquake (Figure 5).

The co-, post- and inter-seismic phases have completely different characteristics in terms of duration and crust behavior. During the inter-seismic phase, faults are locked on the upper crust and the underlying forces act to deform the surface with relatively constant rates of few millimeters per year over large areas. As soon as the accumulated stress exceeds the locking frictional forces, the crust cracks and part of the deformation around the faulted area is elastically and permanently recovered in few seconds. Soon after the earthquake, for a period lasting from few minutes to some years, a further deformation occurs as a consequence of the sudden co-seismic stress release. Deformation rates roughly follow a decreasing exponential law (Figure 5) and can be explained in terms of one or a combination of the following factors: residual dislocation on the ruptured fault (after-slip models), viscoelastic relaxation of the lower crust driven by the coseismic stress change, poro-elastic rebound due to the migration of fluids in the crust ([28] and references therein).

This basic description allows to state how SAR derived data can play a crucial role in the understanding of a seismic cycle. In the co-seismic phase, the expected surface displacement, from tens of centimeters to several meters, and the nearly perfect elastic behavior due to the instantaneous deformation are perfect conditions to model the standard two-pass interferometry with the Okada solutions. [30][31][32] showed that since the 1992 Landers earthquake, the onshore deformation for all the significant earthquakes worldwide (M > 5.5) have been imaged with the standard InSAR interferometry and most of them have been modeled with the Okada model. A common problem found in this approach is that, in general, coseismic interferograms contain a contribution from the post-seismic deformation that can affect some fault parameters [31]. Isolating the co-seismic signal may not be straightforward, but the introduction of continuous measurements, such as CGPS, may be helpful in the post-seismic contribution removal.

During the inter-seismic phase, steady deformation rates are often assumed; however, expected values are of the order of some millimeters per year, so they can be hardly detected with two-pass interferometry: the time needed to accumulate a signal above ordinary InSAR artifacts would be too long to preserve the phase coherence; for this reason, time-series techniques like PS [33], SBAS [34], or image stacking are indicated to get mean velocity maps with millimetric accuracy [35], provided that a consistent number of images is available [36][37]. Elastic models have also been used to fit inter-seismic data [38][39][40], however its use must be carefully considered because the assumption of elasticity can lack of realism.

Another important aspect is the apparent similarity between long-wavelength tectonic signals spreading through a whole SAR image frame and the orbital artifacts [41]. In this context, the use of GPS, unaffected by such artifacts, to constrain the velocity maps can be very effective [36][41].

Post-seismic relaxation has intensity and duration strongly dependent on the earthquake magnitude; short intervals can be easily encompassed between two SAR acquisition [42][43], often cumulated with the co-seismic effects. In this case the two contributions cannot be distinguished unless external continuous observations, as GPS, are introduced [44]. For large earthquakes, the post-seismic effects can last from months to years, and the InSAR time-series approach can be effectively used to describe the crustal displacement time evolution [28]. While the observed displacement is interpreted in terms of after-slip over the seismogenic fault, the Okada solutions can be used to model the signal [45]. However, for long-term deformation, visco-elastic models should be used [42].

Regardless from the seismic cycle phase, deformation modeling with the Okada solution is generally subdivided into two steps: a first non-linear inversion to retrieve all the unconstrained source parameters, followed by a linear inversion to get the distribution of the dislocation over the fault(s). The only measure to adopt before running the linear inversion is the widening of the fault plane obtained via non-linear inversion: the latter represents only a mean source with a mean dislocation value, therefore the fault plane must be enlarged to let the slip vanish to zero.

While most of the signal is already reproduced with uniform slip sources, high frequency spatial fluctuations are recovered adopting the Okada-derived distributed slip sources (Figure 6), solved with the linear system (4).

A last remark is about the possibility of obtaining a more realistic source by relaxing the condition of a flat half-space. The possibility of accounting for the topography in the overall system setup will be described in the next section, where this aspect plays an important role.

## 5. Volcanic modeling

The volcanic activity monitoring is conditioned by our ability to describe and understand the eruption cycle. Several steps can be identified in this cycle: magma generation, melting, storage and ascent, crustal assimilation, degassing, crystallization and surface eruption, not all of which can necessarily occur. In any case, evidences of the incoming eruption can be noticed only late in the cycle and InSAR data are crucial to provide information for the hazard mitigation, even for not erupting stages [46].

In the case of volcanic phenomena, several factors contribute to make the assumption of elasticity less reliable than the co-seismic case; magma intrusion starts below the seismogenic crust, thus involving ductile, high temperature layers able to deform aseismically. Furthermore, inflating and deflating phenomena involve times long enough to activate a visco-elastic behavior. Lastly, the half-space approximation is debatable as well, because of the inevitable significant topography of the investigated areas. Despite this limitations, the aforementioned elastic models have been largely applied in volcanology, with the aim of reproducing almost all surface deformations detected with InSAR [47][3][48][4][49][50].

Magma chamber inflation or deflation is generally modeled with the simple Mogi source, due to the ease of its implementation compared with its effectiveness [47][51]. For magma intrusion in vertical (dykes) or horizontal (sills) cracks, distributed opening sources based on the Okada model are instead adopted [48][50].

Sometimes, complex patterns revealed by InSAR suggest the implementation of a multiple-source system, where also seismic sources can have a role, as shown by the 2005 Afar dyking phase described in [52] (Figure 7). In such contexts, the double step non-linear/linear inversions is adopted to first constrain the source geometries, then to retrieve the slip and opening distributions.

For the long-term crustal deformation, time-series techniques have shown their effectiveness to describe the surface change, even when repeated inflation-deflation cycles are present [53]. In this case, modeling can be carried out by fixing the non-linear source parameters and then fitting the time dependent signal by only varying the linear parameter, i.e. the volume or pressure change.

A further remark, in this context, is about the important role played by topography, since strong elevation variations are expected in the investigated areas; for Mt. Etna, for instance, total relief difference is over 3000 m. To mitigate the assumption of flat half-space, topographic corrections can be applied to the analytical models. This correction consists in the calculation of the source depth not from the zero level of the free surface, but adding the real elevation of the point for which the predicted displacement is being calculated. Such compensation has been compared with finite element models, providing a satisfactory improvement on the modeled data [54][55], as shown in Figure 8.

Finally, the frequent presence of stratified atmosphere, altering the interferogram with fringes due to topography-related radar delays, can be discriminated by comparing independent interferograms, assuming weather conditions uncorrelated in time.

## 6. Data downsampling

Since displacement maps derived from InSAR processing may contain millions of valid pixels, with an high degree of spatial correlation, a way to reduce the data must be adopted in any inversion strategy. Several criteria have been proposed in literature, among which we recall here three different approaches: Quadtree decomposition [56][57], resolution-based [69] and regular mesh [58][59].

A general rule to state which is the best method does not exist, though the Quadtree algorithm is the most used. It is a decomposition algorithm aimed at preserving the “amount of information” in the image, and is in general based on the spatial gradient of the signal; it follows that areas with higher displacement values are sampled at higher spatial frequencies (Figure 9b).

The resolution-based algorithm proposed by [69] is driven by an already known source; this allows to calculate the data resolution matrix [22] used to define where surface data must be sampled to constrain the linear source parameters.

Finally, regular sampling is also widely adopted, since its ease of implementation and its effectiveness in imposing a sampling density independent from displacement values. In fact, the InSAR data resolving power, i.e. the maximum detail level achievable on a source, strongly depends on the location of the observed points, as shown in [60] and not on the displacement field itself. The sampling can be manually customized by defining areas with different sampling density; this also allows to have a good control on the number of observed data to handle in the inversion (Figure 9b).

## 7. Uncertainty and trade-offs

InSAR measurements are always affected by different sources of uncertainty, as shown by the ample literature on this topic (see [23] and references therein). Here we discuss the strategies generally adopted to propagate the data uncertainty to the source parameters in the non-linear and linear inversions.

For convenience, linear inversion is firstly analyzed, where ordinary rules for the error propagation can be applied, as:

where the full variance/covariance matrix** Σ**_{d} of the observed data is used as input to find the full variance/covariance matrix **Σ**_{m} of the model parameters. To build **Σ**_{d}, the power spectra analysis of a displacement map containing only noise and artefacts is used. From this analysis, the covariance vs. distance scatter plot, fitted by means of some simple *ad hoc* function (Figure 10), as described in [23], is retrieved. One of the simplest way to express an InSAR covariogram is the exponential function of the type

where C_{0} is the covariance at zero distance, generally lower than the overall *σ*^{2} variance, and *α* controls the distance at which data can be considered uncorrelated. More sophisticated functions can be found in literature, accounting also for a possible spatial anti-correlation [61].

The full variance/covariance matrix** Σ**_{d} can be obtained by setting every off-diagonal position* σ*_{i,j} to the covariance value obtained with (8), setting *r*_{i,j}, as distance between the *i*-th and the j-th point; diagonal values are set equal to *σ*^{2}. After that, equation (7) can be applied.

In the non-linear case, a formal expression of the uncertainty propagation is difficult to obtain and an empirical approach is commonly used: for tens or hundreds of times synthetic noise datasets **dnoise** are generated, added to the observed data** dobs**, then the inversion is performed and the results collected. To construct the **dnoise** dataset, a Cholesky decomposition of** Σ**_{d} must be preliminary carried out, such that

The **L** matrix is then multiplied by a vector **dGauss** of uncorrelated Gaussian noise

to get the synthetic noise dataset **dnoise** characterized by the same InSAR covariogram, as explained in detail in [62]. After the inversions have been carried out, the results can be efficiently visualized as a grid of scatter plots and histograms, describing trade-offs between coupled parameters and single parameter uncertainty (Figure 11).

## 8. Interactions between sources

Models derived from InSAR data give important hints for the hazard assessment, as discussed later in this chapter, and in this respect an increasingly considered aspect is the way sources interact with each other. An earthquake occurs when the internal strength is exceeded by the surrounding stress, loaded during the interseismic phase. We generally do not know the absolute stress value for a given crust volume, primarily because the loading phase spans through centuries or millennia. On the contrary, we can quantitatively calculate the stress variation induced by a fault dislocation.

The stress released during a seismic event perturbs an area extending far from the source itself, where surrounding locked faults are likely to be present. The rearranged stress conditions following an earthquake may increase or decrease the current (unknown) shear stress level acting on a fault surface: we can therefore calculate if such a variation moved the receiver source closer to its failure.

The analytical solutions proposed in [63] allow to calculate the internal deformations induced by a dislocation over a fault plane; internal deformation can be then easily converted into stress variations, using the Coulomb Failure Function variation (ΔCFF), described in [64], and defined as

where Δτ is the shear stress change over the fault, calculated for a given slip direction, *μ* is the friction coefficient, Δσ_{n} is the normal stress variation and Δp the pore pressure change [65]. The latter is usually unknown and thus neglected in stress transfer calculation.

Despite its apparent simplicity, the application of the stress change analysis must be considered with care, because of the intrinsic unknown of the already existing background stress, the mislocation of possible receiver sources, the uncertainty of the triggering source parameters (derived from inversion) and the presence of spatial inhomogeneities. They all introduce a high level of uncertainty, as discussed in [66] and references therein.

The stress variation analysis has been successfully adopted also in a volcanic context, to study the interaction among sources of different nature: magma chambers, dykes, faults [67][68][50]. Though these analyses are always conducted *a posteriori*, after the event occurrence, they contribute to support the stress-transfer reliability.

However the warning issued by [64] is still holding: “much work remains before we can understand the complete story of how earthquakes work”. The CFF analysis is a powerful tool to describe the interaction between sources, but it is still not adequate to deterministically state how close to the failure a receiver source has been pushed by any triggering event.

## 9. From science to operational risk management

The use of SAR data and techniques has greatly stimulated the progress of the Solid Earth science. Globally, over three hundreds deformation fields related to the earthquake cycle (including inter-, co-, and post-seismic deformation), and over a thousand volcanic edifices, have been studied using SAR data since the beginning of the “InSAR era” [30] [31] [32][70][71]. Important new knowledge has been acquired on processes such as fault dislocation, fault segmentation, magma and gas migration, volcanic spreading, stress transfer, strain accumulation and release, poro-elastic diffusion, and visco-elastic relaxation. Better descriptions of the seismic and volcanic cycles are today available thanks to these studies.

Soon after InSAR started to demonstrate its potential to sustain new scientific developments in crustal deformation studies, practitioners started to investigate the possible use of this new information in the risk management activities [72][73]. It was rapidly realized that the new geodetic "imaging" capabilities provided by satellite SAR sensors could strongly support two main components of the disaster risk management, namely the assessment/prevention and the response components.

However, the space and ground segments of the first satellite SAR systems (i.e. ERS, ENVISAT, JERS, ALOS) were targeted mainly to scientific use and, while they could provide delayed data to carry out pre- or post-disaster scientific analyses eventually supporting the hazard assessment component [74], they lacked the near real time capabilities needed for use in the response phase. Even after the launch of the first commercial SAR satellite in 1995 (Radarsat-1) the use of SAR data in disaster risk management did not flourish, due to the lack of a constant repeat pass, global coverage, and to the high costs required to maintain updated archives.

This situation will change radically in 2014, when the European Space Agency Sentinel-1 operational satellite (and its companion Sentinel-1 B in 2015) will start to provide a full InSAR coverage of nearly all land areas with pre-defined, constant repeat pass [75]. The Sentinel-1 data will be delivered in near real time to selected service providers to generate support products during natural disasters or emergencies. Sentinel-1 minimum revisit time will be 12 days, improving to 6 days after the launch of the second satellite, and the mission continuity will be guaranteed for many years [76].

While the "flat" data flow rate and the "full and open data access" policy of the Sentinels will certainly represent a breakthrough for the use of remote sensing data for operational risk management, other SAR satellite constellations have already started to demonstrate the potential for operational emergency response. The Italian Space Agency COSMO-SkyMed four-satellite constellation was in fact expressly developed to support the monitoring and assessment of natural and anthropogenic disasters [77], although it is a dual-use mission, also employed for defense purposes. This constellation presently allows much shorter revisit times than possible with Sentinel-1, down to 1 day depending on satellite.

In the following sections, with reference to Table 1, we will show examples and examine the operational capabilities of InSAR data, and of geophysical models they can constrain, to support activities in the two main components of seismic risk management: the risk assessment/prevention, and the response to a seismic crisis.

### 9.1. SAR-derived products to support seismic hazard assessment

Seismic Hazard Assessment (SHA) is the process of calculating, for a given area, the probability of the occurrence of a certain ground shaking level within a defined period of time. The typical SHA synthetic result for a region is a map showing the spatial variations of the horizontal Peak Ground Acceleration which have a probability of exceedance of 10% within a 50 or 75 year time frame. These maps (complemented by more detailed, site-specific seismic hazard curves) are generated using information on the existing seismic sources, their activity rates and maximum earthquake magnitude, and by estimating the relations between epicentral distance and ground motion for a given earthquake magnitude and fault type (attenuation relations).

InSAR data can provide important information regarding the earthquake source (Table 1). One field of analysis deals with the actual detection of faults by means of their geomorphological signature (e.g. linear scarps, triangular facets on mountain fronts, displaced terraces, drainage network offsets, etc.). This is a classical application of photo-interpretation techniques (structural and geomorphic), in which B&W SAR intensity images give a contribution comparable to optical images. The satellite image analysis can provide especially valuable support to the mapping of active faults and earthquake sources, in areas which are of difficult access and for which detailed geological data do not exist.

While intensity image analysis reveals the fault presence by investigating peculiar landforms created by surface deformation cumulated during geological times, multitemporal InSAR processing can provide a quantitative measurement of the ongoing deformation rates (and their spatial patterns) used to characterize the fault behavior. Many scientific results obtained in the last ten years [32] have contributed precious information for the parameterization of the seismic sources [79][80][81][82][83][84], but also for the definition of the present deformation rates in areas where multiple sources are present [37][85] [86][87][88][89], for the partitioning of strain among different faults [90][91], for the improvement of tectonic models in seismogenic areas [92][93][94].

User needs | SAR-derived supporting information | SAR data analysis techniques | ||

RISK ASSESSMENT AND PREVENTION | Seismic hazard | Identification of active faults. | Structural maps. | Amplitude image analysis. |

Parameterization of activity rates. | Long-term ground displacement rates maps (inter-seismic ground velocity) at low spatial resolution. Fault models and long-term slip rates. | Time-series InSAR techniques as Persistent scatterers and Small Baseline. Non-linear inversion modeling. | ||

Definition of maximum magnitude earthquake. | Fault model geometry and kinematics. long-term slip rates. | Non-linear inversion modeling of deformation data. | ||

Induced hazard | Pre-event identification of gravitational slope deformations. | Geomorphological analysis. Long-term ground displacement rates at high spatial resolution. | Amplitude image interpretation and soil moisture analysis. Time-series InSAR techniques as Persistent scatterers and Small Baseline. Geomechanical modeling | |

Seismic vulnerability | Identification of structural weaknesses in man-made structures. | Long-term ground displacement rates at very high spatial resolution. | Time-series InSAR techniques as Persistent scatterers. | |

DISASTER RESPONSE | Loss estimation | Rapid (1-2 days) spatial assessment of damage to man-made structures | Maps of damage classes at the district scale. Maps of collapsed structures at the single building scale. Very high resolution co-seismic ground displacements. | Intensity-based change detection and classification. Coherence-based change detection. |

Environmental damage estimation | Rapid spatial assessment of environmental effects of earthquake: fault scarps, diffuse ground displacement, reactivated landslides, drainage reversal or interruption, soil liquefaction, sinkhole collapse, etc. | Co-seismic ground displacement maps at high spatial resolution. Geomorphological and structural analysis. | SAR Interferometry, pixel offset tracking, Multiple Aperture Interferometry, intensity image interpretation, coherence analysis. | |

Event scenario | Rapid assessment of earthquake sources and possible evolution of the aftershock sequence. | Co-seismic ground displacement maps at medium-high spatial resolution. Geometrical and kinematic parameters of the earthquake sources. Stress increase on nearby faults. Post-seismic ground displacement for the next few months after the mainshock. | SAR Interferometry, pixel offset tracking, Multiple Aperture Interferometry. Time-series InSAR techniques. Non-linear and linear inversion modeling of deformation data. Coulomb stress analysis. |

Among the most important SAR-derived inter-seismic source parameters directly contributing to hazard estimates is the long-term slip rate. This is the yearly rate of slip which is modeled as occurring on the deep part of an active fault plane, below a given “locking depth”. The latter is defined as the depth separating the upper brittle crust, where the fault plane is locked, from the lower visco-elastic crust where the fault is instead slowly creeping, as explained in Section 4. If the region tectonics is dominated by a single large fault, as often is the case in plate margin contexts, the inversion of InSAR ground velocity measurements can be used to estimate the inter-seismic slip rate at depth.

The fault creep below the locking depth is in turn related to the stress build up in the upper crust, stress which will be eventually released during the seismic dislocation. Considering an ideal earthquake cycle, and assuming that the friction along the fault plane does not vary much through time, each successive fault rupture should occur when the accumulated shear stress overcomes a similar level of fault strength. Thus, the knowledge of the inter-seismic slip rate and of the seismic history of the fault should allow to estimate a recurrence time for moderate to large earthquakes on a given fault. These recurrence interval estimates have in some cases large uncertainties (fault strength and slip rates are indeed not constant through time), nonetheless the inter-seismic slip rates obtained by geodetic data (mainly InSAR and GPS), integrated with paleoseismological and seismological data, are considered important parameters to support seismic hazard calculations [95][96][97][98]. However, while GNSS data have been more extensively used (see for instance the SHARE EC project, the Working Group on California Earthquake Probabilities, 2011; [97]; [99]) to constrain the occurrence models for operational hazard assessment, InSAR data have so far been employed only marginally.

The reason why the inversion modeling of InSAR ground velocities is still not operationally used for the estimation of inter-seismic slip rates and other important fault parameters (as maximum magnitude earthquake), has to do with deficiencies in the data and with gaps in the underlying science. The inadequacy of present SAR system to provide useful data to measure small ground velocities over large spatial wavelengths has been addressed above, and it is expected to be partially resolved by the Sentinel-1 satellites. The scientific issues concern instead the incomplete knowledge of the processes driving the stress accumulation and release in the crust, which imply that the results of the models used to estimate active fault parameters through the inversion of inter-seismic ground velocities are are dependent on scientific judgment. In SHA, the uncertainties arising from the incomplete knowledge of the earthquake processes are addressed by ensuring that the generation of seismic hazard maps is based on discussions within the scientific community, aiming at developing the widest possible consensus on input data, methods, and practices [100]. We expect that in less than a decade the continuous InSAR data flow from the Sentinel-1 operational mission will have promoted the development of new inter-seismic source models, and better procedures for their significance and uncertainty assessment, effectively spreading the InSAR data use in SHA.

Two further important activities in seismic risk assessment/prevention in which SAR data can give a valuable contribution are (Table 1): 1) the identification of gravitational slope deformations which can undergo reactivation or catastrophic collapse during seismic shaking, and 2) the evaluation of the vulnerability level of buildings or infrastructures.

In the first case SAR imagery is used for two different purposes. The classical geomorphic photo-interpretation carried out also on optical images to detect the landforms indicating gravitational slope deformations, can be enhanced by the capacity of SAR intensity images to provide estimates of soil moisture [101]. In general these techniques are used to provide a spatial mapping of the landslide, with little information on its activity [102]. Another use relies instead on high resolution multi-temporal InSAR data to investigate the ongoing rates of gravitational mass movements in high seismic risk areas, which may be used to evaluate more specific hazard levels due to seismically-triggered landslide collapse [103]. The additional hazards induced by landslides reactivated by seismic shaking is very important especially in areas with high topographic relief: in the Wenchuan, 2008 earthquake over 20,000 casualties were attributed to landslide collapse [104].

Finally, a new field of use for InSAR-derived, high resolution ground deformation maps in the practice of seismic risk prevention, is to support a detailed vulnerability analysis in areas affected also by other deformation phenomena with high spatial frequencies (i.e. subsidence, sinkholes). If these phenomena occur in high seismic risk areas, the classification of the severity of the static deformation of man-made structures can allow a more accurate assessment of their resilience to future dynamic actions caused by seismic shaking.

### 9.2. SAR-derived products to support earthquake emergency response

The response phase of seismic risk management concerns all activities needed to promptly respond to the effects of a damaging earthquake. It can be divided in two main, temporally linked sub-phases, the Immediate response and the Sustained response [105].

During the Immediate Response phase, which usually lasts from few to several days, depending on the dimensions of the disaster, the main priorities are search and rescue actions and the emplacement of immediate preliminary measures to save lives and protect the population (evacuate insecure buildings and districts, provide emergency health services and food, install temporary shelters, etc.). A critical element for the management of this phase is the situational awareness, including all information on the extent of the phenomena (the event scenario), of its consequences (the loss scenario), and possibly future developments of both. Satellite SAR data can provide very important information in this phase, although, as we will see, the temporal requirements are difficult to fulfill.

After this initial phase the response actions are directed towards the return to an acceptable state of operation of human activities (repair or reinstall utility networks and infrastructures, provide more comfortable housing structures and working environments, provide temporary social services, etc.). This Sustained Response sub-phase may last from few to many months, and continuous information to monitor this lengthy process is required.

In the Response phase SAR data can be employed to generate two different families of products: those concerning the observation and quantitative measurement of the disaster effects on the human environment, and those providing information on the geophysical processes related to the phenomena (i.e. the earthquake and the triggered effects).

One of the main products of the first family is the damage assessment map, either provided at the district scale or at the single building scale [106]. While damage maps are also generated using data from optical satellites, which can provide higher resolution than SAR systems, the all-weather capability of the active microwave sensors provides an important advantage for rapid mapping in some regions [107]. Both types of data would require pre-and post-event imagery for an accurate change detection (post-event data only have been used, but with degraded accuracy), with SAR data having the additional constraint of the same looking geometry for the two acquisitions (mandatory for coherence-based detection techniques). Unfortunately, for very high resolution SAR and optical data, global coverage is neither continuous nor constant, and high resolution damage maps cannot be routinely generated. However, if ad hoc monitoring actions are started soon after the mainshock, the post-event images can be effectively used to incrementally map the further damage which may be caused by large aftershocks during the next weeks or months.

Where pre- and post-event, same-geometry SAR data do exist, InSAR techniques can provide accurate maps of co-seismic ground displacements. If very high resolution data are available (<3 m) these maps could be used to detect very localized damage to infrastructures (especially large, linear ones). The most common usage however is for the large scale mapping of ground movements directly related to the seismic dislocation (continuous ground displacement field and surface fault scarp), and for the mapping of local phenomena triggered by the seismic shaking (typically gravitational movements). During the Immediate response phase this information is extremely important to develop the situational awareness, for instance to direct rescue teams or implement safety measures; its value is however inversely related to the its delivery time after the mainshock. Depending on disaster size and location, a synoptic satellite map of a fault scarp or of the reactivated landslides may be outdated by more precise aerial or field surveys in 1-4 days. The co-seismic ground displacement map maintains instead its unique capacity to provide a high resolution image of the ground movement patterns which could not be provided by other means, critical for instance for the accurate mapping of gravitational movements, especially when earthquake-triggered accelerated motions are small and do not result in an immediate catastrophic collapse [5].

An additional important element of the situational awareness is the event scenario, which contains an analysis based on detailed information on the various geophysical variables (historical and instrumental seismicity, inter-seismic deformation, co-seismic deformation, local amplification effects, ground acceleration levels, stress transfer levels, etc.) and objects (earthquake source, nearby active faults, geological units prone to landsliding or liquefaction, etc.) which have an impact on the response actions [108].

As we have seen in section 3, the SAR-derived static co-seismic ground displacements are one of the most important datasets to constrain accurate models of the earthquake source [30]. The standard modeling techniques described previously can be implemented to operationally generate source models supporting event scenario development. Several source models (and event scenarios) may then be progressively generated during a seismic crisis, their quality improving as new seismological, geological, and InSAR observations become available and are jointly used to constrain the inversions.

As clearly shown by some recent cases (Emilia - Italy, 2012; Canterbury -New Zealand, 2010-2011, Balochistan - Pakistan, 2008; Umbria Marche - Italy, 1997; Southern California-USA, 1992), aftershocks can sometimes reach higher magnitude levels and cause stronger damage than the mainshocks [80][109][110][111]. As mentioned in Section 8 the triggering of large aftershocks can often be correlated to an increase of stress on nearby faults due the stress released by the mainshock dislocation and redistributed in the crustal volume. The co-seismic fault slip distribution generated by the linear inversion of InSAR data is used to calculate the variations of the Coulomb failure stress induced on the nearby active faults by each large earthquake occurring in a seismic crisis [80][109]. Even if the level of positive stress transfer cannot be used to deterministically quantify the clock-advance of failure on these faults, the knowledge of stress variations can be used to generate quantitative forecasts of aftershock productivity [66], providing important information for risk management during the response phase.

Another information product of interest for the Sustained response phase is provided by the InSAR monitoring of the post-seismic deformation, in particular that generated by poro-elastic diffusion and fault after-slip (Section 4). These ground movements may amount to few tens of percent of the co-seismic ones, and are expressed either as gradual slip increments occurring along the fault plane (at depth or at the surface), and as slow diffuse variations of the co-seismic ground displacement. They could increase the damage levels on man-made structures (e.g. linear rigid structures as viaducts or utility infrastructures), and should be monitored until they become negligible, usually within 7-12 months.

Given the rapid decay of this type of deformation, a dense InSAR temporal sampling is critical for its effective monitoring. For instance, using the full temporal sampling capacity of the COSMO-SkyMed constellation (~7 images per month) after 1.3 months the first time-series can be generated using multitemporal InSAR techniques (minimum dataset: 10 images), while the initial period can be monitored by classical two–pass InSAR. For lower sampling frequencies, the time-series analysis could begin too late, when much of the ground movements have already occurred: using Radarsat-2 for instance (1.25 image/month) 8 months are needed to accumulate a 10-image dataset, but even using a single Sentinel-1 satellite (2.5 image/month) the first deformation time series could be generated only after 4 months.

### 9.3. Provision of risk management services through the SIGRIS system

During the last 15 years, the importance of satellite Earth observation in disaster risk management has been acknowledged also through specific large scale technological programs, as the European Copernicus program (formerly Global Monitoring for Environment and Security - GMES), of which the new operational Sentinel satellites are the pillars. In this framework, an emergency management service has been developed since 2012 to provide mapping products for global-scale natural or man-made emergencies (GIO-EMS, http://emergency.copernicus.eu/). For earthquake emergencies GIO-EMS can be activated to provide rapid information products, as reference maps and damage assessment maps generated using optical and SAR data, to institutional users and national Civil Protection bodies.

A wider range of products, including all those described above, was demonstrated by the Italian SIGRIS system [108]. This is an Earth Observation monitoring system developed to generate information products based on various types of satellite imagery, to support the management of the seismic risk. The system development was promoted by the Italian Space Agency to exploit the potential of the Italian COSMO-SkyMed SAR satellite constellation, although other SAR and optical data are also used. The system requirements were provided by the Italian Civil Protection, which is presently the main user of the SIGRIS services. The system is now maintained and operated by INGV, a leading Italian geophysical research institute, which manages all the national ground-based networks for the monitoring of earthquake and volcanic phenomena and is also responsible for the production of the national hazard maps.

The SIGRIS system was conceived to provide both Assessment/Prevention and Crisis Response services, in support of various activities of the National Civil Protection Service, of which INGV is an integral part. It uses state-of-the-art InSAR and optical data processing and geophysical modeling algorithms, as well as validation/verification, reporting, and dissemination procedures, most of which are executed through a GIS-based interface.

The SIGRIS Assessment/Prevention service implemented the processing chain to generate most of the SAR-derived information products described in section 9.1, and the products were demonstrated in several test cases. Examples of the results derived from these products are shown in Figures 12 and 13.

The demonstration of the SIGRIS products for Assessment/Prevention confirmed that further improvements in the SAR data and in the geophysical modeling capacities (as mentioned in section 9.1) are needed before robust results can be generated on a routine basis. The new data provided by Sentinel-1, as well as by the ALOS-2, L-band SAR system, will be an important step towards both directions.

For the Response phase, SIGRIS contains processing chains and procedures to generate, validate and deliver the information products described in section 9.2. During the development phase, SIGRIS products were demonstrated for different areas of the world (www.sigris.it), but since 2011 the system is mainly activated for national earthquake emergencies.

Very positive has been the user evaluation of the quantitative, InSAR-based assessment of the co-seismic and post-seismic SIGRIS products, and especially for the earthquake source models. The system was used to obtain timely information products for the response phase of the L'Aquila, 2009, Emilia, 2012, Pollino, 2012, and Lunigiana 2013, Italian earthquakes. The response products were delivered in successive versions with increasing information content. Usually the initial source models are based on fast, standard procedures for inversion and uncertainty assessment [59], while later models are constrained by a larger number of datasets and may involve the use of non-standard modeling techniques [27].

For all the mentioned events the SAR data were instrumental to the definition of the earthquake source, as GPS data were limited (L'Aquila) or almost completely missing (Emilia, Pollino, Lunigiana), and provided minimum constraints for the modeling. Acknowledging the usefulness of COSMO-SkyMed InSAR data for Italian emergencies, in 2010 a routine acquisition plan was devised by the Italian Space Agency. The MapItaly plan is now operative to cover all of the Italian territory with a new acquisition every 16 days, providing the necessary, continuously updated archive, for all InSAR applications.

## 10. Conclusions

The use of SAR data has now become common practice among geophysicists involved in the monitoring and understanding of Solid Earth phenomena. By far the most important use of SAR data is for the measurement of surface ground displacements and for constraining models of crustal deformation.

Various scientific and commercial software packages are available for the data processing and for the modeling of the ground displacements, and in the last years a particular attention has been given to the integration of displacement data from InSAR and GPS, to cope with some of the inherent limitations of SAR systems, and augment the number of applications.

After the termination of ENVISAT-Asar, and ALOS-Palsar, in 2011, geophysical applications have suffered the lack of continuously updated archives, which the other national missions (Radarsat-2, TerraSAR X, COSMO-SkyMed) do not provide routinely over large areas. Great expectations are placed in the Sentinel-1 operational mission, which will provide at least two decades of high quality SAR data with constant and improved repeat pass over most of the emerged lands.

The enhanced characteristics of the Sentinel-1 system, as the larger swath, more precise orbit control, shorter repeat pass, improved resolution, open data access, will provide better and free data for all, increasing the diffusion of SAR data use in Solid Earth Geophysics. Thus the next decades are bound to see new science, better interpretation methods, and more effective operational applications based on Synthetic Aperture Radar.

## Notes

- The Mogi equations can be indifferently written in terms of pressure change or volume change; however, the volume/pressure conversion requires the definition of a virtual source radius.
- We adopt the misleading definition “orbital source” only because it is mathematically handled just like the Okada or Mogi equations.
- Simulated Annealing is largely used in non-linear inversions due to the ease of its implementation, but never with the cooling schedule that guarantees to find the global minimum, that would require a nearly endless computation time.