1 Introduction

Selinunte, one of the first Greek colonies located in south-west Sicily, is known for its large temple ruins. It was founded in the seventh century B.C.E. and expanded in the following two centuries (Reinganum 1827). After occupation by the Carthaginians in 409 B.C.E., periods of settlement—as well as abandonment of the city—have been documented (e.g., Guidoboni et al. 2002). Today, the remains of Selinunte form the largest archeological park in Europe. Two rivers divide the landscape morphologically in three parts: a western part (or western hill), the Acropolis, and an eastern part (or eastern hill) (Fig. 1).

Fig. 1
figure 1

Digital terrain model (DTM) of Selinunte and the surrounding area. Overlying is the geological map based on Amadori et al. (1992). The black rectangles indicate three measurement sites next to Triolo Temple in the western part (W1), on the Acropolis (A1), and north and south of Temple G in the eastern part (E1 and E2) of Selinunte Archaeological Park (shown in detail in Figs. 2, 3, and 4). Locations of earlier seismic profiles measured by Rabbel et al. (2014) are shown by blue points crossing the river valleys. Coordinates are in UTM (33S); the arrow in the insert in the upper left points to the study area in western Sicily

The excavation of the archeological site started in the nineteenth century and is still ongoing (e.g., Montana et al. 2018). The destruction of the temples due to earthquakes was first assumed by Hulot and Fougères (1910). Archaeoseismological researches have more recently confirmed the earthquake hypothesis (Guidoboni et al. 2002; Bottari et al. 2009). They identified two earthquakes, the first one that occurred between 370 and 300 B.C.E. and seems causative for the destruction of the temples in the western part of Selinunte and probably also Temple R located on the Acropolis (Marconi personal communication). The second earthquake between 330 and 500 C.E. destroyed the temples on the Acropolis and the eastern side (Guidoboni et al. 2002; Bottari et al. 2009). However, it is still debated whether all temples of Selinunte show evidence of earthquake damage (Bottari et al. 2009). Although there is no definite information about the conditions of the temples when the earthquakes occurred, according to Guidoboni et al. (2002) and Bottari et al. (2009), there is no doubt that two damaging earthquakes hit the site. Historical and instrumental seismicity of the southwestern area of Sicily is low (e.g., less than 20 earthquakes from M 4.0 up to M 6.4 recorded in the last four centuries according to Italian seismic catalogs; Boschi et al. 2000; Guidoboni et al. 2007; Rovida et al. 2016). The strongest seismic event recorded in western Sicily was the 1968 Belìce earthquake sequence including several shocks with M~ 5, the strongest of which was M 6.4 (Rovida et al. 2016) at an epicentral distance of 23 km from Selinunte. Before the 1968 Belìce seismic sequence, this portion of the island was considered as seismically quiescent region. Focal solutions provide possible faulting mechanisms that range from thrusting on a WSW–ENE striking plane to right lateral transpression on a NNW–SSE striking plane (Bottari 1973; Gasparini et al. 1982; Anderson and Jackson 1987; Frepoli and Amato 2000). However, in the light of the earthquake damage in Selinunte, more information about the ancient earthquakes could advance the seismic hazard model of SW Sicily. Quantitative archaeoseismological studies (e.g., Hinzen 2005; Hinzen et al. 2013; Hinzen et al. 2018) have shown that damage observed in archeological remains can help to constrain the strength of ancient earthquakes. An important step in the modeling of ancient ground motions is the determination of seismic site effects at the time of the earthquake. Hinzen and Weiner (2009), Hinojosa-Prieto and Hinzen (2015), and Hinzen et al. (2016) have given examples of site effect evaluation in archaeoseismological context. For further investigations on the potential sources of the two earthquakes that affected Selinunte (Guidoboni et al. 2002; Bottari et al. 2009) the determination of site effects, which is the purpose if this paper, is an important step.

As the geologic stratigraphy of Selinunte supports the potential existence of low-velocity layers, we used classical refraction seismic combined with multichannel analysis of surface waves (MASW) (Park et al. 1999; Xia et al. 1999; Miller et al. 1999) and ambient noise measurements applying the wireless array analysis (WARAN) system (Ohrnberger et al. 2006; Hannemann et al. 2014). The frequency wavenumber (f-k) analysis (Robinson 1969; Neidell and Taner 1971; Lacoss et al. 1969; Foti et al. 2000; Ohrnberger et al. 2004b) was used to determine dispersion curves for four measuring sites in Selinunte: one at the western part (W1), one at the Acropolis (A1), and two at the eastern part (E1, E2) (Fig. 1). Suitable shear wave velocity, vS, profiles were derived by inversion of the dispersion curves including fundamental and higher modes. The results were controlled by forward calculations of synthetic seismograms which were compared to the measured MASW profiles. Furthermore, horizontal-to-vertical spectral ratios (HVSRs) of ambient noise from all stations were calculated and compared to the inversion results. The study was done in the framework of the SELINUS (Seismological Evaluation of Lasting Induced Natural Undoing of Sanctuaries) project, a cooperation between Cologne University, Instituto Nationale di Geofisica e Vulcanologia (INGV), the Università di Napoli “Federico II,” and the University of Pennsylvania.

2 Geologic setting of the site

