Spatial distribution of site-effects and wave propagation properties in Thessaloniki (N. Greece) using a 3D finite difference method

Scientists from the Geophysical Laboratory (Department of Geophysics, School of Geology of the Aristotle Univ. of Thessaloniki) have studied the site effects of seismic motion in the metropolitan area of the city of Thessaloniki (Northern Greece) for various seismic earthquake scenarios with a 3D finite-difference modeling approach, using the HellasGrid Infrastructure and the EGI with the support of the Scientific Computing Center at A.U.Th.

The city of Thessaloniki (Northern Greece) was selected since it is located in a moderate-to-high seismicity region (Papazachos et al., 1983), with the Servomacedonian massif and Northern Aegean through areas exhibiting the highest seismicity (Figure 1). The city suffered several large earthquakes throughout its history, many of them causing significant damages and human losses (Papazachos and Papazachou, 2002).

Thessaloniki Earthquake Map

Figure 1: Map of known earthquakes with M≥3.0 which occurred in the broader area of central–northern Greece from the historical times (550 BC) till 2007 (Figure after Skarlatoudis et al., 2011a).


An explicit 3D 4th-order velocity-stress finite-difference scheme with discontinuous spatial grid was used to produce synthetic waveforms with numerical simulations. The scheme solves the equation of motion and Hooke’s law for viscoelastic medium with rheology described by the generalized Maxwell body model. Details on the scheme, its grid and material parameterization are provided by Moczo et al. (2002), Kristek & Moczo (2003), Moczo & Kristek (2005), Moczo et al. (2007) and Kristek et al. (2009b).

The computational model used for the simulations is based on the geophysical-geotechnical model and the dynamic characteristics of the soil formations proposed by Anastasiadis et al. (2001) and covers an area of 22 x 16 Km2 (dotted rectangle in Figure 1) (Skarlatoudis et al., 2007; 2008b; Skarlatoudis et al., 2010).

Numerical simulations were performed for six seismic scenarios, corresponding to three different hypocentral locations and two different focal mechanisms for each one. Seismic scenarios with E-W trending normal faults are referred as scenarios (a), while the ones with NW-SE trending normal faults as scenarios (b) (Figure 2). Both types of normal faults (E-W and NW-SE) are the dominant types of faults in the vicinity of the broader Thessaloniki area (e.g. Vamvakaris et al., 2006). Synthetic waveforms were produced for a coarse grid of receivers, in order to study the spatial variation of site-effects on seismic motion in the broader metropolitan area of Thessaloniki (Figure 2).

Earthquake Simulation Scenarios          

Figure 2: Earthquake locations used for the examined seismic scenarios (red stars) and the focal mechanisms used for each scenario. The coarser grid of receivers used for studying the spatial variation of various waveform and site-effect parameters for the six earthquake scenarios is also shown (black diamonds). The location of site OBS, used as a reference station in computations, is denoted with a yellow triangle (Figure after Skarlatoudis et al., 2011a).


The application that implements the 3DFD method is using the MPI libraries for inter-process communications, namely the mpich2 implementation. The compilation and execution of the code was tested in different types of machines and different Fortran90 compilers (commercial and free). The most accurate results and the minimum execution time in each system were achieved with the use of the commercial compiler Pathscale (version 3.0) (Skarlatoudis et al., 2008a). The execution of the 3DFD code is demanding in terms both of CPU power and computer memory. For the aforementioned computational model the memory demands reached 20 GB and the time of computation (per model) was approximately 15 on the HellasGrid Infrastructure with the synchronous usage of 40 Intel Xeon processors.

The implemented workflow relies mainly on gLite middleware (Figure 3). Also a large number of test runs for checking the compatibility of the results on the Grid with the ones obtained from other computational infrastructures have been performed. Moreover the scaling of the execution of the code on the HellasGrid Infrastructure was examined (Skarlatoudis et al., 2008a).

3D FDTD Application Workflow          

Figure 3: Schematic representation of the workflow in HellasGrid infrastructure (Figure after Skarlatoudis et al., 2008a)


