Earth modeling and estimation of the local seismic ground motion due to site geology in complex volcanoclastic areas

Volcanic areas often show complex behaviour as far as seismic waves propagation and seismic motion at surface are concerned. In fact, the finite lateral extent of surface layers such as lava flows, blocks, differential welding and/or zeolitization within pyroclastic deposits, introduces in the propagation of seismic waves effects such as the generation of surface waves at the edge, resonance in lateral direction, diffractions and scattering of energy, which tend to modify the amplitude as well as the duration of the ground motion. The irregular topographic surface, typical of volcanic areas, also strongly influences the seismic site response. Despite this heterogeneity, it is unfortunately a common geophysical and engineering practice to evaluate even in volcanic environments the subsurface velocity field with monodimensional investigation method (i.e. geognostic soundings, refraction survey, down-hole, etc.) prior to the seismic site response computation which in a such cases is obviously also made with 1D algorithms. This approach often leads to highly inaccurate results. In this paper we use a different approach, i.e. a fully 2D P-wave «turning ray» tomographic survey followed by 2D seismic site response modeling. We report here the results of this approach in three sites located at short distance from Mt. Vesuvius and Campi Flegrei and characterized by overburdens constituted by volcanoclastic deposits with large lateral and vertical variations of their elastic properties. Comparison between 1D and 2D Dynamic Amplification Factor shows in all reported cases entirely different results, both in terms of peak period and spectral contents, as expected from the clear bidimensionality of the geological section. Therefore, these studies suggest evaluating carefully the subsoil geological structures in areas characterized by possible large lateral and vertical variations of the elastic properties in order to reach correct seismic site response curves to be used for engineering projects.


Introduction
Volcanic areas are often strongly heterogeneous.They may present large vertical and horizontal variations both in geometry and in their physical characteristics due to their peculiar sedimentation mechanism and to sin-depositional and post-depositional chemical processes.In particular, surge or pyroclastic flow deposits may show strong horizontal and vertical heterogeneities due to both variable granulometry and variable ambient deposition, temperature and chemistry.Variations of these parameters cause lithification processes to proceed differently even at a very small scale ( few or tens of meters).
Despite this heterogeneity, it is common practice, even in a geologically complex environment, like a volcanic area, to evaluate the seismic site response by means of the seismic Dynamic Amplification Factor (or DAF, see later) with monodimensional investigation methods.This approach may lead to highly inaccurate results.A correct (2D or even 3D) DAF estimation requires in fact an accurate reconstruction of the subsurface velocity field with a detail that is impossible to obtain in complex areas using conventional seismic methods.Such seismic surveyings (i.e.downhole, seismic refraction or even seismic reflection) are specifically directed towards obtaining layered models with sharp interfaces and constant material properties.On the other hand, pyroclastic deposits may show heterogeneities at all scales as demonstrated by geology or welllogs.Each 'layer' may contains numerous small scale variations of the material properties and these inhomogeneities are so irregularly distributed that they can hardly be imaged with traditional techniques.We have utilized a 2D tomographic technique to correctly define the subsurface structures.Then, by means of a 2D DAF methodology, the seismic site response was evaluated in such geological complex environments.