Figure 1 shows the area of Selinunte situated on Pleistocene sedimentary deposits and terraced calcarenites (Barreca et al. 2014). The typical stratigraphic section consists of a near-surface calcarenite layer with an underlying sand layer followed by a thick clay deposit (Ruggieri and Sprovieri 1977).

The three parts of the archeological park Selinunte are morphologically separated by the valleys of the Cottone and the Modione Rivers (Fig. 1). In the 1990s, several geophysical investigations such as reflection seismic surveys and geoelectric measurements (Balia 1992; Amadori et al. 1992) were carried out mostly in the eastern part of Selinunte (Brizzolari et al. 1992). Amadori et al. (1992) concluded that here the maximum thickness of the calcarenite layer is 6–7 m below the temples, whereas a thickness of 3 to 5 m was postulated for the Acropolis. No previous geophysical investigations were made at the western part. In a more recent study, Rabbel et al. (2014) measured two reflection seismic profiles (Fig. 1) across the river valleys motivated by the search for the position of the ancient harbor.

According to those previous investigations and the geological map (Fig. 1), the deposits of the Calcarenite di Marsala Formation (Ruggieri and Sprovieri 1977) is expected to be present at least at the eastern site and at the Acropolis. Based on the lower altitude of the temples at the western site it is questionable whether a high velocity layer exists.

3 Dispersion curves

The dispersive behavior of surface waves has gained importance in the past few decades in the seismic site studies, particularly for the determination of vS profiles. An established method to measure dispersion curves is the MASW (Park et al. 1999; Xia et al. 1999; Miller et al. 1999), an advancement of the spectral analysis of surface wave (SASW) method (Nazarian and Stokoe 1984). This active seismic method works best at locations with homogenous layers, a fact which can become problematic in case of long measuring profiles. Bignardi et al. (2014) showed that the dispersion curve is not sensible for gentle slopes with less than 15° inclination while layer dipping beyond 30° will cause misleading results. Lateral heterogeneities within the layers can lead to doubtful results anyway (e.g., Lin and Lin 2007; Bignardi et al. 2014).

Previous investigations (e.g., Forbriger 2003a) showed that low-velocity layers produced an inverse trend of the Rayleigh wave phase velocity of the dispersion curves, which can be formed by multiple modes. This phenomenon and how to handle its interpretation has been part of several investigations, e.g., Forbriger (2003a), Maraschini and Foti (2010), and Pan et al. (2013).

In this paper, the dispersion curve for each site is divided into sections assigned to certain modes. For the multiple mode inversion, the program DINVER (Wathelet et al. 2004; Wathelet 2008) is used. It is included in the GEOPSY software package and solves inversion problems by the application of an improved neighboring algorithm (Wathelet et al. 2004; Wathelet 2008). The resulting velocity profiles were rechecked by forward modeling of a MASW section using Qseis software (Wang 1999) which calculates the Green’s functions for the previously derived velocity models.

4 Data acquisition and processing

In the following, the active and passive seismic experiments carried out in Selinunte Archaeological Park as well as the data processing steps are described.

4.1 Active seismic experiments

