Electronic Supplement to
A Two-Fault Model of the 2003 Leucas (Aegean Arc) Earthquake Based on Topological Inversion of GPS Data

by Vasso Saltogianni and Stathis Stiros

Analytical Inversion Methodology—The TOPINV algorithm

Physical Model

Following the Okada (1985) model, faults are approximated by rectangular surfaces with uniform slip in an elastic half-space and are described by $n=9$ scalar variables (parameters, i=1, 2,..., n), summarized as vector x: easting (E) and northing (N) of the center of a rectangular fault with its upper edge parallel to the ground surface; depth of the upper edge (d); strike, dip, and rake angles; length (L) and width (W) of the inclined rectangular fault; and amount of slip s. Surface displacements due to fault activity are related to the fault parameters using the Okada equations.

Basic Algorithm

The TOPological INVersion algorithm can be analyzed in four steps.

1. A system of observation equations is formed. In symbolic form, for an Okada-type fault x, observed Global Positioning System Cartesian displacement components (east, north, up) at a point lead to two or three equations for $f_j$ (j=1, 2,..., m) of the form $$f_j(\mathit{x})-l_j={\upsilon }_j,\mathrm{(S1)}$$ where ${\upsilon }_j$ is an unknown error (misfit) to account for imperfections in the geophysical (Okada) model $f_j$ and for the error in measurement $l_j$ for equation j. Hence, for a single fault and m observations of surface displacement in total, a system of m equations with $n=9$ unknowns is formed. Usually, redundant observations are available; and, because equations (S1) are highly nonlinear, this system of $n=9$ unknowns for a single fault and $n=18$ unknowns for two faults cannot be solved using conventional algebraic techniques. If the physical model (Okada faults) is assumed realistic and outliers are excluded, the absolute value of each ${\upsilon }_j$ would correspond to a bounded set of real numbers, and these bounds can be approximately assumed proportional to the standard error ${\sigma }_j$ of observation j. Hence, we may define a variable $k>0$ so that $$|f_j-l_j|=\left|{\upsilon }_j\right|<k{\sigma }_j,$$ which yields the observation inequalities $$|f_j-l_j|<k{\sigma }_j\mathrm{.}\mathrm{(S2)}$$
2. A searching area, grid G, in the $R^n$ space is defined on the basis of seismological or other geophysical constraints. For each of the $n=9$ or 18 variables defining a single- or double-fault Okada model, the possible range of the values of each variable is defined on the basis of seismological or geological evidence, and a spacing between them is selected. This transforms the continuous hyperspace into discrete points, and, by definition, at least one of these gridpoints P is close to the unknown solution(s) of the system of equations, if any exists in the space defined by grid G (i.e. if our geophysical assumptions/constraints are correct). So the grid G will be a close set of gridpoints in $R^n$ space, defined by the condition $$x_i^{\min }\le x_i\le x_i^{\max }\mathrm{.}\mathrm{(S3)}$$ And, if $p_i$ gridpoints correspond to variable $x_i$ , grid G will consist of p gridpoints, $$p={\mathrm{\Pi }}_i{p}_i\mathrm{.}\mathrm{(S4)}$$ Grid G can be regarded also as a matrix with n columns and p rows. Each row corresponds to a candidate fault x or a possible solution of the system of equations according to the initial hypothesis (Saltogianni and Stiros, 2012).
3. The algorithm does not try to identify the unknown gridpoints P closest to the unknown solution(s); that is, it d.oes not focus on maxima or minima, but only identifies one or more sets (clusters) of gridpoints that define closed spaces S containing P. This process permits identification of multiple solutions and is based on the scanning of all grid points of grid G in order to identify those that satisfy each of the inequalities (S2) using forward computations and a Boolean test (yes/no). Then the final solution S, which satisfies all the inequalities (shown in S2), is computed as the common section of the sets S_j; that is, $$S={\cap }_j{S}_j\mathrm{.}\mathrm{(S5)}$$ Set S should correspond to a closed n-dimensional space. If not, space S is spread in more than one closed space, each corresponding to a different solution (Saltogianni and Stiros 2013a) The process tends also to optimize S (i.e., to identify S*, which corresponds to the minimum value of k=k*). In practice, at first a value of k around 1 usually is adopted, and then the process is repeated for larger or smaller values of k until the minimum value of k=k* defining an optimal set (cluster) of gridpoints $S^{*}\subseteq G$ is reached. The term “optimal” indicates sets S* with at least two distinct values for each variable; this is necessary to calculate variances and covariances matrices (VCM). This is an essentially deterministic approach, different from sampling-based approaches (Sambridge and Mosegaard, 2002) and from those adopted in other cluster-search methods in gridded hyperspaces (Tasoulis and Vrahatis, 2006), and the only stochastic input is the standard errors of measurements playing the role of weights in analogy to conventional least squares adjustment.
4. The algorithm computes the center of gravity of set (or sets, clusters, closed spaces) S*. This last step leads to a stochastic minimum variance solution of the unknown vector x; such solutions depend only on the initial conditions imposed by grid G and the adopted error weights of observations. From the calculated center of gravity in the Rⁿ space, the full variance covariance matrix (VCM) of the solution can be obtained (i.e., the first and second statistical moments). A weighted mean misfit error ${\chi }_{\upsilon }^2$ between data (observations) and predicted data can also be computed.

