© Author(s) 2010. CC Attribution 3.0 License. Natural Hazards and Earth System Sciences

Abstract. We examine the possible non-linear behaviour of potentially liquefiable layers at selected sites located within the expansion area of the town of Nafplion, East Peloponnese, Greece. Input motion is computed for three scenario earthquakes, selected on the basis of historical seismicity data, using a stochastic strong ground motion simulation technique, which takes into account the finite dimensions of the earthquake sources. Site-specific ground acceleration synthetics and soil profiles are then used to evaluate the liquefaction potential at the sites of interest. The activation scenario of the Iria fault, which is the closest one to Nafplion (M=6.4), is found to be the most hazardous in terms of liquefaction initiation. In this scenario almost all the examined sites exhibit liquefaction features at depths of 6–12 m. For scenario earthquakes at two more distant seismic sources (Epidaurus fault – M6.3; Xylokastro fault – M6.7) strong ground motion amplification phenomena by the shallow soft soil layer are expected to be observed.


Introduction
The prediction of the soil response in the presence of potentially liquefiable layers during an earthquake can importantly contribute to the seismic hazard assessment and the design of building requirements.The fact is that liquefiable layers can considerably influence the characteristics of the ground motion.In particular, due to the soil softening within the liquefiable layers, ground response is altered.The most im-Correspondence to: T. Novikova (tatyana0703@yahoo.com)portant consequence is the change to the spectral content of the ground motion for which the structure has been designed.
In about the last 40 years or so, considerable damage to buildings, roads, bridges, earth dams and other structures was caused by liquefaction events during earthquakes, such as the ones of Wildlife 1987 Superstition Hill, California;Treasure Island, 1989 Loma Prieta, California;Kushiro Port, 1993 Kashiro-Oki, Japan;Port Island, 1995 Kobe, Japan;Yuan Lin, 1999 Chi-Chi, Taiwan; Sakai-Minato city, 2000 Western Tottori, Japan (Seed et al., 1990(Seed et al., , 2000;;Bardet et al., 1995;Comartin et el., 1995;Sitar, 1995;Somerville, 1995;Mori and Sogabe, 2002).These events and also the fact that many major cities around the globe are extensively built on young saturated sediments, alerted scientists about the potential risk associated with the soil liquefaction.
In Greece, a highly seismogenic region, many liquefaction cases due to strong earthquakes have been documented (Papadopoulos and Lefkopoulos, 1993;Papathanassiou et al., 2004).The Nafplion city, which has been selected as the test area in this study (Fig. 1, right), initially was founded on safe ground.However, its lateral expansion in recent years has been made along the coastal zone on ground of questionable safety factor as regards the earthquake and liquefaction risks.In the Argos plain, where Nafplion is situated, at least one liquefaction case has been reported in the past.This was caused by a strong intermediate-depth earthquake (M = 7) which occurred on 2 June 1898, with its epicentre placed at about 30 km NWW from Nafplion (Papadopoulos and Lefkopoulos, 1993).
Moreover, the proximity of the city to active seismogenic faults should be taken seriously into account in the consideration of the seismic hazard assessment.In the present work we focus on the study of the nonlinear ground response effects with particular attention on how potentially liquefiable shallow layers may influence the ground motion characteristics during strong earthquakes at specific sites within the area of expansion of Nafplion.The Nafplion case is suitable to demonstrate how we can effectively predict the events with high probability for liquefaction initiation, following a practice based on the nonlinear analysis, using geological information, geophysical survey results combined together with a number of parameters for hypothetical strong seismic events obtained from stochastic simulation.
The procedure followed has two main steps.The first step is the evaluation of the liquefaction hazard in probabilistic terms by taking into account available geophysical and geotechnical data, i.e. seismic shear-wave velocity (V s ) profiles, Standard Penetration Test (SPT) values to characterize the soil at the sites of interest and utilizing the results of stochastic simulations, the Peak Ground Acceleration (PGA) in combination with earthquake magnitude values.The second step is the evaluation of the influence of the local site conditions on the ground shaking characteristics through a site response analysis in terms of response spectra, acceleration, strain and Excess Pore Water Pressure (EPWP) time histories at different depths.This is accomplished with the aid of nonlinear modelling, which commonly tracks a seismic load through stress-strain space using different assumptions for the details of the stress-strain relationship.
It is also common, to approach such problems with equivalent linear approximation, which requires only few, generally well-known parameters, as shear wave velocity, modulus reduction and damping V s shear strain and soil density.An equivalent-linear approximation is quite simple and assumes that a damped linear model can approximate dynamic soil behaviour through the appropriate choice of material parame-ters.Nevertheless, the respective models do not account for the accumulation of plastic deformation, the strain rate dependence, the influence of stress history on soil stiffness and the evaluation of the pore water pressure, which are points of great importance as regards potentially liquefiable layers response.Yoshida and Iai (1998) pointed out another two weaknesses of the equivalent-linear method.The first is that the damping is constant and independent of frequency.The second is that the resonant frequency amplitudes and shear stresses are overestimated.
In practical applications, equivalent linear analysis has a tendency to give larger peak acceleration and shear stress under large earthquakes and low amplification in high frequency range.Nonlinear analysis provides a more robust characterisation of the true nonlinear soil behaviour.Nevertheless, there are only few applications in practice, since the parameter selection and the code protocols in use are often too complicated.We finally adopted the implementation of the nonlinear analysis since all the required data from the field tests were available including Standard Penetration Testing, laboratory tests on borehole core samples as well as results from geophysical measurements along with the powerful computational tool D-MOD2000.
2 Description of the analysis procedure

Geological site characterization -borehole testing
Realistic assessment of the liquefaction potential in an area requires knowledge of the geotechnical characteristics of its shallow geological formations.In the following we briefly describe available information on the soil formations of the study area stemming from previous studies.
The area under investigation (Fig. 1) is formed by alluvial, mainly lagoonal deposits, overlain to flysch and limestone bedrock.As regards the groundwater regime, within the Quaternary deposits successive groundwater aquifers are developed, being under intensive exploitation by well boring.This caused considerable sea water intrusion in recent years.Within some parts of the investigation area, at an altitude laying a few metres above the sea level, a weak unconfined coastal aquifer is developed at a small depth near surface, which is underlain by deeper confined aquifers.This shallow unconfined aquifer belongs to a local marshland.
Detailed borehole testing was carried out at two sites of the investigation area (G1 and G2 in Fig. 1), which included lithology examination, standard penetration testing, permeability tests and sampling for laboratory tests.The samples were granulometrically examined in the laboratory and the percentages of the liquefaction prone materials, mainly sand and silt, were accurately defined.In addition, a set of geotechnical parameters, such as Atterberg limits, natural moisture, special weight etc. was also acquired.
The lithological recognition of the sedimentary formations showed that the investigation area was mainly covered by sand-silty sediments down to 25 m depth (Apostolidis and Koutsouveli, 2007).In particular, for the G2 borehole, the concentration of the sand became very important at depths of 4.9-8.55 m.The laboratory analysis of samples from that borehole showed that sand concentration, at the depths of 4.90-5.55m and 7.90-8.55m, is about 29% and 66%, respectively.In both boreholes, however, the clay concentration did not exceed the one-digit values for all the samples from depths shallower than 25 m.It is notable that in the shallow depth samples the concentration values of clay was measured to vary from 1% to 5%.

Establishment of V s soil profile by the use of geophysical survey
In the framework of a seismic hazard assessment study in the investigation area, an extensive geophysical survey with gravity and seismic methods was conducted aiming mainly to the detection of possible hidden faults but also to the characterization of the soil formations.The geophysical survey, the methodologies used and the results are presented in detail by Karastathis et al. (2010), however here we are only giving the necessary description of the methods and results utilized for the liquefaction risk assessment.Figure 1 depicts the sites of the geophysical investigations.The seismic shear wave velocity (V s ) needed in our methodology for the calculation of the strong motion amplification due to soil conditions, is usually estimated by seismic methodologies carried out at surface or alternatively through boreholes.More specifically in the present case, the geophysical survey was based on the active seismic techniques of Multichannel Analysis of Surface Waves (MASW) (Park et al., 1999;Xia et al., 1999) and seismic refraction but also on the passive technique of Microtremor Array Measurements (MAM) (Hayashi and Kita, 2009) as well as on crosshole measurements performed on couples of boreholes.

Active seismic techniques
The Multichannel Analysis of Surface Waves (MASW) method is one of the most frequently used in similar cases in the last decade.The method is based on the dispersion properties of the Rayleigh type surface waves, i.e. the change of the phase velocity with frequency.The fundamental mode dispersion curve (phase velocity vs frequency) is derived by various techniques (McMechan and Yedlin, 1981;Gabriels et al., 1987;Park et al., 1999).As a final stage of the method, an inversion technique is implemented using the measured dispersion curve as reference and by iterative inverse modelling constructs the S-wave velocity model that can justify this dispersion.MASW can be used either for 1-D or 2-D surveys.In our case, what was required was only the onedimensional distribution of V s with depth at several places within the investigation area; thus, the application of 1D-MASW was sufficient.However, if there were abrupt variations of Vs in short distances the application of 2D-MASW would have been considered.Figure 2a presents an example seismic record acquired from a MASW survey conducted at WW site.The seismic source of Accelerated Weight Drop (AWD) provided energy of adequately low frequency for a clear determination of the fundamental mode phase velocity curve on the spectrum V phase (f ) as low as 4 Hz frequency (Fig. 2b).In addition, the selection of low frequency geophones (4.5 Hz) was also contributed to this clear spectrum.After iterative inversion processing we reached to a 1-D S-wave velocity model (Fig. 2c) that could produce a similar synthetic V phase (f ) curve.
The resulted 1-D S-wave velocity models from all MASW investigations conducted at the Nafplion expansion area are gathered and presented in Fig. 3a.The profiles QQ and PP conducted in the old city of Nafplion are presented in Fig. 3b.
At DD , BB , QQ and PP sites the selected lay-out length was 72 m.At these sites beyond MASW technique, the seismic refraction of P and S-waves was also implemented.The additional information from the P-wave survey was critical since it assisted the estimation of the depth of the saturated zone.When the P-wave velocity exceeds the value of 1.5 km s −1 and at the same time the Poisson ratio reaches values as high as 0.49 the presence of a fully saturated layer can be assumed.
Figure 4a and b presents two example records from the seismic survey DD .The record at the top (Fig. 4a) is from P-wave survey and at the bottom (Fig. 4b) from the respective S-wave one.The first arrivals (direct and head waves) can be clearly picked.For the P-wave production the Accelerated Weight Drop was utilized.On the other hand a 7 kg weighed hammer striking a wooden beam used for the generation of horizontally polarized shear waves.The two datasets of the P and S-wave arrivals were processed by seismic refraction tomography and inversion modelling techniques (Hayashi, 2003) resulting to the models shown in Fig. 4b and c, respectively.The stratigraphic model of Fig. 4c and the traveltime curves (Fig. 4d) show that the most of the head wave arrivals come from the water-table interface since the seawater intrusion causes an abrupt increase of the P-wave velocity.Figure 4c shows also the very good fit between the observed and synthetic (calculated) traveltimes for the stratigrafic model of Fig. 4b. Figure 4e and f shows the respective model and traveltime curves for the S-wave profile.The fit between calculated and observed times is also good for the shear waves.
The methodology of seismic refraction although requires monotonic increase of the seismic velocity with depth, is generally a very reliable geophysical tool.The comparison of its results with the ones of MASW method enhanced the reliability of MASW, which is highly dependable on the selection of the fundamental phase velocity curve.The results of MASW at DD site is in full accordance with the model   resulted by seismic refraction technique.Respective models for BB , QQ and PP were resulted after similar way of processing.MASW calculated the 1-D S-wave velocity models for these sites and are presented in Fig. 3.However, since QQ survey was located quite far from the coast its results are worth commenting.At this site no shallow aquifer was detected so this area is in a different condition than the costal one.More specifically the analysis of the P and S-wave arrivals resulted to the model of Fig. 5a with very good fit between the observed and calculated traveltimes and also in a full agreement with the MASW results (Fig. 3b).The P-wave velocity values and also the Poisson ratio do not justify the consideration of an aquifer at a shallower depth than 10 m.However, at greater depths the values can justify the presence of water in the porous media but the subsurface material seem to be quite cohesive.
In the old city of Nafpion a seismic survey was carried out at the site PP in order to compare the results with the ones of the city's expansion area.The refraction investigations detected the limestone bedrock at very shallow depth, down to 15 m (Fig. 5b).The aquifer is located at a depth about 3-4 m below loose sediments.
Although the main aim of the investigations at the sites AA , CC , EE , FF , WW was the fault detection and the seismic reflection method was mainly selected and implemented (Karastathis et. al., 2010), we also acquired data usable for MASW technique.The surface wave energy produced by the Accelerated Weight Drop was adequately strong to be recorded at distance longer than 200 m (Fig. 6a, b).The dispersion curves were clearly defined at the respective spectra (Fig. 6c).The results are shown in Fig. 3.

Crosshole testing
Crosshole testing (ASTM D4428/D4428M-07 standard) was applied to determine the P-wave and S-wave velocity at two sites (G1 and G2, see the map of Figs. 1 and 7).The distance between the boreholes of each pair was 5 m and the depth reached was 40 m.The seismic source used was a sparker.The receiver was a triaxial borehole geophone.The crosshole testing provided direct measurements of the seismic wave velocity in the mid space between the two boreholes.The measurements were combined with the results of SPT and those of the laboratory testing.The high validity of the borehole tests assisted to the evaluation of the validity of the surface seismic surveys.

Microtremor Array Measurements
The passive method of Microtremor Array Measurements (MAM) (Hayashi and Kita, 2009) based on SPAC technique (Okada, 2006), uses microtremor recording by a multichannel array to analyze phase-velocities.In the present study the method assisted to the identification of the fundamental dispersion curve.The MAM survey was based on a rectangular angle shaped layout (L) with a number of geophones up to 11 in passive recording.Figure 8a shows an example of a dispersion curve from a MAM conducted at the site of borehole G1.The spectrum point to a fundamental curve similar with those resulted by MASW (Fig. 6c).The MAM method clearly described the fundamental dispersion curve for frequency as low as 2 Hz yielding a greater investigation depth in respect to active MASW.Nevertheless, the main reason for conducting the passive survey was to investigate the 1-D shear wave velocity at an old 100 m deep borehole (site 3042 on the map of Fig. 2) with known geological log, close to a heavy traffic road.In order to evaluate the precision of the method we compare the results of the passive method with the crosshole tests at the sites G1 and G2.The comparison was very encouraging since MAM resulted to a 1-D shear wave velocity model with  values similar to those of crosshole tests.Figure 8b shows the comparison between the two methods at the site G1.The low velocity zone at depth about 30 m was detected by both methods.Figure 8c shows the 1-D S-wave velocity model as resulted by MAM at the site of the old borehole 3042 in comparison with the results of MASW at the neighbor site FF .The results are in a good agreement in this case as well.

Definition of the soil profile by geotechnical study
The profiles for the investigated sites were established mostly on geophysical measurements but also taking into account SPT and laboratory tests on borehole core samples (Apostolidis and Koutsouveli, 2007).The profiles at borehole testing sites G1 (a) and G2 (b) are shown in Figure 7.

Selection of scenario earthquakes
Scenario earthquakes examined herein were selected based on available information on the past seismicity within and around the area of interest.We searched both historical and instrumental seismicity archives for events that have affected the built environment of the town of Nafplion in the past.We limited our search on shallow earthquakes as most widely applied strong ground motion simulation techniques have been tested almost exclusively on such events.Although the seismicity in Greece is well known to be high, the eastern part of Peloponnese, where Nafplion is situated, has not suffered by many catastrophic earthquakes.In fact, during instrumental times no significant earthquakes have been recorded from this area.However, in historical seismicity catalogs one can find at least three moderate-tolarge magnitude earthquakes (on 20 March 1837, 27 June 1769 and on 31 January 1742) that had impact on the town of Nafplion at the time of their occurrence.According to Papazachos and Papazachou (2003) these events had magnitudes of 6.3, 6.4 and 6.7 and are related to the following fault structures, respectively: 1.The Epidaurus fault.Papazachos and Papazachou, 2003) and.