Four seismic refraction profiles and MASW surveys were measured in spring 2017 in the Selinunte Archaeological Park (Figs. 2, 3, and 4). A 24-channel ABEM Terraloc Mk6 v2 was used to collect the signals of the 4.5-Hz vertical component geophones. The spread was 69 m with an interstation distance of 3 m. Three shot points were used for the refraction survey: at 0 m with the first geophone at 1.5 m, in the middle of the spread length at 34.5 m, and at the end at 69 m; here, the last geophone was positioned at 67.5 m (Fig. 5a). For the MASW spread, an equidistant receiver spacing of 3 m was used with a source offset of about 18 m (e.g., http://www.masw.com last accessed 05.2018) for each side of the spread as shown in Fig. 5a. A 5-kg sledge hammer and a polypropylene hammer-plate were used to generate P-waves. For the measurement, 16,384 samples at a sampling interval of 100 μs with a pre-trigger time of 10 ms were selected. Depending on noise conditions, at least five shots per source position were stacked. All shot and receiver positions were surveyed by differential global positioning system (DGPS).

Fig. 2
figure 2

a Aerial view (GoogleEarth, image, 07.2016, last accessed 02.2018) of the measuring area at the eastern part of Selinunte (Fig. 1). Three major temples are labeled E, F, and G. Overlain are the measuring sites E1 north of Temple G and site E2 north of Temple F. Geophone and shot locations of seismic refraction and MASW profiles are shown by yellow triangles; seismometer positions of the WARAN array measurements are shown as red dots. The dashed blue lines show profiles of earlier seismic reflection measurements (Balia 1992). Coordinates are in UTM (33S). b View from Temple G towards north with measuring site E1 in an olive grove. The refraction line was set up parallel to the row of trees. c View from Temple G towards south. In the foreground is the measuring site E2; behind it are the ruins of Temple F and the anastylosis of Temple E, Mediterranean Sea in the background. The refraction line was set up across the grassland in the foreground (photos K.-G. Hinzen)

Fig. 3
figure 3

a Aerial view (GoogleEarth, image, 07.2016, last accessed 02.2018) of the measuring area at the Acropolis of Selinunte denoted A1 (Fig. 1). Geophone and shot point positions of the seismic refraction and MASW profile are sown by yellow triangles. The red circles show the layout of the seismic array measurement. The blue circles give locations of additional HVSR measurements at Temple A. Coordinates are in UTM (33S). b View from West to East with a small aperture ambient noise WARAN array in the foreground. In the back a row of reconstructed columns of Temple C is visible. The refraction line was set up on the east-west-trending path (photo K.-G. Hinzen)

Fig. 4
figure 4

a Aerial view (GoogleEarth, image, 07.2016, last accessed 02.2018) of the measuring area at the western part of Selinunte denoted W1 (Fig. 1). The geophone and shot point positions of the seismic refraction profile and MASW measurement are indicated by yellow triangles. The red circles show the layout of the seismic array measurement. Blue circles show the location of additional HVSR measurements. Coordinates are in UTM (33S). b View towards North with the Triolo Temple in the foreground. The refraction line is located along the road with the reverse shot point close to the car (photo K.-G. Hinzen)

Fig. 5
figure 5

a Layout for the seismic refraction and MASW survey including geophone and shot point positions. b Seismogram examples from site E1 for the three shot points used for the refraction analysis

4.2 Ambient noise measurements

For the ambient noise measurements, the WARAN system of Potsdam University was used which has been described in detail by Ohrnberger et al. (2006). It consists of eight LE-3D/5s seismometers, Earth Data PS2400 digitizers (http://seismo.geo.uni-potsdam.de/waranproj last accessed 05.2018), and wireless network communications between all stations and the recording computer. All array measurements were sampled with a frequency of 100 Hz. Stations were positioned in circles of increasing radius while one station remained at the center. Thus, the requirement for using spatial auto-correlation (SPAC) is satisfied as well as a sufficiently coverage for modified SPAC is warranted (Aki 1957; Bettig et al. 2001). Recording time was at least 30 min for each circular setup. Depending on the site situation (ruins, vegetation, etc.) three to four circles or elliptical station patterns were realized. Table 1 summarizes the radii, specifications, and the measuring times. In cases where a range is given for the radius, only five of the eight stations were moved. The setup of an array with small radius at site W1 is shown in Fig. 3b. Seismometers were covered by buckets to reduce wind generated noise and to protect the instruments from direct sun. The recorded ambient noise data was also used to analyze HVSRs.

Table 1 Array size, recording duration, and parameter setting of the f-k analysis for the ambient noise measurements at Selinunte using WARAN

4.3 General data processing

The general processing workflow shown in Fig. 6 consisted of seven steps which were followed for each site: (1) analysis of the refraction seismic survey; (2) f-k analysis of the MASW and ambient noise (WARAN) dataset; (3) SPAC of the array measurements; (4) dispersion curve forward calculation with gpdc, a routine included in the GEOPSY software package (http://www.geopsy.org last accessed 08.2017) for a first guess of the subsurface profile; (5) inversion of the dispersion curve; (6) forward calculation of synthetic seismograms of a MASW survey for comparison; and (7) comparison of the forward calculation results from (6) with inversion results from (5) and refinement of inversion results if necessary. According to the stratigraphic section of Selinunte (Ruggieri and Sprovieri 1977), a low-velocity sand layer is expected beneath a shallow calcarenite layer. To avoid misinterpretation due to hidden layers, the refraction seismic measurements (step (1)) were only used to estimate the direct P-wave velocity, vP, of the top layer (recent deposit) as well as the velocity of the first refracted wave supposedly from the interface between recent deposits and the calcarenite layer and the layer thickness. We also tried a full tomographic inversion, but ray coverage at the crucial depths was insufficient.

Fig. 6
figure 6

Schematic workflow for the combined use of refraction seismic, MASW, and ambient noise array measurements (WARAN). Numbers refer to the steps described in the text

In step (2), dispersion curves from the MASW and ambient noise measurements were determined. The MASW data f-k analysis was made with the active seismic experiment tool implemented in GEOPSY. For all sites, a time window of 1.4 s was used. To consider geometrical spreading, a normalization factor of the inverse of the root of the distance was applied to each trace. The bandwidth bw was set to 0.1. By multiplying the considered frequency with the factor (1 ± bw), the limits of each frequency band were determined. The conventional f-k toolbox provided by GEOPSY was used to analyze the WARAN data. The length of the analysis time windows for each array is given in Table 1; short-term disturbances (e.g., from passing vehicles, etc.) were excluded. A 5% overlap of the analysis time windows was chosen. For each frequency node of the dispersion curve, the corresponding period was multiplied by a chosen factor ftw between 20 and 50, depending on the array size, to set the length of the time window; ftw was determined empirically based on tests of resulting resolutions of the dispersion curve (for selected values, see Table 1). The f-k gridding parameters were automatically calculated depending on the size of the array by the program warangps, included in the GEOPSY software package. The output of the f-k analysis for array measurements was used as input for the program max2curve included in the GEOPSY software package (Di Giulio et al. 2006; Wathelet et al. 2008) to calculate the probability density function (pdf) for a range of slowness values. This was done separately for each radius of the circular arrays as shown in Fig. 7 and in the Appendix. The black solid and dashed lines indicate the limitations for the measured dispersion curves due to resolution depending on the array aperture and due to aliasing effects, respectively (e.g., Haubrich 1968; Harjes and Henger 1973). The dispersion curves were picked from the pdf distributions and combined with the dispersion curve derived from the MASW measurements to extend the frequency range to higher values. Due to the inverse trend of at least two dispersion curves (E2 and W1) without clear discriminability of higher modes SPAC was applied in step (3) to test for potential mode transition. As indication, we looked for points of discontinuity within the autocorrelation ratio curves individual of each analyzed ring following Ohrnberger et al. (2004a). In general, with SPAC, often an effective mode influenced by different modes of the dispersion curve can be observed for the coefficients and phase velocities (e.g., Ikeda et al. 2012). However, due to the parallel analysis of multiple radii an indication of the presence of higher mode is possible if a deviation from a ‘smooth’ curve is observed in the same frequency band for different array radii. In Fig. 8 an example is given for the 10 m radius measurement of the W1 site. The red circles mark visible deviation points in the curves. At such frequencies the dispersion curve of the site was divided into sections and used as guide line in the forward calculation (step (4)) with the program gpdc, part of the GEOPSY software package (http://www.geopsy.org last accessed 08.2017), for fitting the fundamental mode as well as higher modes to the dispersion curve. The aim of this step was to get a first guess of suitable vS values and layer thicknesses as a priori information for the inversion process in step (5).

Fig. 7
figure 7

Example of a color coded probability density functions (pdf) for the ambient noise array measurements with the WARAN at site E1 for radii of the circular arrays of 7, 20, and 60 m (top to bottom). Solid and dashed thin black lines indicate the resolution limits due to the array aperture and aliasing effects, respectively. The segments of the dispersion curve used in the inversion process are displayed as thick black lines

Fig. 8
figure 8

Example of autocorrelation ratio curves derived for the 10 m radius of the ambient noise measurements at site W1. Deviation points are marked by red circles

The multimode option of DINVER was used for the inversion calculations. The dispersion curve of each site was divided into sections, with breaks at the frequencies derived in step (3) testing different mode indices for each part. Since the dispersion curve inversion is rather insensitive towards vP and density changes (e.g., Forbriger 2003b), these parameters were fixed during the inversion process. The density was set to 2000 kg/m3. In addition, a range between 0.2 and 0.5 was chosen for the Poisson’s ratio for all inversion calculations. The parameter range of vS and the bottom depth, Db, of the layers were chosen based on the results of the forward calculation in step (4) and afterwards optimized based on inversion results. The parameter ranges used for the final inversion process are listed in Table 2 for each site. The P- and S-wave velocities were connected to each other and thus forced to have common layer boundaries. As an additional constraint for each site, only the third layer was allowed to have a lower velocity than the overlying stratum. For the final vS profiles, five inversion runs using the same parameter space were considered. The iteration process was continued until the misfit of the inversion was stable. For site A1 and E1 10,200 models, for site E2 and W1 15,300 were tested. The best-fitting model of each run was used for an averaged model considered as final inversion result.

Table 2 Parameter ranges used for the inversion of combined dispersion curves for vS and Db; HS-halfspace

In step (6), the resulting velocity profiles were used to calculate synthetic seismograms of a MASW survey with full Green’s functions using the code Qseis (Wang 1999) (Fig. 6). The last step (7) was the comparison of the synthetic MASW and the measured dispersion curves to check the reliability of the velocity profile. In case of an unsatisfactory fit, the velocity profiles, derived by the inversion, were adjusted particularly for the top layer with low resolution of vS and the layer thicknesses; an example is shown in Fig. 9.

Fig. 9
figure 9

Example of data from step (7) of the processing for site E2. a Dispersion curve derived from the MASW profile measurement. b Dispersion relation of synthetic MASW seismograms calculated with the velocity model from the inversion (processing step (6)); same as in c after refinement of the top layer’s velocity

5 Results and interpretation

In the following, the data, results, and interpretation of the four locations at three sites (Fig. 1) will be given.

5.1 Site E1

The site E1 is located in the north of Temple G on the eastern hill (Fig. 2) mainly situated within a grove of olive trees. Due to this recent agricultural use of the terrain, the top soil might have been altered by aerating. Small bushes and high-growing grass were present, particularly dense at the northern half of the profile; numerous medium to large cobble stones covered the ground.

5.1.1 Seismic refraction survey

The seismograms for the three refraction shot positions at site E1 are shown in Fig. 5b and the travel time curves in Fig. 10. First breaks of the direct and the first refracted wave are clear. A second refracted wave can be observed in the data but does not overtake the first one within the spread. Since the presence of low-velocity layers can lead to misinterpretation, the depth of the second refractor with a vP of 2330 m/s cannot be determined. The P-wave velocity of the direct wave is vP0 = 655 m/s and of the first refracted wave vP1 = 1525 m/s. The dip of the refractor is only 0.5° with a depth of 3.6 to 4.1 m at the shot points.

Fig. 10
figure 10

Travel time curves of four refraction seismic profiles; the sites in the Selinunte Archaeological Park are indicated in each graph (see Fig. 1 for location). First breaks are shown as crosses in red for shot-positions at 0 m, reverse shots at 69 m in blue, and green for the shot located in the center of the spread. The travel time curves, which correspond in color to the first break symbols, were adjusted interactively to the first breaks considering a single inclined layer over a halfspace model. For the reverse shot at site W1, the first breaks of the refracted wave fade into the noise at about half the profile distance; the light blue crosses mark later arrivals in the seismogram

A striking feature in the seismograms is a change in the velocity of the ground roll (Fig. 5b). According to investigations by Lin and Lin (2007), this might be an indication of a lateral heterogeneity like a fault. Since the first arrivals of direct and refracted waves are not affected, the lateral change is expected to be below the two top layers. Reconstruction plans of Hulot and Fougères (1910) show evidence of former buildings north of Temple G. The position of these constructions correlates well with the changes in the seismic profile. It is possible that the change in the velocity of the surface waves is an effect of ancient building foundations. However, the publications and investigations from 1910 focused on the temples and did not mention explicitly whether and where they found remains of buildings north of Temple G.

5.1.2 MASW

The dispersion curves derived by f-k analysis (Fig. 11) for the two shot positions at − 18 m and 87 m show differences at frequencies above 15 Hz. This is an indication of a heterogeneous shallow subsurface structure and thus agrees with the results from the seismic refraction survey. Therefore, only the section of the dispersion curve at frequencies below 15 Hz was used in the inversion calculations.

Fig. 11
figure 11

Dispersion curves of the MASW measurements. For each of the four sites, the curves from the shot position at − 18 m and the reverse shot at + 87 m are shown in the left and right columns, respectively. The labels in the upper left corner of the plots indicate the location of the measurements

5.1.3 WARAN

Before inversion, the dispersion curves from the measurements with the three array apertures were combined (Fig. 7). The curve based on the largest aperture (60 m) shows two modes that can be visually separated. The slowness increases with frequency and thus does not directly indicate the presence of the expected low-velocity layer.

5.1.4 Inversion results

For the inversion process, it was assumed the fundamental mode ranges from 1 to 1.4 Hz, and values at higher frequencies were associated with the first higher mode. The result of the inversion is shown in Fig. 12. An uncertainty of about 10% of the slowness values has to be considered when picking the dispersion curves. Therefore, in a frequency range of 7 to 23 Hz, individual modes cannot be distinguished. Test calculations showed no significant difference in the velocity models when a second higher mode was assumed to exist in this frequency range; however, misfit values increased in general. The best-fitting model is given in Table 3. It includes a thin layer of higher velocity which can be interpreted as calcarenite and which correlates well with the result of the refraction experiment.

Fig. 12
figure 12

Results of the inversion of the dispersion curves from the combined WARAN and MASW measurements for the four sites, which are indicated in the upper left corner of the graphs. In the left column, the black bars indicate the parts of the dispersion curves used in the inversion. The color of the calculated dispersion curves indicates the misfit of calculated and measured values as quantified by the color bar on the right. The solid black lines show the modes of the dispersion curve used to estimate vS velocity models which are shown in the right column using the same color scheme

Table 3 Resulting layer thickness, h, and vS, vP and materials for four sites at the Selinunte Archaeological Park (* estimated from vs) 

5.2 Site E2

The site E2 is located south of Temple G and north of Temple F (Fig. 2). The terrain is flat with a solid surface. In the center, a frequently used path crosses the site (Fig. 2c), which due to vehicle traffic is probably more compacted than the surroundings. For the refraction and MASW profile, a path was cut through medium to high dry grass.

5.2.1 Seismic refraction survey

The picked first breaks and travel time curves are shown in Fig. 10. The analysis results in a vP0 = 750 m/s for the top layer and vP1 = 1335 m/s for the refracted wave. The thickness of the top layer is 1.1 m at the shot point in the east and 5.9 m at the reverse shot in the west corresponding to a 4.8° dip towards west.

5.2.2 MASW

The dispersion curves of the two shot points at − 18 m and 87 m are in good agreement with up to a frequency of 60 Hz (Fig. 11) even though the refraction indicated inclined layering. However, Bignardi et al. (2014) could show that in case of a dip smaller than 15°, no significant effects were observed in the dispersion curves.

5.2.3 WARAN

The dispersion curve from the WARAN measurement at site E2 (Fig. 12) is combined from the results of arrays with 3, 10, 22, and 37 m aperture. The data from the smallest array are outside the frequency and slowness limits imposed by the aperture size and not further used in the inversion. The general trend of the dispersion curve matches the one derived by MASW measurements for site E2 supporting the idea of inverting the combination of both.

5.2.4 Inversion results

The inverse trend of the dispersion curve indicates a mode transition caused by the presence of a low-velocity layer. Results of SPAC were used to identify frequencies for probable mode transitions. Frequencies up to 40 Hz were used in the inversion, and results are shown in Fig. 12. The interpretation of modes at higher frequencies seemed not appropriate since the contribution of the individual sections to the dispersion curve would become too small. However, the overall trend for higher frequencies shown in forward calculations based on the inversion result was compared to the measured dispersion curve to find suitable models. Figure 9 shows the dispersion curve derived by the MASW survey and the comparison to the dispersion curve resulting from forward calculations using the velocity model from the inversion process (Fig. 9b). In Fig. 9c, the dispersion curve using a refined velocity model is shown in which only the top layer vS was changed from 186 to 400 m/s.

5.3 Site A1

Site A1 is located south-west of Temple C and north-west of Temple A. The center of the site is a crossroad surrounded by dense vegetation to the west and remains of ancient buildings and dry grass to the east (Fig. 3b). The DGPS elevation measurements revealed that the surface is dipping from north to south by 3.0°. The paths are used by pedestrians and show strong surface irregularities due to large rounded paving stones particularly on the north-south trending path.

5.3.1 Seismic refraction survey

Due to limits imposed by the existing ruins and vegetation the refraction profile at site A1 was measured on the west-east trending path on the Acropolis (Fig. 3a). The picked first breaks and corresponding travel time curves are shown in Fig. 10. Obviously, lateral inhomogeneity hampers the selection of straight travel time curves. A vP0 = 632 m/s of the top layer and vP1 = 1454 m/s for the refracted wave could be determined. The depth to the refractor is at 10.1 and 11.3 m at the western and eastern profile end, respectively, indicating a small inclination angle of 1.2°.

5.3.2 MASW

The dispersion curves based on the shots at both ends of the MASW profile (Fig. 11) are consistent at frequencies below 20 Hz but diverge at higher frequencies. The reason for this might be the same shallow heterogeneities as observed in the refraction data, probably due to changes in the substructure of the ancient path. A topographic depression of 0.6 m at the middle of the profile where two paths are crossing (Fig. 3) might also have an influence.

5.3.3 WARAN

The combined dispersion curve from the three arrays at site A1 does not show an inverse character for frequencies larger than 10 Hz. It fits well the general trend of the first mode for a subsurface model which includes a low-velocity layer.

5.3.4 Inversion results

For the inversion process, a start velocity model was prepared based on the Selinuntian stratigraphy. As no clear mode transitions were found, a simple model with dominant fundamental mode was used. Figure 12 shows the results of the inversion process.

5.4 Site W 1

Measurements at site W1 were carried out near Triolo Temple in the western part of the archeological park (Fig. 4) which at a height of 4.5 m a.s.l. is significantly lower than site E1, E2 (40 m a.s.l.) and A1 (25 m a.s.l.). The site is covered by dune sand (Fig. 4b). A tarred path in the east delimits the site from a small slope towards the river bed.

5.4.1 Seismic refraction survey

The first breaks and the travel time curves of the refraction profile (Fig. 4a) are shown in Fig. 10. Particularly for the second half of the profile, the uncertainty of picks increases due to first breaks fading into the noise; also, lateral heterogeneities cannot be excluded. However, velocities vP0 = 640 m/s and vP1 = 1132 m/s were determined with a thickness of the top layer of 2 m and no significant inclination.

5.5 MASW

The dispersion curves from the shot and reverse shot MASW experiments agree well at frequencies up to 40 Hz (Fig. 11) which confirms the negligible inclination of the top layer found by the refraction survey. The inverse trend of the dispersion curve again indicates the presence of a low-velocity layer.

5.6 WARAN

The dispersion curve from the WARAN measurement is combined from the results of arrays with 3, 10, and 20 m aperture. The inverse trend of the dispersion curve matches that from the MASW experiment.

5.6.1 Inversion results

In step (3), the SPAC analysis was applied to define frequencies at which potential mode transitions occur since distinct discontinuities of the dispersion curve derived by f-k are missing. The first potential mode transition at 2.5 Hz is outside the frequency range of sufficient resolution. For this reason, two different assumptions were followed in the inversion process: (1) the fundamental mode is not dominant in the frequency range of the measured dispersion curve at all, and (2) a dominant fundamental mode for the first part of the dispersion curve is considered between 5.2 and 11 Hz. Inversion results for both cases show similar layer velocities; differences mainly exist in the layer thickness. Figure 12 shows the inversion result considering case (2) which was also used for the forward calculation of synthetic seismograms.

6 HVSR

For all records from the WARAN, HVSRs were calculated (Nakamura 1989, 2000). Ratios between Fourier spectra of the root mean square (rms) of the two horizontal components and the vertical component were determined from time windows of 90 s length. An anti-trigger algorithm was used to exclude disturbances of the ambient noise records by local sources, e.g.. moving people or vehicles. The number of usable time windows for the seven surveys was between 200 and 1500. Cosine tapers of 5% of the transform window were applied and smoothing with a Konno and Ohmachi (1998) b-value of 20. The frequency range was from 0.8 to 50 Hz. The HVSRs were grouped into the four array measurements, the two profiles at sites W1 and A1, and the additional measurements on the Acropolis at Temple A (Fig. 3).

Figure 13 shows a summary of the HVSRs which were averaged for each of the seven surveys; the number of stations used is indicated in the figure. The profile across the Triolo Temple (W1 Profile) shows HVSR values close to unity for the whole frequency range. All other HVSRs from the survey exhibit the same character: large values or peaks at both, the low-frequency end and at the high-frequency end, and a clear de-amplification with values smaller than unity between frequencies of 0.37 to 0.63 Hz and frequencies of 25 to 40 Hz. This seems somewhat unusual for HVSRs; however, it is evidence for a low-velocity layer in the subsurface (Castellaro and Mulargia 2009). Di Giacomo et al. (2005) analyzed and modeled HVSRs of a case in Central Italy with a clear velocity inversion (4 m of anthropogenic fill on top of a 15 m thick hard conglomerate overlaying soft clays). Their measured and calculated HVSRs showed clear de-amplification between 1 and 8 Hz. In case of the Selinunte measurements, the low-frequency peak is close to the lower end of the spectra and the high-frequency peaks are close to the Nyquist frequency. Anyway, with a simple f = vS/4h estimate, the high-frequency peak correlates with the thickness h of the top layer; for example, in case of the HVSR for array E1, it is within 0.4 m of the inversion result (Table 3). The low-frequency peak might correlate with a 1.9 km deep velocity shear wave increase at the transition from clay to limestone as revealed in reflection seismic and borehole data (Ferranti et al. 2019).

Fig. 13
figure 13

Averaged HVSRs for the four array measurements, two profiles and stations located at Temple A. The legend shows the color code for the individual measurements, numerals in brackets give the number of stations used to calculate the average HVSRs

7 Discussion

The results of the nonlinear process of dispersion curve inversions cannot be expected to be unique even in cases of a good fit of the measured and theoretical dispersion curves (e.g., Forbriger 2003b). Results can be verified by an independent forward calculation of synthetic seismograms to simulate MASW sections based on the inverted velocity models. Such tests were performed using synthetic seismograms computed with the Qseis code by Wang (1999). Since dispersion curves are less sensible to changes in vP and densities (e.g., Forbriger 2003b), the vS models were complemented with the seismic refraction results and considering the geological constraints. The velocity models were adjusted until the dispersion curves from the inversion and the one obtained from the synthetic seismograms showed good agreement. Even in case of a good fit, an uncertainty of some 10% due to the inversion process in the velocities and layer thicknesses has to be considered.

7.1 Site E1

The thickness of the top layer in the vS profile based on the inversion process (1 m) is significantly smaller than indicated by the refraction survey (3.6 to 4.1 m). A self-evident explanation is the averaging effect over a large volume in the process of the dispersion curve measurement. The S-wave velocity for the halfspace (> 100 m depth) is not well resolved (Fig. 12). According to borehole data provided by Piro and Versino (1995), clay is expected at those depths, and thus, the determined velocity seems to be rather high. But a second, not overtaking refracted wave is observed (vP = 2330 m/s) which can be correlated with the measured halfspace velocity even though vS being overestimated in this case. An additional thin calcarenite or limestone layer is conceivable.

Inverse dispersion behavior is not observed although a low-velocity layer supposedly exists (Ruggieri and Sprovieri 1977). The forward calculation of a MASW section based on the velocity profile from step (5) analog to Fig. 9 shows an inverse trend of the dispersion curve and thus is in disagreement with the measuring results. Several reasons for this discrepancy are possible; according to Mi et al. (2018), energy produced due to guided waves in a low-velocity layer can disturb the dispersion curve. In addition, the influence of potential lateral heterogeneities on the dispersion curve has not yet been investigated in detail. Hence, a contamination of the dispersion curve cannot be excluded and the reliability of the obtained velocity profile is questionable. Since the temples are known to be situated on a calcarenite layer (Amadori et al. 1992), results from site E2 seem to be more appropriate to derive a suitable subsurface profile for the eastern part of Selinunte than those from E1, keeping also in mind the potential disturbance by ancient foundations and the agricultural use of the site. Thus, the velocity profile of E1 is not recommended to be used for further site effect studies referring to the temple structures in eastern Selinunte.

7.2 Site E2

While the S-wave velocity for the deeper layers at site E2 are well constraint by the inversion results (Fig. 12), the velocity for the shallow top layer shows a wide range of possible values. The forward calculation of synthetic seismograms and the determination of their dispersion curve helped to narrow down this velocity to 400 m/s. The parameters of the final velocity model are listed in Table 3 and shown in Fig. 14. The thickness of the calcarenite layer comes to 6.7 m which is in good agreement to the previous results of the reflection survey by Amadori et al. (1992). The P-wave velocities are taken from the refraction survey for the upper two layers. Based on the stratigraphy, the third layer is assumed to be sand which explains the low vS of 199 m/s. Suitable Poisson’s ratios μ range between 0.25 < μ < 0.4 for medium dense sand and 0.3 < μ < 0.45 for dense sand (Das 2013). Because any sand below the calcarenite layer is compacted, a Poisson’s ratio of μ = 0.4 was used to estimate vp to 490 m/s. The last layer resolved in the inversion process is referred to clay. The Poisson ratio for unsaturated clay ranges between 0.1 < μ < 0.3 (Bowles 1996). Since the layer is covered by several layers including the calcarenite, an assumed μ of 0.3 results in vP = 730 m/s.

Fig. 14
figure 14

Preferred shear wave velocity profiles for the western part of the Selinunte Archaeological Park (W1), the Acropolis (A1), and the eastern part (E2)

7.3 Site A1

Although the inverted dispersion curve for site A1 fits the measured curve well (Fig. 12), the forward calculation shows an inverse behavior which does not support the trend of the dispersion curve. However, the velocity model is compatible with the one for site E2, a site with similar geology, and a calcarenite thickness of 3.5 m is supported by the results of Amadori et al. (1992). Thus, it is reasonable to assume that this velocity model is representative for the Acropolis. The discrepancy with the forward calculation is possibly caused by the influence of lateral changes in the subsurface. The final velocity model is given in Table 3 and Fig. 14. P-wave velocities of the deeper layers were estimated with the same Poisson’s ratios used for site E2.

7.4 Site W1

The lower velocities derived for the subsurface at site W1 are caused by dune sand in the top layer. vS for the second layer is less than 500 m/s and therefore slower than the calcarenite layers found at sites E1 and A1. This is a strong argument that there is no near-surface calcarenite layer at this site. A profile excavated 2 m below the foundation of the Triolo Temple by Tusa et al. (1986) during the excavation in 1983 showed only sand layers. The differences to the other sites can be explained by the lower altitude of site W1 in comparison to sites A1, E1, and E2 (Fig. 1). Since only frequencies up to 40 Hz could be considered, the forward calculation was used to adjust the trend of the derived dispersion curve and small changes in velocities and layer thicknesses were applied. The final velocity profile is given in Table 3. The P-wave velocities of the top two layers are well constraint by the refraction survey. For the third layer again a Poisson’s ratio of 0.4 for sand is assumed which gives a vP of 350 m/s. The halfspace shows lower velocities than the second layer of site W1 and also lower than the clay layers of site E2 and A1. For sandy clays, 0.2 < μ < 0.3 is proposed by Bowles (1996) and with μ = 0.3 due to the depth gives a vP of the halfspace of 435 m/s.

7.5 HVSR

The HVSRs of the measurements support the existence of a velocity inversion with a clear de-amplification between some 0.37 to 0.63 Hz and 25 to 40 Hz. Di Giacomo et al. (2005) could show that such a de-amplification from a velocity inversion existed only in HVSRs from microtremor measurements and forward modeling but not in case of HVSRs from earthquake records, probably due to the different travel path of the waves from surface sources and sources at depth, respectively. Thus, the de-amplification should not significantly affect site-specific earthquake ground motions.

8 Summary and conclusions

Site effects are an important issue in quantitative archaeoseismic studies (e.g., Hinzen et al. 2016). Nowadays, SW Sicily shows low seismicity (e.g., Boschi et al. 2000; Guidoboni et al. 2007; Rovida et al. 2016), but the great temple structures of Selinunte were damaged by strong earthquake ground motions following Guidoboni et al. (2002) and Bottari et al. (2009). The determination of seismic velocity profiles through the combination of diverse in situ measurement and data processing methods for the three differently situated parts of Selinunte are the first step in site characterization and essential for planned temple response investigations.

At all three locations, the P-wave velocity of the upper two layers were obtained by seismic refraction experiments and the S-wave velocities were derived by inversion of dispersion curves from MASW and ambient noise array measurements. The results were checked by forward modeled synthetic seismograms and adjusted if necessary.

With regard to the subsurface of the temple structures, the velocity profile of site E2 is recommended for the eastern part where the Temples E, F, and G are located. Results from site A1 represent the subsurface conditions for Temples A, B, C, and D and probably also for Temple O. The models from sites E2 and A1 follow the general stratigraphy of Selinunte (Ruggieri and Sprovieri 1977). The calcarenite thickness of 6.7 m for site E2 and 3.5 m for site A1 are in agreement with previous studies by Amadori et al. (1992). A de-amplification of HVSRs with values clearly below unity support the existence of a velocity inversion at the base of the calcarenite layer. As the temples are presumably positioned mainly on the firm calcarenite (Amadori et al. 1992), the overlaying deposits should be neglected when synthetic seismograms for dynamic studies of the temples are calculated.

For the western part of Selinunte the final velocity profile of site W1 shows no shallow calcarenite layer and thus different site effects can be expected. The final S-wave velocity profiles are shown in Fig. 14. While the refraction measurements cannot resolve any low-velocity layers, for the S-wave velocity, this became possible by the inversion of the dispersion curves; however, the latter gives less constraints for the P-wave velocities.

The combination of MASW and WARAN measurements allowed inverting a larger frequency range of the dispersion curve. However, heterogeneities in the subsurface lead to complex dispersion curves of the MASW measurements. Further, low-velocity layers cause inverted trends of the dispersion curves with several modes visible which are increasingly difficult to distinguish the higher the frequencies. Nonetheless, the overall trend of the measured inverse dispersion curves was used in combination with forward calculations of a MASW survey to optimize the velocity profiles (Figs. 15, 16, 17, and 18).

The dispersion curve derived by MASW measurements gives an average of the chosen profile. Bignardi et al. (2014) emphasizes that in case of heterogeneities in the subsurface, the MASW result depend on the shot direction. In contrast to MASW, the dispersion curves from the WARAN data give an average of the whole area covered by the array aperture. Thus, small local lateral heterogeneities do have less impact on the end result. In the survey area, it was of advantage to combine MASW and WARAN surveys. Since the inverse trend of the dispersion curves leads to reduced invertible frequency ranges, the results of MASW in combination with forward modeling helped to constrain the top layer velocities which are bound to higher frequencies. The use of a wireless network was a great plus as it allowed to quality control of the data during the ongoing measurements and position stations in rough terrain.

Further MASW profiles parallel to the measured spread would generally improve the resolution of the subsurface and its heterogeneities. However, due to rough terrain as well as archeological remains, this is not practical for all sites, particularly for the Acropolis. Strategically, well-placed boreholes could help to better constrain the inversion calculations.