Electronic Supplement to
Fault Model for the 2015 Leucas (Aegean Arc) Earthquake: Analysis Based on Seismological and Geodetic Observations

by Vasso Saltogianni, Tuncay Taymaz, Seda Yolsal-Çevikbilen, Tuna Eken, Fanis Moschas, and Stathis Stiros

This electronic supplement demonstrates additional information and results of the main article. In particular, it contains a detailed documentation of the methodology adopted for the inversion of the geodetic observations, a table of the Global Positioning System (GPS)-derived displacements and figures of P-wave first-motion polarities, a comparison of earthquake source parameters, GPS time series, 2D projections of geodetic solutions, the variance–covariance matrix of the geodetic solution and the geodetic variable slip model.

---

Inversion of Geodetic Observations—The TOPINV Algorithm

Geodetic observations of slip and elastic dislocation equations lead to redundant, nonlinear systems of equation of the type

in which i = 1, 2,…, n, j = 1, 2,…, m, n < m, x_i the unknown variables, ℓ_j the observations, and υ_j small errors in observations with standard deviations σ_j, or symbolically

Such systems of equations cannot be solved using algebraic techniques (least squares, etc.), and for this reason the solution is based on sampling-based techniques and optimization on the basis of a cost function. In the most elaborate cases (see Cervelli et al., 2001; Sambridge and Mosegaard, 2002; Menke, 2012), an n-dimensional search space is first defined (n is the number of unknown variables), and a small sample, usually consisting of 100–1000 points from the search space, is selected. These algorithms then search for a single point ${\hat{\mathbf{x}}}_{\mathbf{o}}$, for which the misfit between predictions and observations is minimum in simplest form

This point is regarded as the solution of the system of equations. A priori estimates of variables (fixed variables) to reduce the number of unknowns are sometimes introduced.

The TOPological INVersion (TOPINV) algorithm is based on a different strategy for a simultaneous inversion for all the n unknown variables. At first, the possible values (search space) for all n unknown variables are used to construct an n-dimensional search space, which is approximated by a (hyper-) grid G. Then, the observation equations (equation S2) are transformed into inequalities (equation S4), which are functions of a scalar optimization parameter k, with σ indicating symbolically standard deviations playing the role of weights.

An initial value of k is first selected and all grid points are tested uniformly (scanned, exhaustive searching) whether they satisfy all observation inequalities (equation S4). This is a simple Boolean test, and for each grid point the output is yes (successful points) or no, because these inequalities contain no unknown variables. This analysis leads to a cluster (set) of successful grid points that define an n-dimensional space containing the true solution of the system.

If this cluster corresponds to a closed, convex n-dimensional space, it corresponds to a single solution, which is computed using a conventional stochastic approach: calculation of the first statistical moment (solution) and second statistical moment (its variance–covariance matrix); a mean misfit is also computed. On the contrary, if this cluster does not correspond to closed, convex space, it corresponds to two or more solutions, and it is spread in two or more parts/spaces each leading to different solution. This process is repeated for different (usually gradually smaller) values of k until an optimal (minimum variance–minimum misfit) stochastic solution is obtained. It is important to notice that uncertainties (variance–covariance matrix) are obtained from the analysis of a whole population of grid points leading to the solution as first and second statistical moments (Mikhail, 1976), whereas in sampling-based approaches uncertainties derive from sampling in the vicinity of the solution (Cervelli et al., 2001).

The overall analysis has a high computational cost, and there is an upper limit for common computers (usually total number of grid points of the order 10¹⁰). To overcome this problem, at first a broad search space with smaller density (grid G1) is adopted, and the algorithm identifies subspaces that contain the solution(s). Then, around these subspaces (coarse, preliminary solutions) nested grids G2, G3, etc., with smaller coverage but higher density are selected, and the algorithm is repeated until an optimal solution is obtained.

The analytical description of the TOPINV algorithm, adopted in the inversion of the geodetic data in the case of the 2015 Leucas earthquake, is more explicitly documented in the electronic supplement of Saltogianni and Stiros (2015).

---

Table

Table S1. Computed and modeled coseismic displacements of the Leucas 2015 earthquake.

---

Figures

Figure S1. Distributions of P-wave first-motion polarities of the 2015 Leucas earthquake. Up and down P-wave polarities are indicated by solid and open circles, respectively. Epicentral distance (Δ) of each station is also shown.

Figure S2. Comparison of selected observed waveforms (solid lines) with synthetic ones (dashed line) computed based on our preferred long-period minimum misfit solution of the Leucas earthquake (top row) with the source parameters reported by the present study with focal depth to be fixed at 5 and 12 km, respectively (second and third rows), Sokos et al. (2016) (fourth row), Global Centroid Moment Tensor (Global CMT; fifth row), and U.S. Geological Survey (USGS)-W phase (sixth row). The stations are identified at the top of each column, with the type of waveform marked by P and SH and followed by the instrument type. At the start of each row are the P- and SH-focal spheres for the focal parameters represented by the five numbers (strike, dip, rake, depth, and seismic moment), showing the positions on the focal spheres of the stations chosen. X shows matches of observed to synthetic waveforms that are worse than in the minimum misfit solution (see also Fig. 2 in the main article).

Figure S3. Time series of coordinate changes of permanent GPS stations during the day of the earthquake (17 November 2015), relative to arbitrary selected origins; kinematic precise point postioning (PPP) solutions using 30 s data. Part of the analyzed records, ~15 hrs long, is shown. A homogeneous pattern in the fluctuations of the coordinates is observed. The mainshock and an aftershock of magnitude 5.0, also associated with some smaller-scale displacements, are indicated by dashed lines. Shading indicates the part of the time series used to compute postseismic displacements; a red line indicates the average preseismic coordinates derived from 7.5-hr-long records. The location of the sites is shown in Figure 1 of the main article.

Figure S4. Daily time series of coordinate changes of stations KIPO and ASSO from late 2010–2015. Red lines indicate the trends that derived from the analysis of the available 3-yr-long record (Fig. S5). Coseismic displacements were computed as differences between trend-predicted coordinates before and after the earthquake.

Figure S5. GPS velocities in International Terrestrial Reference Frame (ITRF) system for the time interval 2011–2013 for the (a) horizontal and (b) vertical component. Black arrows correspond to stations that were used to estimate the velocities in station ASSO (red arrow).

Figure S6. 2D projections of ~170 grid points (small dots) corresponding to the refined solution. The 36 panels correspond to projections of the n-dimensional grid G1 into all possible combinations of 2D planes. Embedded rectangles correspond to the area of the grid G2 (Table 1 in the main article), and the large red dots indicate the preferred (refined) high-angle geodetic solution. Insets in top indicate enlarged the panels of strike–dip and dip–rake with superimposed the seismological point source solution and its error bars.

Figure S7. Graphic representation of the absolute values of the elements of variance–covariance matrix for the final solution representing in Figure S6. Cross correlations are small, because of the simultaneous inversion in nine variables, highly reducing the trade-off characterizing other types of inversion.

Figure S8. Variable slip model of the geodetic model of Figure 5 of the main article. The arrows denote the direction of slip on the surface (i.e., the rake angle). A maximum slip patch of 250 cm appears at the northern area of the fault at ~5 km depth. Original solution not corrected for zero slip at fault borders.

Figure S9. Alternative solutions for the fault model derived from grid G1 (Fig. 4 and Table 1 in the main article) for k = 2.5. The 14 grid points shown in the panels (all 36 possible combinations of the 2D projections of the grid points) correspond to different low-angle fault models, not further analyzed in this study.

---

[ Back ]