ARTICLE IN PRESS

Computers & Geosciences 34 (2008) 967–977

www.elsevier.com/locate/cageo

ISOLA a Fortran code and a Matlab GUI to perform

multiple-point source inversion of seismic data$
Efthimios N. Sokosa,, Jiri Zahradnikb
a
Laboratory of Seismology, Geology Department, University of Patras, Rio 26504, Greece
b
Faculty of Mathematics and Physics, Charles University in Prague, Czech Republic
Received 22 December 2006; received in revised form 3 May 2007; accepted 27 July 2007

Abstract

In this paper, a software package for multiple- or single-point source inversion is presented. The package consists of
ISOLA-GUI, a user-friendly MATLAB-based interface, and the ISOLA Fortran code, which is the computational core of
the application. The methodology used is similar to iterative deconvolution technique, often used in teleseismic studies, but
here adjusted for regional and local distances. The advantage of the software is the graphical interface that provides the
user with an easy to use environment, rich in graphics and data handling routines, while at the same time the speed of
Fortran code is retained. Besides that, the software allows the results to be exported to popular software packages, like
Generic Mapping Tools, while at the same time utilizing them for quality plots of the results. The modular design of
ISOLA-GUI can be used by users for the addition of supplementary routines in all the stages of processing. An example of
the method’s ability to obtain a quick insight into the complexity of an earthquake is presented, using records from a
moderate size event.
r 2008 Elsevier Ltd. All rights reserved.

Keywords: Seismology; Moment tensor determination; Iterative deconvolution; Extended source; Graphical user interface

1. Introduction studies have to include complete seismograms.

These kinds of studies can provide further con-
Precise knowledge of seismic source parameters is straints on the hypocentral depth and earthquake
a prerequisite for understanding physical processes moment, providing valuable information on a
on faults during earthquakes. Although important region’s deformation style. Therefore, modern
information can be extracted even from the first few seismological networks can provide, in an auto-
seconds of the seismic record (such as earthquake mated way, the moment tensor of the event in
location, its fault-plane solution, etc.), detailed addition to the hypocentre and magnitude. On the
other hand, a variety of methods to invert wave-
$
Code available from server at http://www.iamg.org/ forms into spatio-temporal evolution of the rupture
CGEditor/index.htm or http://seismo.geology.upatras.gr/isola/. process and resulting slip distribution have been
Corresponding author. Tel.: +30 2610 997556;
developed. While the growing amount of global,
fax: +30 2610 990639.
E-mail addresses: esokos@upatras.gr (E.N. Sokos), regional and local data, and increasing facility to
jz@karel.troja.mff.cuni.cz (J. Zahradnik). retrieve them in near-real time through the Internet,

0098-3004/$ - see front matter r 2008 Elsevier Ltd. All rights reserved.
doi:10.1016/j.cageo.2007.07.005
 ARTICLE IN PRESS
968 E.N. Sokos, J. Zahradnik / Computers & Geosciences 34 (2008) 967–977