Brief review of pyroclastic deposits properties
A pyroclastic deposit is the result of the sedimentation of clastic material produced by one or more pyroclastic eruptions.Pyroclastic products can be classified according to the mechanism of emplacement in: a) pyroclastic flow deposits; b) surges deposits; c) fallout deposits.
Pyroclastic flows are heavier-than-air cloudflows that travel across the ground at velocities ranging from 10 m/s to 300 m/s.They can attain high temperatures and range from high density flows that move down valleys to dilute flows that extend over mountains (Fisher, 1966;Fisher and Smincke,1984).Pyroclastic flows are gravity flows and some properties (density and viscosity, for example) can change during transport and consequently influence the depositional properties.Pyroclastic surges are relatively low-concentration and highly expanded density currents with particles supported mainly by turbulence (Fisher, 1971) and particle concentration increasing within the lower part of the surge (Valentine, 1987).Pyroclastic surges may be frequently considered equal to «dilute pyroclastic flow».Lateral facies transitions in pyroclastic flow deposits show that pyroclastic flows and surges commonly evolve one from the other.Flow and surge, hence, may move as turbulent flows, in which the fragments with the higher settling velocities tend to migrate toward the base of the flow, gradually increasing fragment concentration and causing a vertical gradual variation in bulk density (Fisher, 1966(Fisher, , 1971(Fisher, , 1983)).Lateral variation may also be found for example in the ignimbritic deposits (i.e.rich pumice and glass pyroclastic flow deposits) depending on emplacement temperature; ignimbrites range in fact from unconsolidated to welded ignimbrites (Fisher and Smincke, 1984).
On the other hand, pyroclastic fallout consists of fragments that have been ejected from vents and have travelled through the atmosphere before falling (frequently almost cold) down to the earth.Due to their depositional mechanism and the consequent granulometric sorting with distance and homogenous deposition temperature, fallout deposits are less chaotic than surge and flow deposits at the scale of the engineering aimed seismic prospecting.
The interaction of depositional mechanism with topographic irregularities is also an important factor in the development of the pyroclastic properties.Whereas fall out falling down to the earth, drapes over the landscape and pyroclastic surges, due to their dilute characters can override the sides of a valley and may mantle topography similarly to fallout, flow deposits tend to fill depressions (Fisher, 1966).Therefore, topography strongly affects geometry and thickness of pyroclastic deposits.

Seismic imaging with traveltime tomography
Seismic refraction traveltime data for imaging the subsurface is a standard, easy and lowcost technique.However, conventional processing of refraction data generally uses simplified modeling.On the contrary, the information on the geometry and mechanical characteristics of the shallow subsurface is essential for seismic site response when in the presence of heterogeneous subsurfaces like pyroclastic flow and surge deposit.
Traveltime tomography may furnish a better detailed subsurface velocity model from picked traveltimes measured from seismic data associated with a variety of source-receiver configurations.Traveltimes are then used to infer the velocity distribution of the medium through which the seismic energy has propagated.
It should be noted that the «traveltime tomography» reconstructs earth velocity models with lower resolution than the other well know technique called «waveform tomography» (Lo and Inderwiesen, 1994), but, on the other hand, traveltime tomography is, in general, easier to apply and computationally faster.
The tomographic method employed in this study (see Stefani,1995) is the Turning Ray Tomography (TRT).TRT is based on first arrivals and this fact is appealing because first arrival times can be generally clearly identified and represent the best data available for near-surface velocity estimation.As said, TRT is useful for estimating near surface velocity structure in areas where conventional refraction techniques fail because of the lack of sharp velocity discontinuity and works well in those sites where the overburden is heterogeneous.The method utilizes first-arrival turning rays (continuously refracted direct rays) and head waves from conventional acquisition geometries to iteratively solve for P wave velocity in the near surface between sources and receivers.Therefore, TRT parametrization includes only velocity gradients without sharp surface discontinuities and utilizes an inversion algorithm based on the «backprojection method» (Langan et al., 1984;Stork and Clayton, 1991).
Traveltime tomography and particularly TRT, have been used mainly for near-surface velocity determination, e.g., for the purpose of «statics» solutions for seismic reflections (e.g., in the areas with large lateral velocity variations; see for example the paper by Bell et al. (1994) relative to the Mississippi Delta area characterized by mudlumps), with the aim of improving the stacking and migration of the deep seismic reflection image, to image salt domes or to check the area of dams (Yilmaz, 1988;Hale et al., 1992;Marsden, 1993;Frasheri et al., 1999).Other applications can be found in Zhu and McMechan (1989), Zhu et al. (2000).These studies suggest great caution when assuming a simple and regular stratification and homogeneity in the overburden as could be in a first approximation assumed for sake of simplification of the computation and mainly for surveying cost reduction.This approximation can only be assumed for a very limited range of divergence from the 1D modeling.

Dynamical Amplification Factor: monodimensional versus bidimensional evaluation
It is well known that the correct evaluation of the seismic site amplification is fundamental in seismic hazard areas (e.g., Italy and Japan), for the prevention of damage.Such evaluation may be complicated by the complexity of the subsoil which can strongly influence the path, the amplitude and the spectral content of the seismic waves at the surface.
Many efforts have been made up to now to derive algorithms that take into account all the parameters linked with the seismic wave propagation in the overburden volume of engineering interest (Aki et al., 1988;Bravo et al., 1988;Lee et al., 1989aLee et al., ,b, 1990;;Fishman and Ahmad, 1995;Rassem et al., 1995;Sherif et al., 1996;Bard, 1998;Field et al., 2000).
The ground motion (G (t )), can be expressed as the convolution operation between the overburden transfer function (T (t )) and the input motion at the basement (I (t )).T(t ) is related to the dynamical and geometrical parameters of the overburden while I(t) depends on the earthquake source, the path from the epicenter and the elastic parameters of all the crossed rock bodies.
For a continuous system G(t) is given by the following relation: where t is the time and τ is the lag time of the convolution operator.
In frequency terms, the relation (4.1) can be written as (4.2) where ω is the angular frequency.
To evaluate the seismic site response, we compute the Dynamic Amplification Factor (DAF) which represents the spectral ratio between the above said ground motion and input motion (4.3)When the overburden is homogeneous and geometrically regularly layered, the estimation of I(ω) is fairly simple and the DAF can be computed by means of 1D modeling (see, for example, Idriss and Seed, 1968;Rapolla and Mele, 1996).One of the well known computing techniques is based on a generalized, non linear, viscoelastic model (Idriss et al., 1992).This 1D method solves the propagation of harmonic vertical shear waves in a one-dimensional system in the frequency domain.The input physical model for the 1D code can consist of n infinite horizontal layers lying on the bedrock.Each layer, considered as homogeneous and isotropic, is characterised by given values of density (ρ), thickness, and shear modulus (G).
The equation utilized in 1D computing is the Kelvin-Voigt linearized visco-elastic equation where u is the displacement, ξ is the viscosity coefficient, t is the time and x is the position.
DAF computation becomes instead much more complicated when we deal with irregular structures and/or heterogeneous material.Many authors have studied the influence of irregular subsoil conditions and/or topography on seismic site response.Aki et al. (1988) saw that the finite lateral extent of surface layers and/or their irregular geometry can introduce additional effects such as resonance and focusing of the surface and body waves.Bravo et al. (1988) applied the boundary analysis method to study the response of a two dimensional horizontal stratified deposit of arbitrary shape under the incidence of SH waves.The response of a basin with two layers was analysed and the comparison with results from one-dimensional analysis shows significant differences.Lee et al. (1989aLee et al. ( ,b, 1990) ) published several studies on the scattering and diffraction of SV waves, plane SH waves, and P waves by circular cylindrical canyons with variable depth/width ratios.This study showed that focusing and resonance effects are greater on the valley flanks.
Recently, other studies have investigated problems of the Dynamic Amplification Factor (DAF) evaluation on geological structure with finite lateral extension.Fishman and Ahmad (1995) modelled alluvial valleys subjected to SH, SV and P waves, and demonstrated the influence on the seismic response of key parameters, such as valley depth, seismic impedance ratio, frequency and angle of incidence on surface ground motion.Rassem et al. (1995) reviewed the simple one-dimensional (1D) and the twodimensional (2D) finite element model concluding that the 1D model approach provides a poor approximation of the free-field motion in a valley with limited width, particularly near the valley edge.Sherif and Lee (1996) studied the diffraction due to plane SH waves around a circular alluvial valley in an elastic wedge-shaped medium.They demonstrated the influence on the surface displacement profile of the angle of the wedge, the frequency of the incident wave, the material properties of the media and the waves angle of incidence.Bruno et al. (1999) showed that in the case of a dipping bedrock if the dip angle is higher than about 10°, the results of DAF spectra computation by 1D and 2D algorithms start to diverge, both in terms of amplitude and peak periods.
2D modeling methods are based on finite differences (Vidale and Helmberger, 1988), or on boundary elements analysis (Kawase, 1988) or also on finite elements algorithms.As regards the latter (more used) methods, Beikae et al. (1994)  damping modulus (D).In addition, their code, incorporates a wave transmitting basement so that the half-space beneath the mesh can be modelled and therefore the need to assume a rigid boundary can be eliminated.
In order to compute the seismic site response, the 2D algorithm solves the following set of equations: where [M] is the mass matrix; [D] is the damping matrix; [K] is the stiffness matrix; u is the nodal displacement vector and I(t) is the earthquake load vector (input motion).
In this non linear system, once the input motion is known, the ground motion estimation is linked to the knowledge of the subsoil geometry and of the elastic parameters.Obviously, the more complex the subsoil is, the more detailed the knowledge of its properties needs to be.

Results
The 2D DAF evaluation approach described previously was applied to several sites in the Neapolitan area, Italy, covered by different pyroclastic deposits in the neighbourhood of Somma Vesuvius and Campi Flegrei volcanoes.The sites are located are within highly urbanized areas and are often close to high risk buildings (schools, city halls, hospitals, etc.).The V p tomographic models were already described by Bais et al. (2002).The other input parameters necessary for DAF estimation, specifically density, Poisson's ratio, and damping ratio were obtained from bibliographic data (Pellegrino, 1967;Nicotera and Lucini, 1967;Carrara et al., 1987;Guadagno et al., 1988;Rapolla and Mele, 1996) and are summarised in table I.In order to avoid lateral reflective effects during 2D modeling, we extended the length of the model of 50 m at both ends.The reference accelerogram for the sites close to Somma Vesuvius (fig.2) was calculated by deriving, with respect to the time, a digitally recorded Vesuvian earthquake of local magitude equal to 4.2.On the contrary, the reference accelerogram used for 1D and 2D DAF modeling analysis for the Campi Flegrei site (fig.3), is a synthetic input motion obtained from the Fourier antitransform of the spectral shape of the seismic motion at the rigid bedrock, given by Italian seismic law.The maximum peak value of the spectral acceleration was set according to Italian seismic law to be 0.04 g and in the Fourier domain the signal was recomposed by considering minimum phase shift for all frequencies of interest.
In this paper, only three typical examples are reported.The first site is Lagno Maddalena, near the town of S. Anastasia (fig.1).The tomographic model (fig.4a) shows an irregular morphology with a central low velocity structure that might probably reflect an erosion feature (fig.4a), possibly related to a paleochannel.Lagno Maddalena is indeed located over an alluvial fan derived from the reworking of pyroclastic soils.Cross sections cutting the fan (Cioni et al., 1999) show old anastomosing drainage patterns with- Table I.Average seismic of main lithotypes that characterize the subsoil of the Neapolitan volcanic areas.γ is the density, V p the p-wave velocity, V S the s-wave velocity, G the dynamic shear modulus, ν the Poisson's ratio, D the damping ratio (Pellegrino, 1967;Nicotera and Lucini, 1967;Carrara et al., 1987;Guadagno et al., 1988;Rapolla and Mele, 1996). in the lower pyroclastic units that were subsequently buried by more recent pyroclastic flows and by reworked volcanoclastic sediments.Figure 4a also shows the 2D finite element mesh modeling the subsoil of Lagno Maddalena.The model uses 235 nodal points and 199 cells.Figure 4b plots the 1D (black line) and 2D (red line) DAF (for both horizontal components) calculated at three surface nodal points.Comparison between 1D and 2D DAF shows, in all reported cases, different results both in terms of peak period and spectral contents, as was expected.Generally speaking, we can see that on 2D curves referred to all tested nodal points, the energy content shifts toward high frequencies with respect to the monodimensional curves.NS and EW DAF curves calculated at the nodal point #74, show a main peak at about 0.28 s on the 2D spectrum that is not visible on the 1D spectrum.
At the central nodal point #130, the main peak does not differ substantially from the 2D main peak.Both spectra have a main period at 0.18 s.This result is probably due to the fact that the nodal point #130 is in the middle of a symmetrical paleochannel.Symmetry possibly may cause reciprocal destructive interference of the 2D effects.At nodal point #155 the 1D and 2D DAF curves again become different.In fact, the 2D DAF curve has several peaks while the 1D DAF curve has only one peak at 0.1 s.Node 155 is located above the right flank of the paleochannel and the difference between the 1D and 2D response may again here be due to the fact that the monodimensional model does not take into account focusing effects.
The second example refers to Beneduce site, a place also located near Sant'Anastasia town (fig.1).Generally speaking, this site is characterized by a gradual lateral variation of the velocity field without strong contrasts (fig.5a).The final tomographic model (fig.5a) is characterized by a high velocity structure located between position 0 and 40 m. Figure 5a shows the 2D finite element mesh modeling which uses 154 nodal points and 130 cells.DAF computations were made on the right side of the TRT investigated section (fig.5d).DAF curves are shown at nodal points 21, 71 and 111.At nodal point #21 a peak at 0.26 s on the 2D spectra is not visible on the 1D curve.Focusing effects may be again invoked to explain this difference.At nodal point #71, the main peak of the 2D DAF curves is at 0.3 s while that of the 1D DAF curves is at 0.18 s.At nodal point 111, away from the flanks of the high velocity structure, the 1D and 2D DAF curves instead become almost similar.
The third example refers to a site near Giugliano (fig.1).The TRT model (fig.6a) shows an articulate bedrock morphology that is  probably due to erosion of the basement during the emplacement of the upper pyroclastic unit.Cole and Scarpati (1993) report that in this area during the deposition of Member B of the so called Neapolitan Yellow Tuff, representing the overburden of this site, the top floor of the lithic Member A (the basement) was strongly eroded.The 2D finite element mesh modeling is shown in fig.6a.The TRT output model was parametrized using 208 nodal points and 182 cells.DAFs were computed (fig.6d) at nodal points 29, 109 and 136.Comparison between 1D and 2D DAF modeling related to this site (fig.6b) shows significant differences at the nodal points 109 and 136, that is in correspondence of the maximum complexity of the geometry of bedrock high velocity structure while at the nodal point #29, where the subsurface morphology only slightly deviates from the monodimen-sional case, the 1D and 2D DAF curves are more similar.

Conclusions
The problem of the evaluation of real surface seismic motions by DAF analysis in volcanic areas characterized by a thick heterogeneous overburden have been discussed in the present paper.2D TRT surveys followed by a 2D non linear, viscoelastic DAF modeling allowed us to obtain a bidimensional seismic site response in many urban and suburban sites located at short distance from Mt. Vesuvius and Campi Flegrei.For all sites, DAF was also estimated with 1D non linear viscoelastic monodimensional algorithms.Comparison between 2D and 1D DAF estimations shows, as can be seen in all reported cases, results significantly different both in terms of peak period and spectral contents, as could have been expected from the clear bidimensionality of the geological profile as deduced from the 2D TRT survey.
As a matter of fact, many authors have shown that the simplification of the geological problems for DAF estimation may lead to unacceptable results.We have confirmed this statement showing that estimation of a 2D structure with 1D algorithms leads to gross mistakes.Obviously better results could be obtained by using 3D modeling, where the geometrical boundary conditions require it.
In volcanic areas, where large lithological heterogeneities may be present, it is of fundamental importance for engineering purposes to account for subsoil complexity.Our results suggest careful evaluation of the subsoil geological structures in these areas characterized by possible large lateral and vertical variations of the elastic properties, to reach correct seismic site response curves to be used for engineering projects.
solved the motion equations at the nodal points of a discrete grid being able to predict a strong generation of Rayleigh waves at dipping bedrock when SV waves are incident.Their 2D program is based on a dynamic, time domain equivalent linear algorithm and each nodal point is characterised by a value of the stiffness (K), density (ρ), thickness (h) shear modulus (

Fig. 1 .
Fig. 1.Map of the investigated Neapolitan urban and sub-urban areas.The study areas were located in Sant'Anastasia and Giugliano.

Fig. 3 .
Fig. 3.The input motion used for 1D and 2D DAF modeling in the test sites of Giugliano (below).Its spectral contents are reported in the upper part of the figure.It is a synthetic accelerogram obtained from the Fourier reconstruction of the spectral response of the Italian seismic low.