The Iria fault (or Argos fault in
3. The Xylokastro fault in the southeastern cost of the Gulf of Corinth.
In this paper we choose to study scenario earthquakes that are analogue to these three, significant for the town of Nafplion, historical events.Therefore, the scenario magnitudes and sources of the scenario earthquakes are the same as those that have been assigned to the three historical events.Details on the relative simulation parameters and appropriate references will be given in the following.

Simulation method and application
To simulate the strong ground acceleration at the sites of interest we applied the stochastic method for finite sources as has been proposed by Beresnev andAtkinson (1997, 1998).The specific method was chosen due to its simplicity in application since available information on the earthquake sources and the regional and local seismic wave attenuation characteristics do not permit a more deterministic approach.Furthermore, the applicability of the method has been shown in past studies to be quite successful throughout a wide range of frequencies and in various seismotectonic environments (e.g.Beresnev and Atkinson, 1998a, b, 1999, 2002;Hartzell et al., 1999;Castro et al., 2001;Hough et al., 2003;Roumelioti and Kiratzi, 2002;Roumelioti et al., 2004).
The applied method involves discretisation of the fault plane into a certain number of subfaults, each of which is assigned an ω −2 spectrum.Each subfault is triggered when the rupture front reaches it.Contributions from all subfaults are empirically attenuated to the observation site and appropriately summed up to produce the synthetic accelerogram.For a detailed description of the applied method and the parameters involved in its application (FINSIM code) readers are referred to the original work of Beresnev andAtkinson (1997, 1998a).
Basic simulation parameters in the application of FINSIM refer to the earthquake source, the path of the seismic waves and the site effect.
As mentioned above, the earthquake source is modelled as a plane and one needs to assign its location, dimensions, strike, dip and depth of its upper edge from the surface.For all three scenarios, locations and geometry of the sources were based on information given by Papazachos and Papazachou (2003).The depth of the upper edge of the fault planes was empirically set to 1 km, i.e. deep enough to avoid placing large values of slip at the ground surface but also shallow enough to be close to empirical data in Greece which suggest that earthquakes of magnitudes similar to those examined in this work may or may not rupture the ground surface (Pavlides and Caputo, 2004).Fault dimensions were computed based on the empirical relations proposed by Wells and Coppersmith (1994), which relate the length (L) and width (w) of a fault plane to the moment magnitude (M w ) of the earthquake.Finally, the number of subfaults into which each fault plane was divided was also determined using an appropriate empirical relation that has been suggested by Beresnev and Atkinson (1999).
Parameterisation of the propagation path includes empirical description of the geometric attenuation and the anelastic attenuation of the seismic waves.For the geometric attenuation we applied a geometric spreading operator of 1/R, where R is the distance from the seismic source, while the anelastic attenuation was described through a frequency-dependent quality factor, Q(f ) = 100f 0.8 , applicable to the broader Aegean area (Hatzidimitriou, 1993(Hatzidimitriou, , 1995)).Fig. 9. Seismic source models adopted in the stochastic strong ground motion simulation method for the three examined earthquake scenarios.Model fault surfaces are devided into subfaults.Yellow asterisks denote the locations of the assumed rupture initiation points (assumed to coincide with the hypocenter).
Near-surface attenuation of the seismic waves was modelled using the kappa, κ, operator of Anderson and Hough (1984) and diminishing the simulated spectra by exp(−π κf ).Simulations at all sites of interest were performed assuming rock site conditions at the ground surface.Site-specific ground characteristics were taken into account in the subsequent step of the liquefaction risk analysis.The value for the kappa operator for rock adopted herein is the one proposed by Margaris and Boore (1998).
Parameters adopted in each of the three sets of simulations are summarized in Table 1.Surface projections of the three fault models are shown in Fig. 9.In this figure each plane appears divided into subfaults and the asterisks on the surface projections denote the assumed hypocenters.Hypocenter locations are, of course, not known for historical earthquakes and therefore their choice for the applications presented in this work is not based on any kind of data.We did try however to simulate "worst-case" conditions in terms of rupture directivity i.e. in cases where the rupture of the fault surface could direct toward the town of Nafplion (scenario earthquakes on Iria and Epidaurus faults) we selected the hypocenter location that would maximize the phenomenon.Regarding the slip on the adopted fault models we used random distributions as we cannot a priori know the details of rupture of a future event.The details of slip play an important role in the distribution of strong ground motion within the surface projection of a fault but this effect fades out with distance from the fault.In our simulations all sites of interest  lie way out the projections of the adopted fault models and thus the distribution of slip is not expected to affect simulation results in a significant way.The product of the stochastic simulations, which was used as input in subsequent steps of our analysis, was a set of synthetic acceleration time histories (S-wave part only) at selective sites within the study area.
In previous studies (e.g.Stewart and Tileylioglu, 2007) it was shown that input motion at the bottom of the 1-D site profile in nonlinear analysis should be specified as "outcropping" (i.e., equivalent free-surface motions) with an elastic base or "within" (i.e., motion recorded at depth in a vertical array) a rigid base.In our study we have employed the "outcropping" one.
Synthetic PGA values from the three scenario earthquakes and at each one of the seven examined sites are included in Table 2.These values are a prerequisite for the subsequent steps of the liquefaction potential assessment.

Evaluation of liquefaction risk
The evaluation of the liquefaction potential first is approached by calculating the factor of safety (FS) against liquefaction.FS is usually expressed as a ratio of cyclic resistance (CRR, soil "strength") based on in-situ test data (i.e.SPT, BPT, CPT) to the average cyclic stress (CSR, earthquake "load") induced in the soil by an earthquake: FS = CRR/CSR.The earthquake demand (CSR) is calculated by the so-called Seed's method (Seed and Idriss, 1971).The equation is as follows: where the constant value 0.65 is a weighting factor, introduced by Seed and Idriss (1971), to calculate the number of uniform stress cycles required to produce the same pore water pressure increase as an irregular earthquake ground motion; σ 0 is the total vertical overburden stress; σ 0 is the effective vertical overburden stress based on water table during earthquake; α max is the Peak Horizontal Ground Acceleration, PGA (in g units); r d is the stress reduction coefficient, varied with depth.The value of r d at the depth of z can be calculated using the following equation (NCEER, 1997): (2) To determine the dynamic geotechnical properties and evaluate the liquefaction resistance (CRR) of soil, we followed the shear wave velocity (V s )-based approach as V s profiles were available from previous geophysical studies in the area (seismic methodologies, mostly MASW and crosshole testing, Karastathis et al., 2008).The V s -based determination method is particularly useful in soils where the sampling is difficult, e.g.gravelly soils, or the penetration tests may be unreliable, and also in case where boring may not be permitted, like in capped landfills.Moreover, another advantage of V s approach is that it directly accommodates profile parametric uncertainty in a statistically rigorous manner: the equation for determining the CRR from V s is empirical and based on case history studies at sites that did and did not liquefy during earthquakes (Andrus and Stoke, 2000): where, MSF = (M w /7.5) −2.56 is the magnitude scaling factor; V S1 is the stress-corrected V S and defined as V s 1 = V s (Pa/σ v ) 0.25 ; where V s is the measured shear-wave velocity (m/s), Pa is the reference stress (100 kPa), σ v is initial effective overburden stress (kPa); V S1C is a correction factor that depends on fines content; K C is a correction factor for cementation and aging; M w is the earthquake magnitude.
In present study we took the value of K c equal to 1, since there is currently no widely accepted method for estimating K c and its variability across the category areas (Andrus and Stokoe, 2000).The actual fines content varies with depth and location.However, in the analysis introduced by Andrus and Stokoe (2000), Eq. ( 3) was derived from sands and gravels that were classified into three broad categories with regards to fines content: ≤5%, 6%-34%, and >35%.In our case studies, according to the data of standard penetration testing, the content of fines of the liquefiable sandy silt ranges from 16% to 60%.
FS is the ultimate result to receive for liquefaction analysis: for FS≥1 there is no potential of liquefaction but for FS<1 liquefaction is expected to occur.However, it is common practice to also use the probability of liquefaction, P L .We, thus, quantified the potential for liquefaction by expressing the factors of safety in terms of probability.One of the important advantages of that method is that the probability of liquefaction is information required for the risk-based design decisions.Probability of liquefaction is obtained using the equation developed by Juang (2001): This analysis helps us to distinct the layers with high liquefaction potential within the structure and to evaluate further their nonlinear response.

Nonlinear response of potentially liquefiable layers
The high potential for liquefaction initiation at some selected sites of Nafplion city during scenario earthquakes could be associated with nonlinear response of the ground.Therefore, the nonlinear stress-strain behaviour of the soils should be taken into account in the site response analysis.There is a variety of codes capable of performing fully nonlinear, effective stress based ground response analysis.Among them such codes as DynaFlow (Prevost, 1983), SwanDyne (Madabhushi and Zeng, 1998) that are codes for the static and transient response of linear and nonlinear 2-and 3-D systems; CYCLIC 1D and Cyberquake (CyberQuake, 1998) which are quite simple 1-D computer programs for simulation of earthquake site response.We accomplish our goal with the aid of the fully nonlinear 1-D code D-MOD2000 (Matasovic and Ordonez, 2007), which is the most recent version of D-MOD 2 (Matasovic, 1993(Matasovic, , 2006)).The algorithm provides a built-in effective-stress analysis option (cyclic degradation of material properties with hydraulic interaction between layers, i.e. dynamic response plus pore water pressure generation plus pore water pressure dissipation and redistribution).
It must be noted that the algorithm requires the availability of nonlinear material parameters, such as index properties and density.The availability of those parameters allowed implementation of the algorithm in our study case.
At each time step, D-MOD2000 refers to the given stressstrain relation to obtain a value of shear stress, τ , from the current value of shear strain, γ .Because the shear stress is computed directly as a function of shear strain, the slope of the stress strain relationship, that is the secant shear modulus, G sec , does not have to be obtained explicitly.The new of the shear modulus is then used in solving the equation of motion at the next time step.With this scheme mixed boundary conditions can be easily handled as the stresses, displacements or velocities at the boundaries are simply set to the required values.For the one dimensional shear wave propagation problem, the boundary conditions are different for the top and bottom layers but are otherwise the same for all intermediate layers.The pore water pressure generation is controlled by semi-empirical models for sand (Dobry et al., 1985;Vucetic andDobry, 1986, 1988) and by a model for clay (Matasovic and Vucetic, 1995), with a computed pore water pressure used to degrade the backbone curve in a way that represents the softening and weakening of soils expected during strong ground shaking.
In comparison with the codes that account for the nonlinear behaviour of the soil using an iterative procedure, D-MOD2000 has the following three advantages: (i) there are no restrictions on the input acceleration level; in fact, D-MOD2000 incorporates a nonlinear constitutive model, while codes which are based on an equivalent-linear model that typically extends up to 1% shear strain level.In large earthquakes and soft soils, shear strains in excess of 1.0% can be induced; (ii) D-MOD2000 can directly calculate soil liquefaction potential given that D-MOD2000 is an effective stress program; (iii) D-MOD2000 accommodates for pore water pressure induced "soil softening" (reduction in shear modulus and strength).
The program D-MOD was extensively tested and validated by numerous researchers.The examples of its application include direct comparison with the closed form solutions for simple excitation cases (Kwok et al., 2007), blind prediction of the site response in the 2004 Parkfield earthquake (Kwok et al., 2006), and revisit to the La Cienega and Treasure Island case histories (Stewart et al., 2008).
The implementation of the D MOD2000 code, for the case of nonlinear ground response analysis in the area of Nafplion, required to define a-priori the three fundamental sets of soil properties that would be used further to create analytical soil profiles for the investigated sites.The first was the shear wave velocity profile, obtained here by seismic surveying, described in previous section.The second was basically the set of curves describing the nonlinear relationship, shear stress vs. shear strain, and hysteretic damping ratio vs. shear strain.The former is expressed in normalized form as modulus reduction curves (G/G max curves) and the later as damping versus strain curves.In the present analysis we selected this set of curves for each model structure employing the EPRI curves (1993), which are based on the depthdistributed velocity.The third one was a set of values for density, ρ.
According to the results of borehole testing, the investigation area is mainly covered by sand-silty sediments down to 25 m depth and below this depth clayey silt material was recognised.Based on this information, in the present study three identical material types were used within profiles for every site: sandy silt with small percentage of clay and fine gravel (1); sandy silt with increased percentage of gravel and pebbles (2); and clayey silt (3).The range of depths for every material type, as well as Vs distribution differs from site to site (see Figs. 2 and 3).The density values were obtained from the results of laboratory test on borehole samples.For material type 1 it is equal to 17.82 kN m −3 ; for material type 2 the value corresponds to 18.15 kN m −3 and for material type 3 is 19.89 kN m −3 .Density variation was not considered within the layers of every material type.
The initial tangent shear modulus G m0 has been determined from the average shear wave velocity profiles and density of the soil layer.The values of τ m0 are determined from the modulus reduction curves from every case study.Coefficients β and s are determined by fitting modulus reduction curves of each material.The value of 1 was chosen for the coefficient υ, that represents the original stress degradation law (Matasovic and Vucetic, 1993).
Coefficients K 2 = 0.0025, m = 0.43 and n = 0.62, used to calculate the rebound modulus, were selected on the basis of recommendations given in Martin et al. (1975).Rayleigh damping coefficients α R and β R are determined for constant viscous damping ξ , and the period of oscillation of T .The values of ξ are estimated for each profile with the aid of calibration procedure described below and the values of oscillation period are determined on the basis of shear beam theory and the average shear wave velocity and mass density of the soil in the profile.
To avoid the error which can be caused by under-damping of the system, we first calibrate D-MOD in the following manner: run a SHAKE2000 analysis to obtain response spectra at the surface; then use the same input information of SHAKE to perform a total stress analysis with D-MOD; subsequently compare SHAKE and D-MOD response spectra at the surface, and change the damping coefficients until we get a reasonable match; finally proceed to the nonlinear effective stress analysis employing D-MOD.
The values of Dobry's pore water pressure model parameters p, F , s and f we selected based on the recommendations given by Dobry et al. (1985).

Estimation of liquefaction probability
The combined application of geophysical, borehole and laboratory tests resulted to the identification of the layers exhibiting nonlinear behaviour and are, thus, liable to the initiation of liquefaction in the test area.More specifically, the mechanical parameters of the materials, such as V s and N 60 of SPT, were actually associated to the lithology, the composition of layers being considered according to the sandsilt-clay concentration.The aquifer was detected at shallow depths.The near surface sand-silty sediments were fully saturated as resulted from the joined application of the seismic surveys of P-and S-waves.More precisely, when the V p /V s ratio is so high to lead the Poisson ratio up to 0.49, we can consider the full saturation of these sediments.Figure 3 provides a good example of this by showing the variation of the dynamic Poisson ratio.It is worth noting also that Vs has very low values in the near surface saturated sand-silty layers, which perhaps indicates high liquefaction potential.We also give an example of using microtremor data (1-D velocity profile) in liquefaction analysis.Although we carried such experiment for three locations (3042, G1, G2), we found out that using MASW methodology is a lot more effective.Thus, here we perform the results only for location 3042.
The results of the liquefaction analysis performed in the frame of the present work are summarised in Table 3.The table includes information on the depth of the liquefiable layers, the Factor of Safety (FS) and the Probability for Liquefaction (P L ) at the investigated sites for the three examined earthquake scenarios.
At first, the depth distribution of FS and P L was calculated, for the scenario earthquake of Epidaurus fault (M = 6.3).As described in a previous section, the input motion was derived by stochastic modelling for bedrock outcropping site conditions.At site CC, between 6 and 10 m depth, P L reaches the value of 60%.The respective value for AA site is about 40%.P L is up to 50% and 25% at the sites of the boreholes G2 and G1, respectively.The probability of liquefaction at the remaining sites (DD, EE, BB, WW , 3042 and FF ) was found to be negligible (0-15%).
Very high P L values were estimated for the scenario earthquake on the Iria fault (M = 6.4).More specifically, at both borehole sites, G1 and G2, P L reached almost 100% at depths between 6 and 8 m.In addition, at the CC site, P L was esti-  mated to range between 40 and 100% at 2-10 m depth.Sites FF and AA also appears to be a hazardous sites, since P L reaches 70%-80% at depths between 4 and 12 m.Liquefaction is less possible at DD site (40% at 4-12 m), BB (10-30% at 6-12 m), EE (20-40% at 8-16 m), WW (20-48% at 8 m) and 3042 (10-20% at 7.5-10 m).
The scenario earthquake of M = 6.7 on the Xylokastro fault is less hazardous as regards the liquefaction initiation.At the two borehole sites G1 and G2, P L is of the order of 40-45% at 6-8 m depth.At site CC P L is 67% at 6-10 m depth and at sites FF , AA the respective value did not exceeded the value of 25-30% at 4-12 m depth.At the rest of the sites the probability for liquefaction for this scenario was negligible.

Evaluation of 1-D nonlinear ground response
Eventually, the influence of the local soil conditions on the ground shaking characteristics is evaluated by applying site response analysis in terms of response spectra, EPWP and acceleration time histories as well as stress-strain histories at different ground level.Through the aforementioned we aim to gain a better understanding of the nonlinear soil behaviour and its influence on ground motions.The analysis was based on the estimated liquefaction susceptibility, on the material properties measured and on the input motion stochastically modelled.

Scenario earthquake of the Epidaurus fault (M=6.3)
P L from the linear equivalent analysis at sites AA, CC and G2 is estimated up to 60%.The nonlinear analysis, however, showed that, although excess pore water pressure is generated during this earthquake at some layers, the PWP ratio does not exceed 0.042; i.e. liquefaction cannot be initiated, when the excess pore pressure has reached a value of about 4% of the in situ effective pressure.Maxima of the PWP ratio for each one of the investigated sites are given in Table 4.The soft soils overburden in the investigation area significantly amplified the input motion.The response spectra for all the sites are shown in Fig. 10.The sites varied in amplification from 1.5 (site PP) to 4.5 (site AA).The effect of soil amplification for this particular case is predominant as compared to the one of the soil liquefaction.Such amplification can be disastrous in the case of this particular scenario.

Scenario earthquake of the Iria fault (M=6.4)
The scenario earthquake on the Iria fault was judged from the equivalent-linear analysis as the most hazardous one as far as the probability of liquefaction is concerned.Values up to 80-100% were calculated at sites of AA, CC, FF , G1 and G2.
In its turn, the nonlinear analysis conducted employing D-MOD demonstrated that sites mentioned above will exhibit strongly nonlinear behaviour and at depths of 6-12 m the liquefaction initiation is quite high.Ground acceleration at the surface layer was computed to have similar or even smaller values compared to the input one.Additionally, the period of the ground acceleration is longer at the surface.This unequivocally suggests that soil liquefaction will be initiated to a certain extent at the aforementioned investigated sites during such an earthquake (Zeghal et al., 1994).When the excess pore pressures reaches 60-80% of the in situ effective pressure, then initial liquefaction will occur at the depth of 6-12 m.Site BB has better conditions than the previously examined, since the nonlinear analysis calculated PWP ratio does not exceed the value of 0.28.Detailed results of the analysis are presented in 13,15,17,19,21,23,25,27,29 depict the calculated acceleration and strain distributions at different depths at the investigated sites.As expected, the largest      14a, 16a, 18a, 20a, 22a, 24a, 26a, 28a, and 30a demonstrate the nonlinear hysteretic behaviour of the soil and weakening at the depths between 6-12 m at each examined site due to relatively high pore pressure generation.
The response of the sites with potentially liquefiable layers will be influenced by the build-up of excess pore water pressure during this scenario earthquake.Due to this fact, the characteristic feature is deamplification of surface spectral accelerations (Fig. 10a, c) with respect to the spectra for the bedrock inputs for all investigated sites Fig. 10b, d).