Various measures, estimated from the 3D synthetic waveforms that can provide a more detailed evaluation of site-effects, such as spectral ratios, Peak Ground Velocity (PGV), cumulative kinetic energy and Housner Intensity, were used to probe the site-effects spatial distribution and ground motion variability. In Figure 4 the Peak Ground Velocity (PGV) ratio is shown for the 3D over the corresponding 1D bedrock reference model [(PGV3D)/(PGV1D)], estimated for the coarser grid of receivers and for the two horizontal components of ground motion, for all scenarios studied (Skarlatoudis et al. 2011a). The observed relative PGV distribution from the six scenarios, exhibits high values along the coastal zone, with the highest value (~4) shown in the area near the city harbor for the E-W component. High values of relative PGV are also observed in the western parts of the model for the E-W component.


Figure 4: Spatial variation of the average, from the six seismic scenarios, ratio [(PGV3D/PGV1D)], for the horizontal components of ground motion (Figure after Skarlatoudis et al., 2011a)


The 3D wave propagation characteristics of the 4th July, 1978 aftershock (M5.1) of the 20th June, 1978 strong mainshock (M6.5) that struck the city of Thessaloniki were also studied using the 3D finite-difference approach. In Figure 5 the spatial distribution of damages in the metropolitan area of Thessaloniki after the 1978 mainshock is presented (left figure) (Leventakis, 2003), together with the corresponding distribution of the RotD50 ground motion measure of the (PGV3D)/(PGV1D) ratio, for the frequency band 0.2Hz-3Hz (Skarlatoudis et al., 2011b). According to Leventakis, (2003) the largest damage was recorded in the city harbor area and parts of the eastern area of the Thessaloniki. Despite the various limitations of the comparison, a quite good correlation is observed between the damage distribution and the PGV spatial variation, suggesting that the role of local site amplifications studies here is much more important than other factors (e.g. differences in source radiation pattern, non-linearity, etc.).

Thessaloniki Damage Distribution

Figure 5: (Left) Spatial distribution of damage distribution in Thessaloniki caused by the mainshock of July 1978 according to Leventakis (2003). (Right) Spatial distribution of the RotD50 measure of relative PGV values (amplifications) from filtered (0.2Hz-3Hz) horizontal components (Figure after Skarlatoudis et al., 2011b).


This work has been partly performed in the framework of PENED-2003 (measure 8.3, action 8.3.4 of the 3rd EU Support Programme) and the Greek-Slovak Cooperation Agreement (EPAN 2004-2006). Most of the computations were realized on the EGI and HellasGrid infrastructureσ with the support of the Scientific Computing Center at the Aristotle University of Thessaloniki (AUTH). A significant part of the results presented here have been published in peer-review journals (see inline references) and/or presented in national and international conferences (see references at the end of this document).


Contact details:

  • Papazachos C.B., Professor, AUTH, kpapaza (at)
  • Skarlatoudis A.A, Dr. Seismologist, AUTH, askarlat (at)
  • Scientific Computing Center, AUTH, contact (at)



  1. Papazachos, B. C., Tsapanos, T. M. and Panagiotopoulos, D., (1983). The time, magnitude and space distribution of the 1978 Thessaloniki seismic sequence. The Thessaloniki northern Greece earthquake of June 20, 1978 and its seismic sequence. Technical chamber of Greece, section of central Macedonia, 117-131, 1983.
  2. Skarlatoudis A.A., C.B. Papazachos, P. Moczo, J. Kristek, N. Theodoulidis and P. Apostolidis, (2007). Evaluation of ground motions simulations for the city of Thessaloniki, Greece using the FD method: the role of site effects and focal mechanism at short epicentral distances, European Geosciences Union (EGU) General Assembly, Vienna, Austria.
  3. Skarlatoudis A.A., Korosoglou, P., Kanellopoulos, C. and Papazachos C.B, (2008a). Interaction of a 3D finite-difference application for computing synthetic waveforms with the HellasGrid infrastructure, 1st HellasGrid User Forum, Athens, Greece, 3rd EGEE User Forum, Clermont-Ferrand, France.
  4. Skarlatoudis A.A., C.B. Papazachos, P. Moczo, J. Kristek and N. Theodoulidis, (2008b). Ground motions simulations for the city of Thessaloniki, Greece, using a 3-D Finite-Difference wave propagation method, European Geosciences Union (EGU) General Assembly, Vienna, Austria and 31st European General Assembly of the European Seismological Commission, Chania, Greece.
  5. Skarlatoudis A.A., C.B. Papazachos and N. Theodoulidis, (2011b). Site response study of the city of Thessaloniki (N. Greece), for the 04/07/1978 (M5.1) aftershock, using a 3D Finite-Difference wave propagation method, accepted for publication in Bull. Seism. Soc. Am.