are challenging enough to steadily update the the Kikuchi and Kanamori method to regional and
modelling methods. local distances, which efficiently combines the
Several codes are available to invert the tele- Fortran and Matlab skills. Fortran for the power
seismic data into point- or multiple-point source and speed for compute-intensive wave propagation
moment-tensor models, including the slip distribu- modelling and source inversion (Green’s functions)
tion models (Dreger and Helmberger, 1993; Kikuchi and Matlab for quick and user-friendly data
and Kanamori, 1991; McCaffrey et al., 1991). handling, control of the inversion process and
Regional data can be inverted by, for example, the plotting of the results. The name of the package is
codes of Dreger and Helmberger (1993), Ammon ISOLA, since the retrieved point-source moment
(2001).1 Complexity of regional and local studies, tensors (called subevents) represent isolated aspe-
compared to the teleseismic ones, is obvious: at rities, or slip patches of the fault. The package has
distances of the order of 100–103 km seismic signals been developed since 2003 and although some
are not solely composed of a few simple phases, high preliminary versions of the Fortran part were
frequencies are abundant, interference waves (such released on the Internet (Zahradnik et al., 2005),2
as, for example, Lg) may predominate on records, this is the first presentation in a journal. The
and near-field effects can be significant. Ray package depends on Matlab for the GUI display,
methods, applicable at teleseismic distances, have data pre-processing and plotting. The Fortran codes
to be replaced by full-wave methods, representing exist as standalone executable files, which are called
complete wavefields. Several frequency–wavenum- by the ISOLA-GUI, in order to perform tasks like
ber methods can be used, such as, for example, the Green’s function calculation and inversion. Users
discrete wavenumber method (Bouchon, 1981). As a may decide to use just the Fortran codes but the use
rule, these demand far more computer power than of the Matlab-based ISOLA-GUI simplifies and
the ray methods. The other reason for computer accelerates the whole procedure, substantially. Be-
speed is due to the fact that the waveform inversion sides Matlab graphics, GMT (Wessel and Smith,
has to include grid search. Indeed, although seismic 1991) is used in all the stages of the procedure, in
moment tensor is related to records via linear order to produce publication-quality figures of the
relations, optimization of the source position and retrieved moment tensors. The M_map software
timing (not precisely known from standard location) package3 is also used mainly for plotting maps
makes the inversion non-linear. All these factors call within the Matlab environment.
for the use of a fast computer language, e.g. Fortran. ISOLA-GUI offers access to all the inversion
However, other factors call for tools with extended procedures through modules that perform specific
processing and plotting capabilities, like Matlab. It is tasks, e.g. data import, filter, inversion, plotting, etc.
because source retrieval is a complex process, consist- Modules can run in parallel and most of them
ing of a series of related stages, requiring an interactive ‘‘remember’’ user’s response or try to suggest
approach and great deal of graphical representation. suitable values of the parameters, in order to speed
Typically, to avoid the need for detailed (unknown) up the processing. Although most of the input
structural models, relatively low frequencies should be parameters are self-explanatory, ‘‘tooltips’’ have
extracted from the data, but, long-period signals may been added for some of them. However these are
fall below the (natural or instrumental) noise level. not a substitute of the program’s manual, which
Thus selection of usable records and their careful contains an extensive description of the method
processing is crucial, and, together with the selection with examples. The user is urged to study the
of a suitable frequency band, requires extensive manual prior to the application of the inversion.
interpretational experience. Matching records with The objective of this paper is to inform the
synthetics, while simultaneously solving problems of research community about the package and its
non-uniqueness of the inversion (e.g. recognizing poor broad variety of options. Specific details of the
resolution of some of the parameters), is another software use are presented in the manual accom-
reason, for an interactive graphical mode. panying the program package. Earthquake physics
Therefore, in this paper, we present a new in this paper is reduced to a necessary minimum and
software package, that is based on an extension of the reader is referred to the paper describing

1 2
Web address: http://eqseis.geosc.psu.edu/cammon/HTML/ Web address: http://seis30.karlov.mff.cuni.cz/.
3
MTinvDocs/mtinv01.html. Web address: http://www.eos.ubc.ca/rich/map.html.
 ARTICLE IN PRESS
E.N. Sokos, J. Zahradnik / Computers & Geosciences 34 (2008) 967–977 969

