2. VARIOUS TOPICS IN
SEISMOLOGY – TECTONICS PERTAINING TO EQ PREDICTION
There is a general notion, in the seismological community, that
earthquakes can occur in a random way in every place in a large area
that exhibits seismicity. This is corroborated by the existing
theories with regard to the generation mechanisms of earthquakes and
the strain charge conditions held in such places. The Self Organized
Criticality (SOC) theory, proposed by Bak (1966); Bak; Tag; and
Wiesenfeld (1988), indicates that physical systems, like the
seismogenic areas, are at critical state or instability and
therefore, any small earthquake has some probability to cascade into
a large earthquake. Other researchers (Otsuka 1972a, b; Vere-Jones
1976; Bak and Tang 1989; Ito and Matzutaki 1990) concluded that the
generation of a strong earthquake depends on very small variations
of stored elastic energy, fault strength or variations of the
elastic properties of the seismogenic region, therefore,
statistically, a strong earthquake can take place anywhere.
Since all the latter theories are combined to the fact that the
Earth, as it is deduced from the geological, tectonic maps, compiled
by the geologists, is highly fractured, then it is justified to
conclude that strong earthquakes can, in a random way, occur almost
anywhere. However, as the following analysis will demonstrate, the
localization of strong earthquakes follows a scheme which is
prescribed by the deep tectonics and the intense fracturing of the
lithosphere, as far as it concerns the shallow earthquakes which
occur in the crust.
2.1 Spatial distribution of strong EQs.
Before any attempt is to be made for any earthquake prediction, it
is required to be known at least the general tectonic setting for
the area under study. This is justified by the fact that EQs are the
results of intense, tectonic processes which take place in the
lithosphere. In this work a different approach, in relation to EQs,
was adopted. Applied geophysical (gravitational) methods have been
used in an attempt to explain the spatial distribution of strong EQs
and to prove that these occur only on predefined narrow faulting -
fracture zones of the lithosphere, detectable in a deterministic
scheme. The theory of lithospheric plates and their motion is well
known and has been adopted by the seismological community. A
generalized picture of the Earth's lithospheric plates is presented
in figure (2.1.1).
Fig. 2.1.1. The Earth's lithospheric plates (after USGS).
The main lithospheric Earth's plates are subdivided in smaller
ones. McKenzie (1972, 1978) presented a small-scale lithospheric
model for the particular case of the Greek territory. This is
presented in the following figure (2.1.2).
Fig. 2.1.2. McKenzie's (1972, 1978) lithospheric, plate model,
proposed for the Greek territory.
As these plates are in continuous motion, due to applied forces
from the surrounding larger plates, extension or compression is
applied on them as a result. Therefore, extensional or compressional
fracture zones are formed. The kinematics of the Greek territory
tectonic plates was studied by Papazachos et al. (1996) and is
presented in the following figure (2.1.3). It is understood that
strong seismic events are expected to occur where large tectonic
activity takes place.
Kinematics of the Greek territory tectonic plates proposed, by
Papazachos et al. (1996).
The main features of this map are the areas which are characterized
by compressional, extensional and strike-slip motion. As it will be
explained later on, this plays an important role in the kinematics
of the Aegean area and moreover facilitates the calculation of its
angular, rotational velocity.
2.1.1 Mapping of major, seismic fracture zones –
Mapping of fracture zones - faults is of great importance in
seismological studies. These are the places where earlier EQs took
place and most probably future EQs will occur. Therefore, the
accurate mapping of these fracture zones - faults is of highest
priority. The problem that arises is: what kind of faults - fracture
zones must be mapped and is it always possible?
Earlier studies have revealed rather broad, seismic zones of
increased, seismic risk, based, on the Greek seismicity history. In
next figure (22.214.171.124) De Bremaecker et al. (1982) divided the
regional, Greek area into smaller blocks which generally move
towards southwest with different velocities in respect to the plate
Fig. 126.96.36.199. Finite element network with boundary condition
velocities, regarding European plate as proposed by De Bremaecker et al. (1982) for the Greek
In another study, Papazachos (1988) defined the most active seismic
areas in Greece, based on the spatial distribution of the shallow
earthquakes. The latter is shown in next figure (188.8.131.52).
Seismic zoning of shallow earthquakes in the Greek and surrounding areas (Papazachos, 1988).
Furthermore, the Greek territory was divided in specific zones
according to the seismicity as it is expressed by the yearly
released, seismic moment per 10000 Km2 (Papazachos, 1989). This is
presented in the figure (184.108.40.206).
Fig. 220.127.116.11. Seismicity of the
seismic zones of the Greek territory as it is expressed by the yearly, released, seismic moment per 10000 Km2 (Papazachos,
In a more recent study, Hanus and Vanek (1993) presented another
mapping of the active, seismic zones in the Greek territory which is
presented in the following figure (18.104.22.168).
Fig. 22.214.171.124 Mapping of the
active, seismic zones in the Greek territory as presented by Hanus and Vanek, (1993).
Comparing the latter three figures (126.96.36.199, 3, 4) it is shown that
the mapping of the active, seismic zones in the Greek territory
(concerning shallow seismicity) is very subjective to each
researcher. The characterization of an area, as an active, seismic
one, is based on subjective criteria of each scientist, who presents
his work on this topic. The only common, broad feature between them
is the absence of active, seismic zones in the central Aegean area (Cyclades).
An objective, active, seismic zoning of the Greek territory can be
obtained by taking into account the deep lithospheric fracture zones
As far as it concerns the detailed mapping of seismic fracture
zones and faults, classical geological surface observations cannot
detect all of them. This is especially true when no surface trace
has been produced by any deep, subsurface, tectonic activity.
Furthermore, it is difficult to assign a seismological significance,
to a detected on surface fault, without any other deep tectonic
knowledge. The previous problem has been overcome, partly, by
studying the delineation of the seismic events over a certain area.
An example of this procedure is presented in the following figure
of a deep fracture zone – fault from its seismic activity,
mapped over a certain period of time (June,
2001). The seismic activity for the assigned, spatial window
(dashed square) was mapped in June 2001. It is obvious the
delineation of the seismically activated fault. As a result
of this activity a strong seismic event occurred at the
lower left part of this fault (Psara EQ, Ms = 5.6R, June
2.1.2. Mapping of major, seismic fracture zones -
faults from the study
of the Earth's
The creation of fracture zones and faulting of the Earth is due to
the present stress field conditions which are applied in a
particular area. Fracturing takes place in various modes (Mattauer,
1973), when extensional or compressional forces are applied on it.
In the following figure (188.8.131.52) the black arrows indicate the
applied forces, while the lines in the solid block indicate the
various generated modes of fracturing. In case (A) there is block
movement which creates a typical thrust fault combined with internal
micro-fracturing typically normal and parallel to the direction of
the applied stress field. In case (B) the generated, internal,
micro-fracturing is interconnected and delineated diagonally in
respect to the applied stress field direction, while in case (C), a
characteristic strike-slip fault has been formed, due to exceptional
type (B) stress field conditions.
184.108.40.206. Fracturing –
faulting which is generated by applied stress (Mattauer,
1973). The arrows indicate the applied stress direction (compressional,
The form of a typical fault is presented in the following figure
(220.127.116.11), regardless of its cause of origin. Typically, the ground
is split into two parts through the fault plane, which exhibit a
differential mode of relative movement to each other. The
slip-vector of the generated fault, characterizes the latter.
sketch of a down-throw fault. The thick vector represents
the slip vector; the thin vectors indicate the relative
movement of each ground block to the other one. The
down-throw block has slipped along the fault plane.
The slip-vector, that characterizes a fault, can be generally
considered as the combination of its orthogonal components, when it
is projected to an orthogonal coordinating system X, Y, Z. This is
presented in the following figure (18.104.22.168).
In each case, the interrelation of the magnitudes between these
components can generate the various types of faulting which are
observed in nature. Moreover, by accepting that the fault plane
coincides with the Z-X plane, then the SV-Y component equals to
zero. In such a case, (when SV-Z = 0), no thrust, neither normal
faulting can exist. Since the slip-vector lays in the fault plane,
then the only possible, ground block movement is along the strike of
the fault. The latter is the case of a strike-slip fault and is
characterized as sinister or dextral one depending on the direction
of the block movement. In case SV-X = 0, then normal faulting or
thrust faulting is generated, depending on the sign of the SV-Z
component. In nature, faulting that exhibits both types of faulting
modes, is, typically, generated. The type of generated faulting is
characterized by the dominant and much larger component of the fault
22.214.171.124. Analysis of a
vector into three orthogonal components. The fault
slip-vector is analyzed into its orthogonal components SV-X,
SV-Y and SV-Z.
In a most general case, the slip-vector is analyzed into orthogonal
components, along the X-Y horizontal plane and the vertical to the
ground Z-axis. In such a case the vertical component VS-Z plays a
very important role in the overall tectonic and the stratigraphic
conditions which characterize an area. A relative, vertical movement
of the two ground blocks of the fault has, as a consequence, the
stratigraphic modification of the overlaying geological formations.
The latter has as an effect, the generation of ground lateral
density discontinuities, and therefore, are generated changes in the
intensity of the Earth’s gravity field. This is the physical basis
upon which the gravity methodologies for the study of the geology
tectonics and stratigraphy of an area were founded.
Large tectonic fault / fracturing zones of the Earth are mainly
reflected on the morphology of its gravity field. The wavelength of
these features depends on the dimensions of the corresponding fault
/ fracture zones. Therefore, by studying the large wavelengths of
the gravity field, we can obtain information for the deep large
faulting / fractures zones of the Earth.
Moreover, by studying the horizontal gradient of the corresponding
gravity field, we can locate the precise location of a fault /
fracture zone quite easily. In the following figure (126.96.36.199), (a)
represents a geological fault of the Earth, (b) is the corresponding
gravity field, while in (c) is presented the corresponding
horizontal gradient of it. Details on this topic can be found in any
geophysical gravity textbook.
188.8.131.52. (a) =
Geological fault model, (b) = Corresponding gravity field,
(c) = Corresponding gravity horizontal gradient.
Generally, the gravity method of Applied Geophysics is capable of
tracing lateral density discontinuities of the underground,
2.1.3 Application of the methodology in Greece.
On 13th May, 1995 a devastating earthquake (M=6.1R, Lat=400.18,
Lon=210.71) occurred in "Grevena" area, northern Greece and caused
many damages at the nearby towns and villages (some deaths were
What is interesting, in this earthquake, is the fact of total
absence of any statistical seismo-logical indication that a strong
EQ could take place, in an area like this, which was characterized
as of very low seismic risk (zone 1, out of 4) and of total absence
of any known, large, geological tectonic features to justify any
possible strong seismic event.
The earthquake which occurred, the absence of statistical and
geological - tectonic data which could reveal the seismically
dangerous character of the "Grevena" area and the required,
deterministic cause of this earthquake motivated the analysis to
The very same scenario was replicated on 7th September 1999 in
Athens, the capital city of Greece. Although Athens is located in an
area which is characterized as "zone 2", as far as it concerns its
seismic risk, the complete absence of knowledge of deep tectonics
had postulated the notion that the capital of Greece is built over
"safe" ground. Unfortunately, for the people (over a hundred) who
died, apart from the very large damages in buildings, during this
earthquake, this was proved to be totally wrong, after this
So the question that arises is: Where do the strong earthquakes
The following analysis may be the answer to the above question.
It is well known that strong earthquakes occur when a fracture zone
/ fault releases its stress load and rupture occurs along a large
part of it. The longest the fault / fracture zone has been
activated, the largest the magnitude of the earthquake is.
The problem therefore, is formulated as follows: Is it possible to
map the deep fracture zones, where strong EQs occur? The study of
strong earthquakes in Greece, statistically treated, may vaguely
reveal the large, deep, tectonic features that are prerequisites for
its occurrence. The same problem, of the detection of deep tectonics
of an area, faced from the view of Applied Geophysics, is rather
simple, as it will be explained later on.
The key point of this study is that strong earthquakes must
coincide in location with large tectonic faulting / fracturing
systems of the lithosphere. Consequently, the specific target of
this analysis is the location of these tectonic features, which
locations can't be traced by other geological methods.
The following map (fig. 184.108.40.206), where is presented the location
of the strong, (M>=6R), seismic events of the Greek territory, which
occurred during the period 1900 – 1997, has been compiled, in order
to have a clear view on this topic. The circle, around each
location, is the area defined by the (conservative) location
estimation error (+/- 25Km). It must be pointed out that the
location of each earthquake can be anywhere in the error circle,
without decreasing the validity of the map.
with Ms>=6R which occurred in Greece during the period 1901
- 1997. The red circle indicates the error, accepted, for
the location determination that equals to +/- 25Km. The map
correlates in size with the Bouguer Anomaly map of Greece (Makris,
The gravity field transformation approach into its gradient can be
applied in two-dimensional form on the gravity field of Greece. A
simplified presentation of the Greek gravity field, Bouguer anomaly
(Makris, 1981) is presented in the following figure (220.127.116.11).
18.104.22.168. Simplified (Thanassoulas,
1998) gravity (Bouguer anomaly) map of Greece (Makris,
1981). The scale is in mgals.
The following scheme has been applied in order to utilize the
horizontal gradient transforma-tion.
Gradx,y(G) = G(x,y)*TL(x,y)
Where: G(x,y) is the gravity field at point x, y
TL(x,y) is the operator, applied, for the transformation to be
Gradx,y(G) is the resulted, horizontal gradient
TL has been calculated analytically by the use of a polynomial
surface of second order, fitted, to a progressively sliding window
over the gravity field data.
In order to avoid near surface tectonic features, and targeting to
the deeper structure of the Earth, a 2-D "window" of 20 x 20Km was
used in all these transformations. The results of this operation are
presented as photo relief map, illuminated, from NE direction. In
the next figure (22.214.171.124) is shown the result of the transformation
of the Bouguer anomaly map of Greece into a gradient one.
gradient gravity map of Greece (Thanassoulas, 1998).
The latter map is presented in 3-D photo relief (fig. 126.96.36.199), so
that the high gradient value areas are better visualized.
Deep fracture zones express (horizontal gravity gradient ridges)
themselves as the bounda-ries between dark black and white zones.
188.8.131.52. 3-D photo
relief compiled gradient map of Greece. Each scale unit is
10Km. Lighting from NE (Thanassoulas, 1988).
It is very interesting to note the deep fracture zone, crossing in
a NW - SE direction the Attiki area (X=32, Y=34), where Athens is
located and the strong EQ of Athens took place. This feature,
unfortunately, reveals the dangerous and seismically risky position,
where the capital city of Greece is located.
The detailed fracture zones-faults (thick black lines) locations,
determined by the latter method, are presented in the following
lithospheric fracture zones (thick black lines) determined
from the horizontal gradient map of Greece (Thanassoulas,
1998). The solid dots indicate the locations of strong
(Ms>6R) EQs, which occurred during the period 1901 – 1997.
2.1.4 Strong EQs location map of Greece (1901- 1997).
It was stated earlier that strong earthquakes have their origin in
deep and large faults / fracture zones of the Earth. In the
following figure (184.108.40.206), is made a comparison between the
location of the strong earthquakes, determined by seismological
methods, and the location of the deep faults / fracture zones
determined by the transformation of the gravity field into its
horizontal gradient. The EQs catalog, which is used, is the one
provided by the National Observatory of Athens (NOA).
Therefore, it is expected that at least a fault zone crosses the
error location circle of each strong earthquake, which is mapped.
This fault zone is, presumably, the actual fault zone that was
activated and generated the corresponding earthquake. Consequently,
the "most probable true" location, of the corresponding earthquake,
is not the one suggested by seismological methods, but the one which
resulted by the projection of the "seismological" location of the
epicenter, to the nearest fault / fracture zone.
220.127.116.11. Location of
strong, seismic events (Ms>=6R, solid red dots) in Greece,
for the period 1901 – 1997, in relation to the location of
the fracture zones - faults (thick black lines), determined
by the transformation of the corresponding gravity field.
A detailed examination of this map indicates that the vast majority
of the earthquakes are located in a distance shorter than 25Km from
the nearest fracture zone. Therefore, if the location of each
earthquake is moved within its location error radius towards the
collocated fracture zone, then the validity of the EQ location map
does not change, while at the same time there is a very good
coincidence of the EQ locations with the location of the fracture
zones (fig. 18.104.22.168). A very small number of EQs, mainly in the
western part of Greece, do not follow this rule. This could be
attributed either to missing existent fracture zones / faults which
the transformation of the gravity field didn't succeed to identify
or to the wrongly calculated, by the seismological methods, location
of the deviating earthquakes. In either case the degree of
correlation (statistically) is very high.
location, of the presented EQs, according to the calculated
fracture zones - faults from gravity map transformation. The
EQs data file contains data from 1901 to 1997, when the
original study was presented (Thanassoulas, 1997). The EQs,
outside the gravity map boundaries, have been excluded.
A detailed example, showing the validity of the methodology, is
presented in the following figure (22.214.171.124).
of a deep fracture zone - fault from its mapped, seismic
activity over a certain period, compared with the one,
calculated, by the gravity field transformation. A (left) =
fracture zones, calculated, through gravity field
transformation, B (right) = fracture zone, mapped, from its
seismicity, which took place during June 2001. Comparison is
made between the inset frames.
It is very interesting to see how the seismically, mapped, fracture
zone in fig. (126.96.36.199 - B) coincides with the one which is mapped by
the gravity method and is presented in figure (188.8.131.52 - A). Further
more the specific fracture zone extends towards NE in the area of
Western Turkey. This fracture zone has already produced two strong
EQs with Ms >6R in the past.
2.1.5. Verification of fault zones by strong EQs
which occurred during the period 1998-2006.
Accepting the validity of the presented methodology, it is applied
as a test on the earthquakes which took place during the period
1998-2006 in Greece. These earthquakes must coincide with the
corresponding, deep, lithospheric fracture zones / faults.
In the following figure (184.108.40.206), are presented for validation
purposes the strong (Ms>=6R) earthquakes (solid red circles) which
took place in Greece, after the completion of this study (Thanassoulas,
1998), that is the period 1998-2006, along with the determined fault
/ fracture zones.
of the location of recent strong (1998 - 2006, Ms>=6R)
seismic events (red circles) to the location of the already
calculated fracture zones - faults by the transformation of
the Greek gravity field.
It has been made clear that all the latest, strong, seismic events
coincide with the location of the fracture zones - faults which have
been already calculated.
Two more examples, related to Lefkada EQ (2003/08/14, Ms = 6.4R)
and Kythira EQ (2006/01/08, Ms = 6.9R), follow. Lefkada EQ caused
large damages on Lefkada Island and is the most destructive recent
one. The above mentioned EQ is presented in the following figure
220.127.116.11. Lefkada EQ,
red circle (2003/ 08/14, Ms = 6.4R) in relation to the deep,
lithospheric fracturing (brown lines), calculated by the
gravity field transformation (Thanassoulas, 1998).
Kythira EQ is the strongest one which occurred in Greece, recently.
It was felt almost allover Greece and the neighbor countries, too.
Kythira EQ is presented in figure (18.104.22.168).
22.214.171.124. Kythira EQ,
red circle (2006/01/08, Ms = 6.9R) in relation to the deep,
lithospheric fracturing (brown lines), calculated by the
gravity field transformation (Thanassoulas, 1998).
The presented analysis of strong earthquakes’ location which
earthquakes occurred in Greece during the last 100 years (more or
less) and the detailed examples, too, show clearly that, contrary to
seismological notions that strong EQs can occur, statistically, in
any place, these happened in specific, deep, lithospheric fracture
zones / faults. These can be mapped by the use of applied
geophysical methods, and in particular, through the analysis of the
Earth’s gravity field that is modified by these strong, tectonic,
lithospheric anomalies. The methodology which is presented could be
used to prepare more detailed maps of seismic risk for wider areas,
even though seismological data of such areas are missing or are of
limited extent. This is an advantage, compared to seismological
methods which require a statistically large number of earthquakes to
take place in order to evaluate the seismic risk. The case of the EQ
in Grevena (seismic risk zone I) and the EQ in Athens (seismic risk
zone II) are two characteristic examples, which justify the validity
of this methodology.
In conclusion, taking into consideration the latest, increased,
seismic activity of the Greek territory, it is justified to expect,
strong earthquakes to occur along this faulting / fracture-zoning
network in the future. Therefore, the State Authorities must take
into consideration the fact that the seismic risk is enhanced along
these specific fracturing / faulting zones. The fact that, in some
of them we have not experienced a strong earthquake yet, does not
justify, at all, the notion that it is an aseismic zone. On the
contrary, it must be considered as a seismogenic area, potentially
to be activated in the future. What is open, as a question, that is
only the timing of the future strong EQ.
2.2. Depth distribution of strong earthquakes.
A parameter of an earthquake of great importance is its depth of
occurrence. The vast majority of the seismicity of a seismic region
takes place in the crust (fig. 2.2.1) of the Earth. The crust, as it
is shown in the sketch presentation of the next figure (2.2.1), is
very thin compared with the other discrete parts of the Earth’
interior (mantle, outer-inner core).
2.2.1. The structure
of the Earth (www.physicalgeography.net). The crust, mantle,
outer – inner core of the Earth’s interior (not to scale).
A more detailed presentation of the lithospheric structure is
presented in the following figure (2.2.2). The continental and
oceanic crust is shown, in relation to a larger geological unit,
referred, as “the lithosphere”, which overlays the “plastic
asthenosphere” and the “upper mantle”.
2.2.2. The structure
of the lithosphere (www.physicalgeography.net). The
continental – oceanic crust, the lithosphere, the plastic
asthenosphere as long as their representative depth extents
are shown (not to scale).
The characteristic depth values that are assigned to each of them
have been identified by specific studies upon the change of their
physical properties. The analysis of the P velocity of the seismic
waves as a function of depth (Mueller and Landisman, 1966) has
revealed the existence of a velocity decrease zone at a depth around
10Km, while an abrupt increase is observed at depth of 30Km (Moho
discontinuity). The latter is presented in figure (2.2.3).
2.2.3. Model of P
velocity distribution within the crust, that shows the
postulated low velocity layer in the upper crust (Mueller
and Landisman, 1966) and its abrupt increase at Moho depth.
Toksoz et al. (1967) studied the shear wave velocity as a function
of depth. The results of this study are presented in the following
figure (2.2.4). The main feature of this study, concerning the depth
of occurrence of earthquakes, is the shear wave velocity decrease,
observed, at a depth of 100Km in tectonic regions.
2.2.4. Upper mantle
shear velocity models, for oceanic, continental shield and
tectonic regions, based on dispersion of surface waves. The
profiles are uncertain below 500Km depths, because of
insufficient data (Toksoz et al. 1967).
The majority of the EQs occur in the seismogenic zone of the
lithosphere that is where earthquakes mostly take place. The base of
the seismogenic zone is the top of the more ductile asthenosphere. A
simplified model of an earthquake and its associated rupture surface
is shown in figure (2.2.5).
2.2.5. A simplified
model of an earthquake is shown and its associated rupture
surface (S). The dimensions of the rupture are denoted as W
and L for its depth and length extent accordingly. S =
rupture surface and D = direction of rupture.
The detailed study of a seismogenic area and particularly the study
of the depth distribution of the earthquakes which occur in that
area, for a rather large period of time, reveal the depth extent of
its corresponding, seismogenic zone.
This is illustrated in the following figure (2.2.6), as an example
that refers to the depth distribution of the foci of EQs,
registered, at the Nicoya peninsula, Costa Rica (Avants et al. 2001,
Newman et al. 2002).
distribution of the foci of EQs registered at Nicoya
peninsula, Costa Rica that illustrates along-strike
variability in the seismogenic zone seismicity (Avants et
al. 2001, Newman et al. 2002).
The seismogenic zone of the Nicoya Peninsula revealed from the
up-dip events extends from 10 – 20Km, while the analysis of the
total events (around 2500) indicates that the seismogenic zone
starts from almost 10Km and extends up to 50Km.
The physical parameters of the lithosphere, generally, change along
depth and therefore, the mechanisms that generate earthquake
precursory phenomena depend upon the depth of occurrence of each EQ.
Consequently, earthquakes of different magnitude which occur in
different depths, probably will generate different precursory
phenomena and of variable intensity.
In the Greek territory, the depth distribution of earthquakes is
analyzed in an introductory form as follows:
Number of EQs, as a function of occurrence depth, for a
fixed magnitude, in steps of .5R
(5.5 – 7.5 R).
Number of EQs as a function of magnitude for a fixed
depth, in steps of 5 Km (0 – 10Km).
In both cases the EQ catalog of NOA, which was used, spans from
1901 to 2006. However, it must be pointed out that the first part of
this catalog (1901 – 1950) refers to the early period of earthquake
parameter estimation and therefore, because of technological lack
during this period, a depth of 0Km was assigned to each registered
EQ. Therefore, EQs which refer to a depth of 0Km will be ignored.
2.2.1. Number of EQs as a function of occurrence
depth for a fixed magnitude of
Depth distribution of the EQs in Greek territory with magnitude of
5.5 R is shown in the following figure (126.96.36.199).
distribution of EQs in Greek territory, with magnitude of
Four distinct peaks are characteristic in this diagram. The first
is located at 0Km depth and will be ignored as was mentioned
earlier, the second one is located at a depth of 10Km, the third one
spans from 30 to 40Km, while a forth one is located at 100Km. Single
EQs extend to a depth of 150Km.
2.2.2. Number of EQs as a function of occurrence
depth for a fixed magnitude of
In the next drawing (188.8.131.52) has been used instead a magnitude of
distribution of the EQs in Greek territory, with magnitude
For the magnitude of 6.0 R the depth distribution of the EQs is
similar to the 5.5R one, but with lower peak values. The depth zones
of 10Km, 30-40Km, and 100Km are characteristically the same but an
extra peak is shown, at a depth of 120Km.
2.2.3. Number of EQs as a function of occurrence
depth for a fixed magnitude of
For the next magnitude level of 6.5 R (fig. 184.108.40.206) the No of EQs
is less than the previous case, but the very same depth zones are
revealed and one more at a depth of 65Km, as well.
distribution of the EQs in Greek territory, with magnitude
It must be taken into account that EQs of such a magnitude are rare
in the Greek territory and therefore, even a single seismic event of
such a magnitude is very important for the iden-tification of these
The latter is validated by the depth distribution of EQs with
magnitude of 7.0 R and 7.5 R, shown in the following figures
(220.127.116.11 – 18.104.22.168).
2.2.4. Number of EQs as a function of occurrence
depth for a fixed magnitude of
distribution of the EQs in Greek territory, with magnitude
2.2.5. Number of EQs as a function of occurrence
depth for a fixed magnitude of
distribution of the EQs in Greek territory, with magnitude
The study of the depth distribution of the EQs for the range of
magnitudes from 5.5 R to 7.5 R indicates the existence of
preferential depths in the lithosphere, where the EQs take place.
2.2.6. Number of EQs as a function of magnitude for a
fixed depth of 0, 5, 10Km.
22.214.171.124. Number of EQs
as a function of magnitude for fixed depths of 0, 5 and
In this graph is made a comparison between the No of EQs which
occur at an intermediate depth of 5Km and the depth of 10Km. It is
evident that, in their majority, the EQs present a tendency to occur
at the depth of 10Km of the seismic zone.
2.2.7. Number of EQs as a function of magnitude for a
fixed depth of 100Km.
Finally, a similar analysis has been made for the depth of 100Km.
126.96.36.199. Number of EQs
as a function of magnitude, for a fixed depth of 100Km.
Although the number of EQs is rather small, at the depth of 100Km,
it is evident that mainly (more than 50%) strong events occur at
The seismic zones, identified by the previous diagrams, coincide
with the depth of P seismic (Fortsch and Moho) velocity
discontinuities, observed, by Mueller and Landisman (1966) and with
the S shear, seismic wave velocity discontinuity, observed, at a
depth of around 100 Km, (Toksoz et al. 1967).
In conclusion, it is expected that strong EQs which take place in a
seismogenic area behave, probably, differently, as far as it
concerns the generation of precursory phenomena. Therefore, it is
possible that some times precursory phenomena may not be generated
or these are undetectable with the present state of physical
knowledge and technology, available to the scientists.
2.3. Seismic energy density at distance (x) from the
During the occurrence of a strong EQ a large amount of strain
energy, accumulated and stored in the regional focal area, is
released. That energy propagates from the focal area outwards, in a
more or less spherical surface mode, until it reaches a seismic wave
velocity discontinuity. In such a case it can be reflected,
refracted or diffracted, depending on the specific physical -
tectonic conditions, met in the area of incidence.
The energy, which is released through an earthquake, obeys the
simple laws of physics. The application of these laws in seismology
explains, in a very simple way, some well known notions as: a) a
rocky ground is safer than a sedimentary one, b) different ground
acceleration is observed in very close distances. These phenomena
and observations can be explained by energy manipulation in terms of
physical laws, as it will be explained in the text to follow.
The actual, local, seismic intensity felt, due to a distant
earthquake, depends, in its simplest approach, on the energy /
surface unit that “in-flows” at this area. Assuming homogeneous
Earth, the largest the distance between the focal area and the
affected area is, the less the seismic intensity is felt in this
area. This is demonstrated in the following figure (2.3.1).
energy that outflows at a solid angle from the focal area
and passes through S1 and S2 remote surfaces of different
locations and referring to the same solid angle.
Assuming that the outflow energy in a solid angle is (E) and
S1 = ðR12, S2 = ðR22
I1 = E/S1 and I2 = E/S2
Where I1, 2 are the seismic energy per surface unit, in other words
the intensity felt at each area 1, 2.
Substituting S1, 2 in equation (2.3.2) by the expressions of
equation (2.3.1) we get
I2/I1 = (R1/R2)2
And by taking into account that R = (solid angle/2) * x
Where: x stands for the distance of each area from the hypocenter
of the EQ, then the equation (2.3.3) results into:
I2 = I1*(x1/x2)2
Equation (2.3.4) indicates that the seismic intensity, which is
felt at two different sites from the same EQ, is inversely
proportional to the square of the ratio of their distances from the
hypocenter of the EQ.
The same is true for the acceleration which is observed in an area,
registered, during the occurrence of an earthquake. However,
sometimes the accelerogram indicates larger acceleration than what
is normally expected. This is due to local geological conditions and
is especially caused by lateral, geological discontinuities, when
the seismic energy propagates through material of different
A simple physical - mathematical explanation is as follows:
It is assumed that a discontinuity exists between a geological
formation of d1 density and another one of d2. Adjacent the
discontinuity, we consider the same unit volume M in both geological
formations (fig. 2.3.2).
discontinuity sketch model. Seismic energy flows from medium
of density d2 towards medium of density d1. At the
discontinuity an equal volume M is activated in both sides.
It is assumed that there is no energy loss during this
As long as seismic energy propagates, along the geological
formations, the same quantity of energy flows through the adjacent
unit volumes at the discontinuity. The latter is expressed by the
Seismic energy = ½ d2M(V2)2
Seismic energy = ½ d1M(V1)2
Where V1, 2 are the motion velocities of the two (M) volume masses
of the different density media.
By equating the second part of equations (2.3.5) it results into:
V1 = V2*(d2/d1)1/2
Indicating that, for the same amount of energy transfer through a
lateral discontinuity, the relative motion velocity, in the two
different media, depends on the square root of the ratio of their
densities. This is the physical explanation of the “amplification”
of the seismic waves, observed, over loose materials. It is nothing
more than the utilization of the conservation of energy which flows
from one geological formation to the other, under the assumption
that there are no other significant energy losses of any kind.
Furthermore, the seismic energy is transferred through an
oscillating seismic wave of angular velocity:
Vù = V0*sin(ùt)
Therefore, the ratio Vù1 / Vù2 of the angular velocities of the two
media is expressed by the equation:
Vù1 / Vù2 = [V01* sin(ùt)] / [V02* sin(ùt)]
Vù1 / Vù2 = V01 / V02
Where: V01 and V02 are the maximum values of the velocities,
observed, in the two geological formations.
V01 ≠ V02
dV01/dt ≠ dV02/dt
The latter equation (2.3.10) indicates that the acceleration which
will be observed in the two geological media will be different for
the same activating seismic wave. Moreover, the difference in the
observed acceleration will depend on the ratio of their densities,
as it is indicated by equation (2.3.6).
The already presented topic is related to the seismological notion
which more or less suggests: “a large earthquake which occurs in the
sea is not so dangerous for the nearby towns, since most of the
energy is absorbed by the seawater”. The latter is presented in the
following figure (2.3.3).
diagram that represents the seismic energy propagation from
hypocenter, through its traveling wave fronts, towards a
The seismic energy that arrives to the town region is larger than
what was absorbed by the seawater. This is demonstrated in figure
(2.3.3). The energy that travels along the path EQ hypocenter – town
and is carried by the seismic wave fronts is larger than what can be
absorbed by the thin (compared to the EQ – town distance) seawater
layer. Moreover, the S waves that carry part of the total energy,
cannot propagate in the seawater. The vast majority of the seismic
energy propagates along the referred EQ – town line path through
solid ground. Therefore, all the latter physical laws are valid and
In conclusion, a strong earthquake can be destructive and always
dangerous for the nearby towns, regardless of its epicentral
location in seawater or on solid ground.
2.4. Electrical resistivity lithospheric model.
In section (2.2) the physical parameters of P and S seismic wave
velocity were presented as a function of depth in the crust and
upper mantle. Although these parameters prescribe more or less the
physical behavior of the medium, another property of it, the
resistivity, prescribes its electrical behavior.
It is well known, in the geological – geophysical sciences, that
each geological formation presents a range of resistivity.
Therefore, it is possible, by obtaining knowledge of the depth
distribution of the resistivity of a geological area, to assign to
it the most probable geological formations. In areas which are more
tectonically complicated it is possible to model complex tectonics,
just by using the appropriate electrical methodology. Details on
this topic can be found in any textbook for Applied Geophysics in
the corresponding chapter of the application of electrical methods.
The depths, involved for the study of the crust and the upper
mantle, are of the order of some decades of Kilometers. Therefore,
the magnetotelluric methods were used, initially, in order to
achieve large depth penetration of the current, used, for this kind
of studies. The key feature of these methods is the depth
penetration (skin depth) of an electromagnetic wave in the ground
which depends on its frequency and the resistivity of the ground.
The latter is expressed by the following equation (2.4.1).
Z = 500*(ñ/f)1/2
Where (Z) is the penetration depth in meters, (ñ) is the ground
resistivity in Ohm*m and (f) is the frequency of the electromagnetic
wave in cycles/second. A main drawback of these methodologies is the
limited band pass (T<1hour) of the filters, used for registering the
electric and magnetic field, while at the same time, for larger
periods signals, a very long time is required for an effective
registration. An alternative to this methodology is the traditional
Schlumberger electrical sounding of the AMNB electrode
configuration. In such a kind of electrical arrays the length AB of
the electrical dipole which induces the DC current into the ground,
controls the current depth penetration and therefore, the depth of
Results obtained from such a kind of operations were reported by
Blohm and Flathe (1970) after having used a 150Km AB electrical
dipole, Blohm (1972) and Homilius and Blohm (1973) in the Rhine
Graben area. Other researchers, using a larger AB dipole length (up
to AB = 600Km) in a different tectonic province (Southern Africa)
investigated the resistivity distribution in the Earth to larger
depths. In this frame of work Van Zijl (1969), Van Zijl et al.
(1970) studied the crustal conductivity structure of South Africa,
Van Jijl and Joubert (1975) for the same purpose used AB length of
450Km and 400Km. In 1976 Van Jijl reported results, obtained, from
30 deep electrical soundings (AB spacing = 40Km) carried out at the
Umtali – Pietersburg area (South Africa) while Blohm, Worzyk and
Scriba (1977) presented results, obtained, from a Shlumberger
sounding with AB = 1250Km which is the longest used, so far, and
known in the literature that deals with this geophysical
The unique graph of the apparent resistivity as a function of AB/2,
obtained, from this experiment, is presented in the following figure
resistivity as a function of depth (AB/2), obtained, from
the application of the electrical sounding methodology in
Southern Africa by using the Cabora Bassa power line (Blohm,
Worzyk and Scriba, 1977) as a power source.
The transformation of the apparent resistivity depth function into
the corresponding resistivity layers revealed the following
resistivity layering, presented as model-1 and model-2, which models
are both valid for the area of the study. The first model, model-1,
is based on a four-layer composition of the lithosphere, as it is
presented in the following figure (2.4.2).
resistivity model – 1 (after Blohm et al. 1977).
A more detailed and equally valid model is the following model-2,
which is composed by five layers (fig. 2.4.3).
resistivity model – 2 (after Blohm et al. 1977).
By taking into account the equivalent principle on the
interpretation of the resistivity graph, it is possible to accept
more different “resistivity interpretations”, even composed with
more electrical layers, but these two have the following, general
A resistive zone with a maximum thickness of 24Km is followed by an
intermediate, conductive layer which has a thickness of about 18Km.
The next layer exhibits a high resistivity and a thickness of 120Km.
The resistivity of the last layer is not defined very well, but it
is definitely lower.
The resistivity of the upper crustal layer seems to be rather high
(ñ>=50000 Ohmm), while the conductive layer, found, in the lower
crust, coincides with the results, obtained, by the use of the
magnetotelluric method (Keller, 1971; Adam, 1976). The resistivity
of the upper mantle (8000 Ohmm) is in reasonable agreement with
laboratory results (Brace, 1971), while the thickness of the
lithosphere (160Km) agrees very well with the values 150 – 175Km,
obtained, by teleseismic delay times (Fairhead and Reeves, 1976).
The final layer has low resistivity (ñ<= 50 Ohmm) in accordance to
the results, obtained, by the deep magnetic sounding method (Schmucker,
A comparison between (P) and (S) velocity depth distribution and
the resistivity in the crust indicates that, there is a good
agreement of Moho depth, suggested by the P velocity distribution
and the thickness of the resistive layers of the upper crust.
Specifically, Moho depth is clearly indicated by the bottom of the
conductive layer of the 3rd layer in model-2, while the (S) velocity
depth distribution coincides, within the first 200Km, with the high
resistant package of the crust.
Summarizing these results, it is evident that the seismogenic zone
in the crust behaves, particularly, as a highly resistant medium in
terms of its electrical properties. The latter, combined with its
mechanical properties, is very important in the overall development
of the precursory phenomena which are used in earthquake prediction
A final point which must be mentioned is the generalized electrical
model of the crust. In terms of geological formations, the very few
first kilometers of the crust from the ground surface, consist of
mainly sedimentary formations. This implies that the top few first
kilometers of the crust consist of highly conductive material,
followed, by highly resistant layers up to a depth of about
160-170Km and finally, followed, by a highly conductive basement.
This is presented in the following figure (2.4.4).
crust conductivity model.
This electrical conductivity sequence can be viewed as an
electrical system which is composed on an electrical shield (top and
bottom conductivity) which “shields” the inner resistant conductor
from external, electrical, electromagnetic influences. The latter
will be speculated about in more details in the sections to follow.
2.5. Physical models, used, in this methodology.
The issue of the earthquake prediction has raised strong debates in
the scientific community. One of the arguments, which are presented
against it, is the “absence of a physical model” that will justify
any methodology to this end. In other words, there is no physical
mechanism, universally agreed, for the seismogenesis.
The different physical models, which have been proposed up to date,
such as the “Rebound Model” (Reid, 1911), the “Self Organized
Criticality (SOC) model (Bak et al. 1988, Bak 1996), the Chaotic
Non-linear Systems (Anderson 1990, Kagan 1994, Main 1996), highly
depend on the initial conditions, thus suggest the unpredictability
of EQs. Other studies have shown that a small EQ can grow into a
strong one, depending on very small variations on elastic
properties, fault parameters variations as friction, stored, strain
energy (Otsuka 1972a, b; Vere-Jones1976, Ito et al. 1990).
Moreover, the fault friction plays the most important role in
seismogenesis, in the frame of the tectonic regime of the
lithosphere, while fracturing is a secondary one (Scholz, 1998). The
term “stick-slip frictional instability”, introduced, by Brace and
Byerlee (1966), suggests that strain energy is accumulated during
the “stick” period, while an EQ occurs at the “slip” period.
Actually, the phenomenon of seismogenesis is a mixture of frictional
slip failure and shear fracture (Ohnaka, 2003).
The seismogenesis mechanisms, postulated to date, refer to the
physical-mechanical procedures that take place in the focal area,
before and during an EQ occurs. In a most indirect way they refer to
EQ precursors and therefore, to the capability of prognostic
parameters determination. An example of such a calculation is the
“power law - time to failure model” that correlates the magnitude
(M) of a future strong EQ to the remaining time (tr) towards its
occurrence (Bufe and Varnes 1993, Bowman et al. 1998).
The absence of a valid and robust relation between the different
seismogenic mechanisms and the seismic precursors, required, for the
determination of the prognostic parameters (location, magnitude,
time), led the seismologists to apply statistical methods for the
issue of the earthquake prediction. Methodologies, as the algorithms
M8 (Keilis-Borok et al. 1990, Healey et al. 1992, Romashkova et al.
2002) and CN (Keilis-Borok et al. 1990a), have been applied with
some success, mainly for long-term and medium-term prediction for
rather large areas.
These statistical and any other prognostic algorithms fail to
satisfy the postulated logical equation (1.1). At this point it must
be stressed out that, what is really needed for a successful
earthquake prediction is a physical model upon which the appropriate
calculations will be based in order to analyze the precursory
available data, so that earthquake prediction can be implemented.
The difference between these physical models and the physical,
seismogenic mechanisms is essential. The seismogenic mechanisms
refer to what actually happens in the seismogenic region and
therefore, they are a physical “close-up” view in the seismogenic
area, while the physical models, used, for the calculation of the
prognostic parameters, are a more generalized, physical approach of
the seismogenic region, and they are independent from the actual
seismogenic mechanisms which take place in it.
Moreover, it is anticipated that, a single physical mechanism
cannot provide answers for all earthquake prognostic parameters.
This is evident from the large number of publications, which are
related to the topic of earthquake prediction. The majority of them
refer only to time, or regional area in relation to magnitude. None
of them deals, simultaneously, with all the prognostic parameters.
In the methodology to be presented, a different approach was
followed. During the course of this research (1981 – 2003) it was
found that each prognostic parameter required, a different
generalized physical model, to be used. Each one of them fulfilled
the logical equation (1.1). Therefore, the simultaneous use of them,
applied on the appropriate precursory data, facilitate the
implementation of the sort-term earthquake prediction.
These models are the following:
a. The lithospheric seismic energy flow model. This is a direct
application of the energy conservation law of physics. The physical
system, in concern, is the lithospheric, seismogenic region. This
model is used for the calculation of the magnitude of an imminent,
strong EQ. Furthermore it explains, in energy transfer terms through
the lithosphere, the “accelerated deformation” and “seismic
quiescence” methodology, used, by the seismologists.
b. The oscillating lithospheric plate model. Lithosphere is treated
as a plate, which is forced into oscillation by the tidal forces.
This methodology is in use by geophysicists, for the correction of
gravity measurements, due to errors, induced, by the tidal
oscillation of lithosphere. The latter is used for the time of
occurrence determination in conjunction with the homogeneous ground
c. The homogeneous ground Earth model. This is the start-up model
in the analysis of the electrical methods, used in Applied
Geophysics. It is used for the determination of the epicentral area.
It must be mentioned, here, that all these models have been long
ago, well-known, physical models, applied on Earth, for the
application of geophysical studies for different geological and
The detailed use of each one of these physical models will be
presented in the prognostic parameters calculation.
2.5.1. The lithospheric seismic energy flow model.
188.8.131.52. Theoretical analysis.
It is generally accepted that stress energy built-up, in a focal
area, is a very slow process which, closely, follows the motion of
the lithospheric plates. It takes a long period of time (probably a
large number of years) to reach the point when an earthquake will
occur, because of rock fracturing.
Under normal conditions, the stored energy is discharged through
the background, small magnitude, seismicity of the seismogenic area.
The latter is demonstrated in figure (184.108.40.206.1).
postulated model of seismic energy flow, through a
seismogenic area of lithosphere.
When a strong earthquake is in preparation, before its occurrence,
the normal seismicity either decreases for a certain period of time
and therefore, the seismic “quiescence”, which is detected, is used
as a precursory indicator or increases and therefore, the
“accelerated deformation” is observed.
The mathematical analysis of the postulated model is as follows:
Let us denote as (Ein) the inflow energy over a short period of
time (dt) in a seismogenic area of the lithosphere and as (Eout) the
energy outflow of the same area, due to its seismicity for the same
period of time. (Ein) and (Eout) are time functions, Ein = Ein(t)
and Eout = Eout(t).
The seismogenic area is charged with energy Est = Est(t) as
For any very short period of time (dt), equation (220.127.116.11) takes
the following form:
Est(dt) = d(Est(t))
For successive values of dt, a discrete valued function Y(t) is
defined as follows:
Y(t) = Est(t)
Equation (18.104.22.168) can take the following forms:
a. Y(t) = C = 0
In this case, the stored energy in the seismogenic area, equals to
zero, therefore, it is at a state of stable zero-charge conditions.
b. Y(t) = C < 0
In this case, the stored energy decreases, therefore, the
seismogenic area continuously discharges, towards a very stable
state of uncharged conditions.
c. Y(t) = C >0
In this case the stored energy increases continuously, therefore,
the seismogenic area is charged towards a state of highly unstable
conditions, leading to the occurrence of an earthquake.
The term (C) of equation (22.214.171.124) can generally be either time
dependent or time independent.
- C is time independent:
Equations (126.96.36.199) and (188.8.131.52) can be combined in equation
d(Est(t)) = C >0
The cumulative energy (Ecum), stored in the seismogenic area, can
be calculated as a function of time (Ecum(t)) by integrating both
sides of equation (184.108.40.206), in respect to time (t),
Ecum (t) =
Since term C is constant and time independent, the calculated
function Ecum(t) has the form of:
Ecum (t) = C*t + b
Where, b denotes the integration constant.
The linear equation (220.127.116.11) was firstly introduced by
Thanassoulas et al. (2001), along with the postulated, lithospheric
seismic energy flow model, for the calculation of the maximum
magnitude of an imminent, strong EQ.
- C is time depended:
In this case, equation (18.104.22.168) can be represented by an nth-order
d(Est(t)) = antn + an-1tn-1 + ......... + a0
The cumulative energy Ecum(t), stored, in the
seismogenic area, can be calculated as a function of time by
integrating in time, both sides of equation (22.214.171.124).
ò( antn + an-1tn-1 + ......... + a0
)dt (2.5.1. 9a)
As an n+1 order polynomial,
Ecum(t) = kn+1 tn+1 +kn tn +………k0
Where, ki values represent the polynomial constants.
Summarizing the forms the equation Ecum(t) (that represents the
energy flow) takes, the following cases are possible:
a. linear polynomial - constant energy flow
b. higher order polynomial - real acceleration
c. accelerated for a period of time long before the main, seismic
event, which is followed by a constant energy flow , just before the
occurrence of the strong EQ.
Cases (b) and (c) are represented, by the well-known “time to
failure” function, very often.
Although, mathematically, it is possible, to transform any polynomial to
any arbitrary function, i.e. time to failure function by calculating
the appropriate parameters of the latter, by using LSQ techniques,
still remains the parameter of arbitrarity, as far as it concerns
the validity of physics behind this transformation.
Moreover, the time to failure function depends on two variables.
The first one is the magnitude of the imminent EQ and the second one
is the time to failure, left. In order to overcome the problem of
solving a two parametric equation (infinite number of solutions),
the parameter C (Bowman et al. 1998) was introduced, that is the
ratio of power low fit error over the linear fit error, as far as it
concerns the cumulative seismic energy release. Still, the notion of
this ratio is set completely, arbitrarily.
Therefore, it is suggested that the magnitude and time to failure
of a strong, imminent EQ, calculated, by these methodologies, are
not supported by any validated, physical mechanism and should be
126.96.36.199. Application of the model on real EQ cases.
The application of the lithospheric, seismic energy flow model
requires two parameters to be known, in advance. The first one is an
initial estimation of the area extent of the physical system itself
(seismogenic area). The second one is its seismic history.
The choice of the first parameter is utilized by a) knowing the
epicentral area of the future EQ by another methodology, as it will
be explained later on b) by taking into account the sometimes
observed increased low level seismicity of an under study
seismogenic area, which is some times observed and finally c) by
taking into account the deep, lithospheric fracture zones which are
mapped by the conversion of the regional gravity map of the seismic
region in to a gradient one.
This procedure is explained by the figures to follow.
As long as the seismogenic area has been activated, a small
magnitude seismic activity, associated with the main rupture of the
rock formation, increases in general, in the different order
fracturing branches (nth order ridel) and therefore, this seismic
activity is observed along and very close to the trace of the main
seismic to be activated fault. The latter is presented in the
following figure (188.8.131.52.1).
184.108.40.206.1. The main,
seismic fault branches into smaller ones, as it approaches
the ground surface (Mattauer 1973, Vialon et al. 1976). The
group of fractures that advances up to the ground surface
forms the observed “fracture zone”.
Since small-scale seismicity has started to emerge, the
corresponding foci will be located close to the vicinity of the main
fault, which will be activated in the future. A sketch drawing
indicating this main earthquake precursory activity is presented in
the following figure (220.127.116.11.2).
drawing of the theoretically small magnitude seismicity,
expected, to be observed on ground surface and in the
vicinity of a main, seismic fault, prior to a strong EQ. The
solid dots represent the epicenters of the precursory,
seismic activity, while the thin lines indicate the nth
order “ridel” different faults. The thick, black line
represents the main fault which is going to be activated.
Consequently, the physical system which will be studied must be
chosen, close, to the main fracture, in a way that it takes into
account this specific distribution of the small magnitude EQs. The
second parameter, the seismic history of the seismogenic area, is
taken from the EQ files from the Seismological State Observatories.
In case of Greece, this is the National Observatory of Athens (NOA).
The application of the methodology is shown in the following
figures. It is assumed that, the epicentral area has been
approximately determined by other means. As a first step, is defined
the seismogenic area which will be considered as the corresponding
physical system of which the seismic energy release will be analyzed
approximately, determined seismogenic area is indicated by
the blue polygon. It has been drawn in such a way to center
the associated deep, main, lithospheric fracture zone and is
elongated along it.
The next step to take is to calculate backwards the cumulative
energy release of this area for a long period of time. The graph
which results indicates whether this area has entered the
acceleration phase or not. The latter is presented in the following
seismic energy release calculated, for the seismogenic area,
which is indicated in the previous figure (18.104.22.168.3) for
the period 1998 – 2006.
The study of this graph indicates that this example, seismogenic
area has been set into cumulative, seismic energy release
acceleration mode for the last 8 years. A slightly different
presentation is shown in next figure (22.214.171.124.5) by fitting a 6th
degree polynomial. This facilitates the analytical calculation of
the time gradient of the cumulative seismic energy release.
seismic energy release fits with a 6th degree polynomial.
The black line indicates the cumulative, seismic energy
release data, while the red one indicates the fitted,
It is obvious that the acceleration has been initiated from the end
of 2005, indicating that an earthquake will strike “soon”. Actually,
a strong (M = 6.9R) earthquake took place at a short distance
towards NW on this lithospheric fracture zone at the start of 2006.
Its details will be presented in the examples to follow.
Quite often, the start of the increase of the cumulative energy
release is not defined very sharply, but there is a gradual change
over a rather lengthy period. In such cases, instead of fitting only
a polynomial, as an advanced step, the corresponding gradient of the
polynomial is calculated analytically.
The process of gradient calculation acts as a high-pass filter on
the cumulative seismic energy release data, low-order terms of the
fitted polynomial are excluded by this operation and therefore the
resolving capability for detecting smaller level changes is larger.
This process is shown in the following figure (126.96.36.199.6).
seismic, energy release time-gradient, calculated
analytically, from the polynomial of figure (188.8.131.52.5).
184.108.40.206.1. Application of the theoretical model on
the EQ in Zakynthos (02/12/2002,
Ms = 5.8R).
The first example, based, on the methodology which has been already
presented, is that of Zakynthos EQ. A few months before the main,
seismic event which took place on 02/12/2002, with a magnitude of Ms
= 5.8R, an increased, small-scale, seismic activity had been
observed at Zakynthos area. The area of interest is indicated by a
red circle in the following figure (220.127.116.11.1.1).
18.104.22.168.1.1. The red
circle indicates the regional area of Zakynthos. The blue
circles indicate the small magnitude seismic activity all
over Greece for a time period of a couple of months before
the EQ occurred.
The seismicity which was observed to increase, dictated the
application of the latter methodology. The area, to be investigated,
was selected in such a way, to include the main, lithospheric
fracture zones, which are present in the region. The latter is
presented in the following figure (22.214.171.124.1.2). The blue frame
indicates the area of interest.
126.96.36.199.1.2. Area of
interest (blue frame) for which the lithospheric seismic
energy release methodology will be applied. The brown lines
indicate the deep lithospheric fracture zones.
Next step is to calculate backwards the cumulative seismic energy
release for some years. The latter was utilized by using the EQ
file, which is available online in the web by the NOA, Athens,
Greece. This operation indicated that the area of Zakynthos was set
in seismic acceleration mode from the start of the year 2000. This
is presented in the following figure (188.8.131.52.1.3).
seismic energy released, from the seismogenic area of
Zakynthos for the period 2000 – 2002.
The 6th order polynomial, fitted in the cumulative, seismic energy
release data for the area of Zakynthos, is presented in the
following figure (184.108.40.206.1.4).
220.127.116.11.1.4. The 6th order
polynomial fitted, in the cumulative, seismic, energy
release data for the area of Zakynthos. The brown line
indicates the cumulative, energy data, the blue one
indicates the fitted polynomial. The red arrow indicates the
time when the Zakynthos EQ occurred.
The analytical expression of the polynomial which is obtained
through the LSQ fitting procedure, utilizes the calculation of the
time gradient of the polynomial and therefore, the rate of change of
the cumulative, seismic energy release. This is presented in the
following figure (18.104.22.168.1.5).
22.214.171.124.1.5. Rate of
change in time of the cumulative, seismic energy release
observed, prior to the EQ of Zakynthos. The red arrow
indicates the time when this EQ occurred.
It must be pointed out that all this analysis had been made before
(one month) the EQ occurred and therefore, this EQ was expected to
occur soon. The, expected, EQ did happened within a month and its
location coincides with the location of the deep lithospheric
fracturing which is located in this area (fig. 126.96.36.199.1.6),
188.8.131.52.1.6. Location of
the EQ in Zakynthos in relation to the deep, lithospheric
fracture zones, which are identified in the same area, by
analyzing the gravity field.
thus verifying the validity of the corresponding, deep,
lithospheric fracture zones / faults map (Thanassoulas, 1998).
184.108.40.206.2. Application of the theoretical model on
the EQ in Kythira (08/01/2006,
Ms = 6.9R).
This example is an “a posteriori” one. The validity of the method
was tested against Kythira EQ (08/01/2006, Ms = 6.9R), which is
presented in the following figure (220.127.116.11.2.1).
(blue concentric circles) of the EQ in Kythira (08/01/2006,
Ms = 6.9R), in relation with the location of the deep,
fracture zones / faults of the lithosphere.
The seismogenic region, which is taken into account, is presented
as a blue polygon frame in the figure (18.104.22.168.2.2) below.
22.214.171.124.2.2. Area of
interest (blue frame) for which will be applied the
lithospheric, seismic energy release methodology. Brown
lines indicate the deep, lithospheric fracture zones.
The cumulative, seismic energy release, which is calculated for a
period of 14 years (1992 – 2006), indicates that this seismogenic
region was set in acceleration mode for the last 2 years (2004 –
2006). This is presented in figure (126.96.36.199.2.3).
seismic energy release determined, for the seismogenic
region of Kythira, for the period 1992 – 2006.
The fitted, 6th order polynomial function of cumulative, released,
seismic energy vs. time, indicates a rapid increase “graph knee”,
almost 2 years before the corresponding earthquake occurred. This is
presented in next figure (188.8.131.52.2.4).
seismic energy release is fitted by a 6th degree polynomial.
The red line indicates the cumulative, seismic energy
release data, while the blue one indicates the fitted,
polynomial values. The black arrow indicates the time when
the EQ occurred.
In this case, the time gradient, which is calculated from the
analytical expression of the polynomial, resolves the time, when the
seismogenic area entered the period of intense seismic acceleration
in a much better way. The latter is shown in next figure
184.108.40.206.2.5. Rate of
change, in time, of the cumulative seismic energy release
observed, prior to the EQ in Kythira. The black arrow
indicates when the EQ occurred.
A point that must be clarified, more, is the way it is set the
seismogenic area extent. What is accepted to date in the
seismological community is that, the stronger the pending EQ is, the
larger the area which is overcharged with strain deformation around
the epicenter is. The latter is represented by the following
log R = 0.42M – 0.68
where: (R) is the radius of the strain charged area and (M) is the
magnitude of the pending earthquake (Papazachos et al. 2000,
Papazachos et al. 2001).
A radius of almost 200Km has been estimated, after having adopted
this formula, for the assignment of the seismogenic region which is
going to be used for the calculation of the seismic energy released
for the case of the EQ in Kythira (Ms = 6.9R). As a first approach,
this area is approximated by a surface, enclosed in a, closely,
orthogonal frame of dimensions: 4 degrees (latitude) by 5 degrees
(longitude). The latter is presented in the following figure
area (blue frame) theoretically charged with strain
deformation, required for the generation of the EQ in
Kythira (Ms = 6.9 R).
It is clearly obvious that the seismogenic area, determined, by the
use of equation (220.127.116.11) will be affected, not only by the main
fault where this strong event took place, but also by the adjacent
The latter will contribute, with their seismic energy release, to
the one, calculated for the entire frame. We must avoid it, when
such calculations are made, since the seismic energy release of the
fault which is expected to be activated, is studied. Moreover,
calculations made for the identification of the accelerated
deformation are masked by the presence of faults, which will not be
activated at all. This is demonstrated in the following figures.
Let us assume a 2 by 2 degrees frame (fig. 18.104.22.168.2.7) that
encloses some more main fracture zones, besides the one that was
22.214.171.124.2.7. A 2 by 2
degrees seismogenic area (blue frame), initially, assumed,
for the cumulative, seismic energy release determination.
The cumulative seismic energy release is calculated as a function
of time, for the same period (1992 – 2006) as the EQ of Kythira.
This is presented in the next figure (fig. 126.96.36.199.2.8).
cumulative, seismic energy release calculated, for the same
period (1992 – 2006) as for the EQ in Kythira.
This graph fits quite well a straight line, except for the period
1998-2002, when a moderate magnitude seismic event took place.
There is no evidence to show that any, accelerated, deformation
mechanism has been initiated. The same is obvious from next figure
fitting (6th order), as determined on the cumulative seismic
energy release data.
For the next example, was used a smaller frame (1 by 1 degree),
which encloses only the activated area, but does not take into
account the strike of the seismically, active fault. This frame is
presented in next figure (188.8.131.52.2.10).
184.108.40.206.2.10. A 1 by 1
degrees seismogenic area (blue frame), initially, assumed,
for the cumulative, seismic energy release determination.
The corresponding, cumulative, seismic energy release graph is
shown in the following figure (220.127.116.11.2.11).
seismic energy release, determined for the 1 by 1 degrees
seismogenic region, for the period 1992 – 2006.
The accelerating deformation characteristics of the seismic energy
are initially revealed in this graph. This is more obvious in the
following figure (18.104.22.168.2.12) where was used the polynomial,
fitting (6th order) as determined on the cumulative seismic
energy release data.
The same accelerating deformation character of the seismicity of
the corresponding frame is indicated by the time gradient graph
(fig. 22.214.171.124.2.13) of the cumulative seismic energy release data.
126.96.36.199.2.13. Rate of
change in time, for the corresponding 1 by 1 degrees
seismogenic area, of the cumulative, seismic, energy
The comparison of the three cases presented (2 by 2, 1 by 1, and
frame oriented along fault strike) indicates that, the most
indicative case, for the detection of the accelerating deformation
status of a seismogenic region, is the one, where the assumed
regional seismogenic area (frame) conforms to the tectonic character
of the seismogenic area and the strike of the seismogenic fault.
The latter leads us to the following conclusion. In the wider
regional area, where a strong EQ will take place, the equation:
logR = 0.42M – 0.68
is valid and defines the strain, deformed, area.
However, in the vicinity of the fault, to be activated, the most
drastic, elastic deformations take place and this explains why the
EQ takes place in this fault, during the strain release.
Consequently, this leads directly to the Rebound Theory of Reid.
After having presented the above examples, what is proposed and
probably holds in the seismogenic area is the following:
- Up to a distance (R) from the epicenter location,
calculated by the equation (188.8.131.52),
the preparation of a strong EQ, strain is accumulated.
- Very close to the fault which will be activated, the
elastic deformation prevails and
the dominant model for the earthquake generation mechanism is the
one of the
It is obvious that, the accelerating deformation status of each
deep, lithospheric fracture zone must be studied in detail. Further
more, it is possible to compile maps which will indicate the
acceleration deformation spatial distribution all over a wide area,
and hence the corresponding, seismic risk will be evaluated.
In conclusion, the combined use of the lithospheric fracture zones
and the cumulative, seismic energy flow lithospheric model selected
for detailed analysis can provide us with valuable information for
the seismic energy charge of any area.
2.5.2. The lithospheric plate oscillating model.
The Earth-tides mechanism was early recognized as a potential
trigger for the occurrence of strong earthquakes. This approach was
followed to study the time of occurrence of EQs and their
correlation to Earth-tides. Knopoff (1964), Shlien (1972), Heaton
(1982) and Shirley (1988) suggested the Earth-tides as a triggering
mechanism of strong EQs, Yamazaki (1965, 1967), Rikitake et al.
(1967) studied the oscillatory behavior of strained rocks, due to
Earth-tides, while Ryabl et al. (1968), Mohler (1980) and Sounau et
al. (1982) correlated Earth-tides to local micro-earthquakes and
The following figure (184.108.40.206) demonstrates the mechanism which
generates the Earth-tides. On the left are shown the forces which
are applied on the Earth’s surface by the Moon or the Sun. Because
of these forces, the lithospheric plate of the Earth deforms, as it
is shown, on the right part of the figure.
applied on the Earth’s surface (Stacey, 1969) by the Moon or
the Sun (left) and the resulted, deformation (Garland, 1971)
of the lithospheric plate (right).
Earth-tides exhibit an oscillating mode of behavior. The study of
Earth-tide oscillations has shown that, basically, it consists of
the following main components, presented, in the table below.
----------- --- -----------------
M2 Principal Lunar
S2 Principal Solar
N2 Lunar Ellipticity
O1 Lunar Declination 25.82
P1 Solar Declination 24.07
M1 Moon declination 14 days
Ssa Moon declination 6
Finally, an Earth-tide wave is generated with in a year’s period,
because of the Earth’s motion in an ellipse with the Sun in the
The deformation of the lithosphere follows the oscillatory
character of the Earth-tides. Garland (1971), Stacey (1969), Sazhina
and Grushinsky (1971), study in detail this type of Earth’s
oscillation, not to mention the majority of the Geophysical
The lithosphere, following the Earth-tide oscillatory forces,
self-oscillates with maximum amplitude which varies as follows:
According to Garland (1971), the maximum elevation of the
lithosphere is of the order
According to Sazhina and Grushinsky (1971), the lithosphere
oscillates with a maximum
p – p value of 53.4cm for the
Moon component, while for the Sun component is only 24.6cm.
That means that, the maximum oscillatory amplitude can reach the
value of 78cm p – p, in cases, when the two components are added in
A sample of a day’s tidal oscillation is presented in the following
calculated for the 18th November, 1997, Greece ö = 380,
ë = 240.
The vertical scale is in mgals, while the horizontal one is in
minutes (1 day = 1440 minutes). In this figure it is possible to
identify the K1 (23.93hr), K2 (11.97hr) and S2 (12.00hr) components.
The next figure (220.127.116.11) represents the T = 14 days oscillation of
oscillation with T = 14 days, for the time period of
05/08/2000 – 13/10/2000
Figure (18.104.22.168) represents the yearly oscillation of the
Earth-tide, along with the 14 days period oscillation, superimposed,
22.214.171.124. Yearly and 14
day’s period, Earth-tide oscillation is shown, for the year
The next figure (126.96.36.199) represents the yearly and the 14 days
period oscillation of the Earth-tides.
188.8.131.52. Yearly and 14
day’s period Earth-tide oscillation calculated, for the
period 01/01/2001 – 31/12/2002.
Using the Rudman et al. (1977) methodology, tidal waves have been
generated. As long as the oscillating, Earth-tides forces act upon
the lithospheric plate, the latter is excited in the same mechanical
mode of oscillation.
According to this model, the lithosphere is considered as a rigid
plate which can oscillate, because of applied oscillating, external
forces. As a first approximated, simple model, is used the single
beam, supported at point A and B (fig. 184.108.40.206).
distribution (bottom graph) along a resting beam on two
supports A, B (Papaioannou, 1952).
Then the following equations hold:
p=c (constant weight load along its axis)
Há = Hâ = pl/2
ä÷(÷) = Dx = p(l/2 – x)
Mx = (px/2)(l-x) (bend moment)
Mì = max for Dì = 0
The load distribution of the beam (Papaioannou, 1952) along its
length is presented by the bottom part of the figure (220.127.116.11).
In the case of the lithospheric plate the maximum oscillating load
is applied, locally, by the Earth-tides.
Therefore, it is activated and deformed in an oscillatory mode.
This is presented in the following figure (18.104.22.168). The solid line
parallelogram represents the lithospheric plate at resting position,
while the dashed lines show the two extreme positions of its maximum
amplitude of oscillation.
The lower part of the figure shows the stress load distribution of
lithosphere, at the maximum or minimum of the Earth-tides.
plate at rest (solid line) and maximum amplitude oscillation
(dashed line). Bottom drawing indicates the stress-load,
along the plate, during its oscillation.
The oscillation of the lithospheric plate, due to Earth-tides, has
two severe consequences:
-The first one is that, at the peak amplitude of its
oscillation, when the plate is at nearly critical stress load, it
can reach the extreme conditions which are necessary for an
earthquake to occur (Thanassoulas et al. 2001), provided that, the
focal area itself, has been charged, enough, by the linear increase
of stress, due to plate’s motion. This takes place, especially, when
all the oscillating components of the Earth-tides are “in phase”.
-The second one is that, at the very same stress load
critical conditions, the focal area of an imminent earthquake
generates electrical signals (Thanassoulas 1991, Clint 1999), due to
its large crystal lattice deformation and to the piezoelectric
properties of the quartzite content in it.
An EQ will occur when the stress
load that increases gradually with time, due to plate’s motion, will
reach the failure stress level of the rock at the focal area. This
is represented, schematically, in the following figure (22.214.171.124).
126.96.36.199. The EQ occurs
(red arrow) when the stress-load (in arbitral units) exceeds
its critical level (red line).
The Y-axis represents stress increase in arbitrary units, while the
x-axis represents time in arbitrary units. The horizontal red line
represents the rupture stress level of the rock, at the focal area,
and the blue, inclined, line denotes the state of stress charge of
the rock of the focal area. The red arrow indicates the occurrence
of the earthquake when the stress load of the rock equals the
critical rupture stress level.
The following (fig. 188.8.131.52) stress charge load, increase form,
takes place in the seismogenic region, by combining the two
different stress increase mechanisms, the linear caused, by the
plate’s motion and the oscillatory one, caused, by the Earth-tides.
184.108.40.206. The Linear –
oscillatory mode increase of stress is presented in sketch
As long as the repetition period of strong EQs (red arrows) is of
the order of some years, it is obvious that critical stress levels,
for the occurrence of an EQ, will occur when daily, local, maximum
peaks of the stress-load curve, reach the fracturing stress level.
In the case of a 30 years recharge period of a seismogenic area, the
ratio of the recharge period to the oscillations one, is:
T recharge / 1 day = 30*360/1 = 10800.
That value of 10800 satisfies the criterion of:
T of plate recharge >> T of oscillation
Therefore, theoretically, strong EQs are expected to occur only at
the maximum peaks of the Earth-tide waves, since these will,
firstly, reach the fracturing stress level of the seismogenic area
The practical use of the lithospheric oscillating plate model will
be presented at the chapter that analyzes the method of the
determination of the time of occurrence of a strong earthquake.
2.5.3. The homogeneous ground Earth model.
The upper part of the Earth (lithosphere – crust) consists of
different lithological - geological formations which have been
affected by the different tectonic processes, which took place in
the past. This geological-tectonic complex, although it affects in a
very general approach, the Earth’s physical properties spatial
distribution, it can be considered as non-existent under certain
physical assumptions. In other words Earth is assumed as a
homogeneous ground. That means that its physical properties do not
depend on the coordinates (x, y, z) which describe the Earth
At a first glance, if we take into account the real geological,
tectonic conditions which exist in a seismogenic region, this
assumption appears to be false. However, in the specific case of the
generation of an electric field in the focal area and, especially,
when its wavelength is much more larger than the wavelength of the
geological-tectonic elements, met, in the seismogenic area, then,
for this specific case only, the seismogenic area can be considered
as consisted of homogeneous ground, as far as it concerns its
In physical terms, the X-ray spectroscopy is a similar example, to
what was mentioned before. This methodology uses X-rays because
their wavelength is comparable, in size, to the crystal lattice grid
size and therefore the diffraction of the X-rays is utilized. On the
other hand, the crystal lattice is non-existent for the incident
large wavelength “light” and therefore, in general, it passes
through it without any other consequence.
The model of the homogeneous ground, in terms of its electrical
properties, will be used in the analysis of the precursory,
electrical signals, which will be presented later on. The basic,
theoretical analysis of an electric field that is generated by some
mechanism in a homogeneous material (in this case in the upper part
of the Earth, Telford et al. 1976, Grant and West 1965, Bhattacharya
and Patra 1968) is as follows:
It is assumed that a current (I) flows in an isotropic, homogeneous
medium. If (J) is the current density and (äA) is an element of
surface, then the current (I) that passes through äA is: I = J* äA.
The electric field (E) and the current density (J) are related
through the Ohm’s law:
J = óÅ
Where (ó) is the conductivity of the medium.
The electric field (E) is the gradient of a scalar potential,
Since the current flows outwards radially, in all
directions, from the point current source, the generated
equipotentials which are orthogonal to the current flow lines, will
be spherical surfaces, given by r = constant. This is presented in
the following figure (220.127.116.11).
surfaces generated, by a point current source in an
isotropic, homogeneous medium (Telford et al. 1976).
When the generated, equipotential, spherical surfaces are
intersected by a plane of conductivity discontinuity, such as the
ground surface in the case of the solid Earth, then equipotential
circles are formed at the plane of conductivity discontinuity. This
is illustrated in next figure (18.104.22.168).
22.214.171.124. Plan view of
equipotential circles, generated on the ground surface, by
the intersection of the spherical, equipotential surfaces,
generated by the Point Current Source (PCS) in the ground,
and the ground surface plane.
The latter is the main assumption, which is used throughout the
methodology, and is applied, for the epicentral location
determination of an imminent, strong earthquake.