Nested Grids

A practical limitation of the algorithm is that the number of p rows should be kept in a limit imposed by the computer capacity. In order to avoid an excessively large grid G, grids with less than p=10¹⁰ gridpoints are usually used for common computers (see Table 1 in main article).

For this reason, the analysis is made first for a large, coarse grid (large range and separation for values of each variable), and it is then repeated for essentially nested, finer grids (with smaller ranges and separation for each variable), permitting more precise and accurate approximation of the sets S*, which include the solutions of the system of equations (S1).

The overall process is explained in Figure S1 for two selected variables, the strike and the dip of the two segments for model 1 (described in detail in the main text). The a priori conditions adopted for both segments require the strike to be in the 0°–30° range and the dip in the 50°–90° range, and a spacing of 10° was selected to produce a first, coarse grid G1. Application of the algorithm permitted to discard some of the gridpoints (those marked by open circles) and select those with solid circles. In fact Figure S1 reflects the projection of hypergrid G1 onto two planes, each defined by the dip and strike of the two segments. Subsequent application of the same process with finer grids around the selected points (grids nested to G1, not shown for simplicity) led to the final grid G5 and its projections superimposed on those of G1. Points of G5 marked with open triangles were rejected, and those marked with solid triangles were selected as the cloud of points defining the optimal solution (projection of gridpoints S* on two planes).These points (solid triangles) represent ~500 in total points in the R¹⁸ space, surround and approximate the “true” solution and from their first and second moments, the best estimates of the dip and strike and their variances and covariances shown in Tables 2 and 3 in the main article were computed.

Instead of a nested grid, a solution may indicate the need of an extended grid. For example, in Figure S1, if the solution for G1 showed one single point at strike 10° and dip 50° for one of the segments (i.e. at the margins of the grid), it would be necessary to extend the grid to cover the possibility of dip 40°, if no other restriction exists.

The above approach is implemented into an algorithm using FORTRAN coding.

Validation

The accuracy and the possibility of missing values and of false alarms of the algorithm have been assessed using synthetic data (a “truth-based approach”). Based on a selected geophysical model, a vector x was assumed (i.e., fault characteristics, etc.) and observations of displacement (synthetic observations) were produced. Then the synthetic observations (or synthetic observations contaminated by white noise) were inverted on the basis of the same geophysical model, and an estimation $\hat{\mathit{x}}$ of vector x was computed. Finally the difference between x and $\hat{\mathit{x}}$ was examined and found statistically non-different (Saltogianni and Stiros 2013b; Stiros and Saltogianni 2014).

The reasons for the successful inversion are that the adopted algorithm is based on an objective approach, without any need for a priori fixing any of the variables at any step of the analysis. Input in the analysis were the geodetic observations and their uncertainties, while the adopted search space G for the unknown fault variables (parameters) was designed to exclude physically non-reasonable values of the variables (parameters). The only relatively subjective parameters introduced were those defining grid spacing, which, however, was compatible with the uncertainties of observations (Table 1 in main article). In addition, as can be derived from Figure S1, the distribution of the “successful” grid points used to obtain models 1 and 2 in relation to the whole grid was satisfactory from the statistical point of view (see main article).

---

Figure

Figure S1. Projection of hyper-grid G1 on planes (2-D spaces) defined by the strike and the dip of the Cephalonia and Leucas segments for model 1 (see main article for details). Gridpoints are marked by circles and represent all a priori possible solutions of the system of inequalities (S2) using grid G1 (Table 1 in main article). The algorithm led to the discarding of points marked by open circles, and the solution for grid G1 is marked by solid circles, corresponding to ~200 gridpoints in the R¹⁸ space in total. Repeating the same process using grids nested to G1, it became possible to obtain a solution derived from a final grid G5, indicated by triangles. The final solution is marked by solid triangles, which represent ~500 gridpoints in the R¹⁸ space. Statistics of these points (first and second moments) lead to the optimal solution and its variance–covariance matrix.

---

[ Back ]