previous applications of the code to the 2003 component) an option for format conversion from
Lefkada, Greece, earthquake of magnitude Mw6.3 SAC5 and Guralp Compressed Format to ISOLA
(Zahradnik et al., 2005). The code has proved useful native format, is provided. The conversion is
also for other M5+ earthquakes at distances up to accomplished using Matlab codes (m-files) freely
400 km (Adamova et al., 2006; Tselentis et al., available on the Web6 (Thorne, 2005).7 These m-
2006), but also for weaker M3+ and M4 events at files and the accompanying GUI codes can be used
local distances (e.g. see the moment tensor solutions as a guide for any user who wishes to implement his
routinely provided by the University of Patras specific data format conversion. Alternatively, the
Seismological Laboratory to the European Medi- data format conversion can be done using an
terranean Seismological Centre (EMSC)4). external program. Besides the format conversion,
a few options are provided for data quality
2. Program description assessment prior to moment tensor inversion. The
user can test various filters, perform instrument
ISOLA-GUI presented here was created using the correction, integrate or shift the data, etc. These
guide tool of Matlab and its purpose is to allow easy tasks are of major importance since in many cases
interaction with the Fortran codes (mainly input file the data contain various spurious signals, e.g.
preparation and output file display) as well as fast trends, that can affect the inversion, while at the
data pre-processing, Green’s function calculations same time these are not easily recognizable in the
and graphic output of the results. The software is time series data (Zahradnik and Plesinger, 2005). It
rich in graphics and the various plots are created is therefore recommended that these tests are
using a combination of Matlab’s graphic capabil- applied prior to any inversion attempt.
ities as well as GMT software (Wessel and Smith,
1991; Wessel and Smith, 1998). A lot of care has 2.2. Crustal model
been taken to produce a software package that will
be easy to handle, will help the user judge the The crustal model to be used in all the subsequent
accuracy of the inversion results, and will produce steps of the analysis (like Green’s function calcula-
fast and comprehensive output, accurate results and tion, polarity checking, etc) is defined just once
publication ready plots. The program has been during the ISOLA-GUI execution, using the crust-
created with Matlab 6.5 under Windows operating mod tool. The P and S velocities (Vp, Vs), density,
system but has been tested with newer Matlab and P and S quality factors (Qp, Qs) are typed in
versions (Matlab 7.0) as well as with other operating text boxes or read from text files, while the Matlab
systems (Linux) and can work after minor mod- code takes care of the proper formatting, that is
ifications, mainly in the system call commands and important for the ISOLA Fortran codes. A few
in the text file structure. facilities are available for creation of plots of Vp, Vs
The program starts by running the ISOLA m-file, velocities with depth and calculation of Vs or
which generates the main user interface (Fig. 1). density values using empirical relations. The code
From this main GUI, the following functions are as it stands supports a single crustal model for all
available (a) data format conversion, (b) data the stations but future versions will support multiple
inspection, filtering, shifting, (c) crustal model crustal models.
definition, (d) data pre-processing and source
definition, and (e) Green’s function calculation, 2.3. Data pre-processing—source definition
inversion, and plotting of results. These functions
are described in detail in the following sections. The data pre-processing part (instrument correc-
tion, alignment, integration, decimation, etc.) must
2.1. Data (format conversion—inspection) (utilities be carried out carefully, since any error introduced
I, II) at this stage can cause significant errors in the
inversion step. Using ISOLA-GUI, the user can
Although ISOLA’s native format is quite simple apply all the above corrections on the data and at
(4 column ascii files, i.e., time, NS, EW, Z
5
Web address: http://www.llnl.gov/sac/.
4 6
Web address: http://www.emsc-csem.org/index.php?page= Web address: http://www.guralp.net/.
7
current&sub=qmt. Web address: http://gcc.asu.edu/mthorne/saclab.
 ARTICLE IN PRESS
970 E.N. Sokos, J. Zahradnik / Computers & Geosciences 34 (2008) 967–977

Fig. 1. Main ISOLA Matlab GUI.

Fig. 2. ISOLA-GUI tool for data pre-processing.

the same time have a visual control of the results. under the epicentre, while with the extended source
The user is guided by the program in performing all option, sources on a horizontal line or on a plane
the tasks in the proper order, while the GUI keeps with arbitrary orientation can be defined. The
track of user’s last response and stores files in ‘‘point’’ source option is used for a first approxima-
proper folders. In this way the process is fast and tion of the inversion problem and depends on the
well controlled (Fig. 2). If there is a need for some quality of the given epicentre. Using this type of
special correction or data processing, a new module inversion can indicate the inversion’s preferred
can be easily added in the procedure with minimum depth for the event, but if the epicentre is wrong
Matlab programming. the results could be misleading. Thus in many cases
Definition of trial source positions for the event it may be desirable to rerun the inversion around the
source is done through a GUI code, named epicentre, using trial positions on a plane, or on a
sourcepre, which has been specially designed to few horizontal planes at different depths. In this
automate and simplify this task. Two options (the way the inversion will find the optimum position of
‘‘point’’ or ‘‘extended’’ source) are available to the the event in space. The option of an arbitrary trial
user. By point we mean trial sources that are located plane may be useful when studying the fine structure
 ARTICLE IN PRESS
E.N. Sokos, J. Zahradnik / Computers & Geosciences 34 (2008) 967–977 971

of the source in the nodal planes. The ISOLA-GUI the number of subevents, and the time span of their
creates all the necessary files for the Green’s occurrence. Various helpful features are included,
function calculations—inversion, along with plots like a special option for automatic calculation of
of the source distribution (in geographical and station weights, based on the maximum amplitude
Cartesian coordinates) for user inspection. Under of the displacement record, and a tool for plotting
the Matlab environment the plots are produced correlation and focal mechanism vs source position
using the M_Map tool, but GMT format files are and its time shift (Fig. 4). Besides the correlation,
created too, thus the user can produce customized the double couple percentage of the solution (DC%)
maps latter. can also be plotted, since in many cases it can
provide indication about problematic inversion
results. The correlation plot is quite important and
2.4. Green’s function computation—inversion a major advantage of the software is that the user
can have visual control of the inversion results,
The Green’s function calculation is performed while the inversion progresses, and actually drive
using the frequency wavenumber method (Bouchon, the inversion procedure based on this diagram and
1981). The code is written in Fortran but the user relevant constraints, e.g. polarity measurements,
controls its use through the ISOLA-GUI, using the rupture velocity, the assumed tectonic style in the
greenpre tool. We have tried to limit the number of area, etc. The application of constraints may be
input parameters as much as possible, thus the user desirable depending on the resolution of the
has to specify only the maximum frequency of the inversion parameters. For example, for the case of
Green’s function computation. Then the code low resolution the constrained (fixed) focal mechan-
automatically generates the Green’s functions and ism may prevent ‘‘wild jumping’’ of subevents in
the corresponding elementary seismogram files, space and time. In this sense, a requirement for user
while at the same time it stores everything in the to select his preferred choice of the subevent
proper folders. position and time is quite common and definitely
The inversion is the core of the procedure and is improves the results compared to the fully auto-
performed using another ISOLA-GUI tool called matic run, which is based only on the absolute value
invert, presented in Fig. 3. The user can select the of the correlation. In any case both user’s and
type of inversion (full, deviatoric, double couple automatic selection of subevents are possible during
(DC), fixed), the frequency band of the inversion, the inversion procedure.

Fig. 3. ISOLA GUI used for inversion calculations.

 ARTICLE IN PRESS
972 E.N. Sokos, J. Zahradnik / Computers & Geosciences 34 (2008) 967–977

11
100
90
10
80
70
9 60

DC%
50
8 40
0.7 0.6 0.5 0.4 0.4 0.5 0.6 0.7 0.8 0.8 0.7 0.6 0.5 0.6 0.7 30
7 20
Source position

10
6 0

5 0.4 1.0
0.9
4 0.8
0.8 0.5 0.7
0.9

Correlation
3 0.6
0.5
2 0.4
0.3
0.3 0.2
1
0.1
0.0
0
-8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12
Time (sec)

Fig. 4. Plot of correlation vs time shift, source position and focal mechanism for single source inversion. Largest correlation was obtained
for source 2 (6 km depth) and 3.16 s time shift. Focal mechanism shading changes according to double couple percentage, while for
diagram clarity, we only plot solutions having double couple percentage larger that 70%. Preferred solution is depicted by a larger
beachball.

2.5. Plotting of results 3. Case study, Zakynthos 2006 seismic sequence

The last step is the plotting of the inversion To demonstrate the use of the software, we
results, in particular the moment tensor solution present the moment tensor analysis for a strong
and the best fit between data and synthetics (Fig. 5). earthquake that occurred on April 11, 2006 at
Again, a special tool was designed, named plotres, 17:29 UTC, offshore of Zakynthos island in western
which uses Matlab and GMT graphic capabilities to Greece, (EMSC location Lat:37.71N, Lon:20.91E,
prepare quality plots of the results. The variance Mw ¼ 5.6). The event was recorded by three broad
reduction per tested source and the corresponding band instruments of PSLNET, a new satellite
focal mechanism can be plotted together on a graph telemetry network, starting operation in 2006 and
in order to inspect the inversion convergence belonging to the Seismological Laboratory of the
through the whole computation. An option for University of Patras (Tselentis et al., 2006). This
comparing the computed focal mechanism solution event was the strongest from a series of medium size
with the available first motion polarities is also earthquakes (more than five events of Mw44.5)
available. The plotting tool was designed primarily that affected the island of Zakynthos for a period of
for users with limited knowledge of usage of GMT 2 months causing moderate damage to the island’s
or Matlab graphics. Nevertheless our intention was structures.8 The records of this event, as SAC files,
not to prepare a GMT front-end but help users are included with the program distribution to serve
produce nice maps quickly. Experienced users can as an example. The same event will be used here to
easily produce custom plots using the output text demonstrate the analysis using ISOLA-GUI.
files with the inversion results. Additionally, Matlab Initially, the data are converted to ISOLA native
programmers can easily add their own special format using the specially designed SACimport tool.
features to the plotres tool, such as creation of a Then the quality of the records and the signal to
text file with the moment tensor solution, suitable
8
for email notification. Web address: http://www.itsak.gr/Zakynthos_2006SMR.pdf.
 ARTICLE IN PRESS
E.N. Sokos, J. Zahradnik / Computers & Geosciences 34 (2008) 967–977 973

Fig. 5. (a) Waveform comparison for single source inversion. (b) Moment tensor solution results, output by ISOLA-GUI for single source
inversion.
 ARTICLE IN PRESS
974 E.N. Sokos, J. Zahradnik / Computers & Geosciences 34 (2008) 967–977

noise ratio at various frequency bands are checked 3.2. Multiple source inversion
using the datafilt tool in order to define the proper
frequency band for the inversion. The next step is to Multiple source inversion is another key point of
prepare the crustal model that will be used in the the ISOLA method. Usually this type of inversion is
analysis, in this case it is the model previously performed for large earthquakes using data at
proposed by Haslinger et al. (1999) during local teleseismic distances, but with ISOLA it is possible
earthquake tomography studies in western Greece. to perform the same type of inversion using local or
The model is converted to the proper format using regional data. In this way more details of the seismic
the crustmod tool. Next steps involve input of event source can be revealed.
information, selection of stations to be used in the In order to run a multiple source inversion, a set
inversion, and raw data pre-processing, which are of possible source positions has to be defined; this is
done using the eventinfo, stationselect, and datain accomplished with the sourcepre tool. In Fig. 8, a
tools (Fig. 1). At this point, we can decide on single map of the trial source positions for the studied
or multiple point source inversion and prepare the event is displayed. In this multiple source inversion
necessary trial source position and Green’s function run, we use source positions located on a horizontal
files using the sourcepre and greenpre tools, respec- plane at a depth of 6 km (the depth preferred in the
tively. We will present both single and multiple single source inversion) having a spacing of 12 km;
source inversion for this event, starting with the first we use a grid of 25 trial source positions. Using this
one. set of sources and the appropriate Green’s func-
tions, we start running the inversion and subse-
3.1. Single source inversion quently searching for two elementary sources, using
the invert tool. The selected frequency band for the
We select deviatoric point source inversion for a inversion is extended to higher frequencies
series of trial source positions lying at various (0.02–0.13 Hz, with cosine tapering) in order to
depths below the epicentre position determined by obtain detailed information about the source
EMSC. The total trial sources tested were 10, with a rupture process. During the inversion ISOLA
3 km vertical separation, spanning depths from 3 to pauses and asks the user to accept the automatically
30 km. This type of inversion will indicate the detected optimum source position and time, or to
preferred centroid depth, and it is a logical first step propose a different solution for the first subevent.
during the moment tensor study. Using the invert Creating a plot of all the possible solutions assists
tool, inversion parameters such as the frequency the user to decide if the proposed automatic
band for the inversion, and parameters of the time solution is acceptable. We need to note here that
search can be set. For this single source inversion ISOLA automatically proposes a solution based on
run, we used frequencies extending from 0.037 to the maximum correlation. This solution can differ
0.062 Hz with cosine tapering applied to both ends greatly in the DC% when compared with another
(Fig. 5a). Using the same tool, during the inversion ‘‘nearby’’ solution with a slightly smaller correlation
run, a plot of correlation and focal mechanism vs value. Having a tool that displays this in an efficient
source position and its time shift is constructed. graphical way allows the user to explore all the
Such a plot is presented in Fig. 4 where it is concurrent solutions and choose one, based on
clear, from correlation maximum, that the preferred various assumptions. Such a plot for the first
solution is at a depth of 6 km (source number 2 subevent is presented in Fig. 6. The automatic
corresponds to this depth) for a time shift of solution is marked with a larger beachball, while an
3.16 s. The results of the inversion are plotted later option for plotting results for a single source is also
using the specially designed plotres tool (Fig. 5a, b). available (Fig. 7).
What is interesting from the single source results is After a detailed check and a combined use of the
that the DC% is really low for this event, indicating above diagrams, one can easily identify other
that the possible effect of the source complexity possible solutions that have slightly smaller correla-
can be ‘mapped’ into a large but only apparent tion and thus have not been automatically proposed
CLVD value. Of course this is just one explanation by the program. For example, source 8 at time shift
for the large CLVD value but it motivates us to 0.6 s, or source 13 at time shift 4.6 s, has similar
display the use of ISOLA, for multiple source correlation and focal mechanism, but differ sig-
inversion. nificantly from the automatically chosen solution in
 ARTICLE IN PRESS
E.N. Sokos, J. Zahradnik / Computers & Geosciences 34 (2008) 967–977 975

26
25 100
24 0.6 90
23 80
22 0.4 70
0.2 60
21 0.2

DC%
0.4
20 0.4
0.6 0.4
50
19 40
18 0.6 30
17 0.4 20
0.4
16 0.2 0.2
10
Source position

15 0.6 0.4
0.4
0
14 0.4
0.4
13 0.8 0.6
0.6
12 0.4 0.8
11 0.2
10 0.4 0.4
9 0.4 0.6
0.6 0.6
8 0.6

Correlation
7 0.4 0.4
0.4
6 0.2 0.4 0.4
0.4
5 0.4
4 0. 2
3 0.6 0.6
2 0.4
1 0.4
0.0
0
-8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12
Time (sec)

Fig. 6. Correlation diagram for first sub-event during multiple source inversion. Horizontal axis is time shift, vertical axis is trial source
number, and contours are correlation values. Focal mechanism shading changes according to double couple percentage. Only solutions
having double couple percentage larger that 70% are plotted. Preferred solution is depicted by a larger beachball.

0.9
Source number 13

0.8

0.7
Correlation

0.6
100
0.5 90
80
70
0.4 60
DC%

50
40
0.3 30
20
10
0
-8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12
Time shift (sec)

Fig. 7. Correlation plot for the first sub-event. Horizontal axis is time shift and vertical axis shows correlation; automatic solution
proposed by the program is depicted by an arrow.

their source position and DC%. At this point, based solution instead of the automatically proposed 4 s.
on what has been described already, we select source The focal mechanism is the same for both solutions,
position 13 and time shift of 4.6 s as our preferred but the DC% increases from 60% to 95%, with just
 ARTICLE IN PRESS
976 E.N. Sokos, J. Zahradnik / Computers & Geosciences 34 (2008) 967–977

a 3% difference in correlation. Such a choice is ISOLA-GUI presented in this paper allow an easy
perhaps better compatible with the tectonic regime application of the iterative deconvolution method of
in Zakynthos, since it explains the first sub-event as Kikuchi and Kanamori (1991), for local and
a pure shear source. Following the same procedure regional events, thus the method can be used for
for the second sub-event, we end up with an both single and multiple source geometry. The
explanation of the waveforms as being due to Matlab-based GUI facilitates easy data handling,
rupture of two pure-shear sources (Fig. 8). The while at the same time providing complete user
focal mechanism of the second subevent is different control during all the processing steps and easy
from the first one, suggesting a complicated rupture, graphical display of the results, using both Matlab
and possibly explaining the large CLVD part for the and GMT capabilities. The waveform inversion is
single source solution. Nevertheless, final conclu- done using the Fortran ISOLA code, in order to
sions concerning the rupture history of this event take advantage of language speed, while user
would need special care (Zahradnik et al., 2007), interaction even during the inversion is within the
and a full physical justification of the DC% would Matlab environment. Various tools are available for
require more efforts, well beyond the scope of the the user to explore the data quality, adjust the
present paper. computational parameters, correct the data, and
display inversion results in a publication-ready
4. Conclusions form. The modular design of the GUI allows the
easy upgrade as well as the inclusion of user-created
The moment tensor analysis of local and regional additional modules, using elementary skills in the
events is a standard procedure applied to broad- Matlab programming language. The software, the
band seismic networks worldwide, in many cases manual, and examples files are available for down-
even as an automatic procedure. In all such load at http://seismo.geology.upatras.gr/isola/, for
applications the inversion is performed for a the Windows operating system, while a Unix/Linux
single-point source. The programs ISOLA and version will be available soon.

Fig. 8. Plot of 25 trial source positions, along with final results for multiple source inversion.
 ARTICLE IN PRESS
E.N. Sokos, J. Zahradnik / Computers & Geosciences 34 (2008) 967–977 977

Acknowledgments Kikuchi, M., Kanamori, H., 1991. Inversion of complex body

waves; III. Bulletin of Seismological Society of America 81,
The authors would like to thank Petra Adamova 2335–2350.
McCaffrey, R., Abers, G., Zwick, P., 1991. Inversion of
for extensive testing of the software. Some parts of teleseismic body waves. In: Lee, W.H.K. (Ed.), Digital
the Matlab GUI use functions freely available on Seismogram Analysis and Waveform Inversion. International
the Web, like the rsac and coral distribution. In the Association of Seismology and Physics of the Earth’s Interior,
Fortran section, several subroutines written by El Cerrito, CA., pp. 81–166.
other authors were also used: M. Bouchon, O. Thorne, M.S., 2005. Broadband waveform modelling of deep
mantle structure. Ph.D. Dissertation, Arizona State Univer-
Coutant, J. Sileny, J. Jansky and O. Novotny, J.A. sity, United States, 283pp /http://www.aeic.alaska.edu/
Snoke, P. Reasenberg and D. Oppenheimer. Codes input/mthorne/pubs/index.htmlS.
from Fortran Numerical Recipes were also em- Tselentis, G.-A., Sokos, E., Serpetsidaki, A., Zahradnik, J.,
ployed. This work was supported by the EC project Jansky, J., Plicka, V., 2006. New satellite BB network of
3HAZ-Corinth and the following research projects western Greece and its application to the M5 earthquake
sequence of Zakynthos, April 2006. In: Abstract Book First
in Czech Republic: GAUK 279/2006/BGEO/MFF, European Conference on Earthquake Engineering and
GACR 205/07/0502, MSM0021620860. Seismology, Geneva, Switzerland, p. 273 /http://www.ecees.
org/S.
Wessel, P., Smith, W.H.F., 1991. Free software helps map and
References display data, EOS. Transactions of the American Geophysical
Union 72 (41), 445–446.
Adamova, P., Zahradnik, J., Stavrakakis, G., 2006. Moment Wessel, P., Smith, W.H.F., 1998. New improved version of
tensor inversion of possibly multiple events at regional generic mapping tools, EOS. Transactions of the American
distances. In: Abstract Book First European Conference on Geophysical Union 79 (47), 579.
Earthquake Engineering and Seismology, Geneva, Switzer- Zahradnik, J., Plesinger, A., 2005. Long-period pulses in broad-
land, p. 450 /http://www.ecees.org/S. band records of near earthquakes. Bulletin of the Seismolo-
Bouchon, M., 1981. A simple method to calculate Green’s gical Society of America 95 (5), 1928–1939.
functions for elastic layered media. Bulletin of Seismological Zahradnik, J., Serpetsidaki, A., Sokos, E., Tselentis, G.A., 2005.
Society of America 71, 959–971. Iterative deconvolution of regional waveforms and a double-
Dreger, D.S., Helmberger, D.V., 1993. Determination of source event interpretation of the 2003 Lefkada earthquake, Greece.
parameters at regional distances with single station or sparse Bulletin of the Seismological Society of America 95 (1),
network data. Journal of Geophysical Research 98, 8107–8125. 159–172.
Haslinger, F., Kissling, E., Ansorge, J., Hatzfeld, D., Papadimi- Zahradnik, J., Sokos, E., Tselentis, G.-A., Martakis, N., 2007.
triou, E., Karakostas, V., Makropoulos, K., Kahle, H.G., Non-double-couple mechanism of moderate earthquakes near
Peter, Y., 1999. 3D crustal structure from local earthquake Zakynthos, Greece, April 2006; explanation in terms of
tomography around the gulf of Arta (Ionian region, NW complexity. Geophysical Prospecting 55, 1–16.
Greece). Tectonophysics 304 (3), 201–218.