Scenario earthquake of the Xylokastro fault (M=6.7)
Although PL at the site CC reached 67% in this scenario, the nonlinear analysis calculated PWP lower than 0.02; therefore, liquefaction phenomena cannot be initiated.The maximum values of PWP ratio are included in Table 4 and the response spectra in Fig. 10.In this scenario the amplification factor was estimated to exhibit variations from site to site.The effect of soil amplification was the predominant effect in this case too.
Generally speaking, we can expect that the results obtained with aid of D-MOD are underpredicted most likely due to overdamping produced by poor modelling of soil behaviour due to difficulties in soil parameter evaluation.Therefore, it is reasonable to expect that additional geotechnical data would provide basis for more accurate analyses.

Conclusions and recommendations
In the area where the Nafplion city is under expansion, a layer with high probability to the liquefaction initiation was detected at shallow depths (6-12 m).Nonlinear analysis showed that it can substantially influence the ground motion during future earthquakes with similar characteristics to those of the examined scenarios.Depending on the scenario seismic source, effects of soil amplification and deamplification with consequent liquefaction initiation can be observed jointly or separately.
In particular, the scenario of Iria fault's activation by an M = 6.4 earthquake, is the most hazardous for the study area in terms of liquefaction initiation.In such case all examined sites can exhibit nonlinear behaviour mostly at depths of 6-12 m.Strong ground motion amplification by the shallow soft soil layer can be observed if the Epidaurus and Xylokastro fault rapture at lengths that correspond to M = 6.3 and M = 6.7 earthquakes, respectively.
From a practical point of view, the presence of the potentially liquefiable layer may cause several problems.For example, structure foundations close to the liquefiable layer may lose support, while ground settlement due to liquefaction may cause serious damage in the built environment.On the other hand, the lifelines and the existing railway lines may be destroyed in case of significant lateral displacements.
Since the area of expansion of the town of Nafplion was found to be susceptible to liquefaction initiation, improvement of the in-situ soils should be considered.This could be accomplished by proper modification of the soils to increase their stiffness and strength in response to static loads for many years.
Although there is a variety of improvement techniques (densification, reinforcement, grouting and mixing, drainage technique etc.) the use of densification and reinforcement methods may be more appropriate since they are suitable for the soil types characterizing the investigated area.

Fig. 1 .
Fig. 1.Geological map of the broader study area (onset: location of Nafplion on the general map of Greece).Locations of boreholes (red open squares) and seismic survey lines (black lines) are also mapped.

V
Fig. 2. (a) An example record of the WW survey.The record is dominated by the surface waves.(b) The fundamental Vphase(f) dispersion curve is well defined in the respective spectrum.(c) The resulted 1-D S-wave velocity model at WW site.

Fig. 5 .
Fig. 5.The statigraphic models at the sites QQ (a) and PP (b) as resulted from seismic refraction surveys of P-and S-waves.
Fig. 6.(a) Examples records of the seismic survey at AA and (b) EE sites.(c) The phase velocity spectrum with frequency gives a clear picture about the dispersion curve.

Fig. 8 .
Fig. 8. (a) An example dispersion curve from a Multichannel analysis of microtremors at the site G1.(b) The resulted 1-D velocity model of the passive survey was in a good agreement with the crosshole testing conducted nearby.(c) The 1-D Vs model at the site 3042.

Fig. 18 .Fig. 19 .Fig. 20 .Fig. 21 .Fig. 22 .
Fig. 18.Stress-strain history (a) for the layer with high probability to liquefy; PWP time histories at different ground level (b); PWP and shear strain time histories (c) for the layer with high probability to liquefy.

Fig. 23 .
Fig. 23.Acceleration (left column) and strain (right column) time histories at different ground level.

Fig. 24 .
Fig. 24.Stress-strain history (a) for the layer with high probability to liquefy; PWP time histories at different ground level (b; PWP and shear strain time histories (c) for the layer with high probability to liquefy.

Table 1 .
Parameters adopted in the stochastic simulation of strong ground motion from three scenario earthquakes in the broader area of Nafplion city (see text for further information on the selected scenario events).

Table 2 .
PGA values at the horizontal component of ground motion as computed stochastically for each one of the three examined earthquake scenarios and each examined site.

Table 3 .
Results of the liquefaction risk analysis at the examined sites.Depths of the water table -and of the potentially liquefiable layer are also included.

Table 4 .
Depth distribution of maximum PWP ratio for all investigated sites.