Effects of Underground Cavities on the Frequency Spectrum of Seismic Shear Waves

A numerical method is proposed to study the scattering of seismic shear waves induced by the presence of underground cavities in homogeneous soils. The method is based on the superposition of two solutions: the solution of the free-wave propagation problem in a uniform half-space, easily determined analytically, and the solution of the wave scattering problem due to the cave presence, evaluated numerically by means of an ad hoc code implemented by using the ANSYS Parametric Design Language. In the twodimensional setting, this technique is applied to the case of a single cave, placed at a certain depth from the ground level. The frequency spectrum of the seismic shear oscillation on the ground surface is determined for different dimensions and depths of the cave and compared with the spectrum registered without caves. The influence of the cave dimensions and depth on the spectrum amplification is analyzed and discussed.


Introduction
In the last decades a growing attention has been paid to the assessment of seismic vulnerability of existing ancient buildings and to the design of safe structures for seismic loadings.To this extent, a large variety of technical codes for designers and engineers have been developed, based on dynamical analysis or equivalent statical calculations.In all of them, crucial data to be taken into account is the frequency spectrum of the seismic oscillation registered at the ground level.The frequency content of a seismic signal depends on many factors such as the topography of the considered zone, the characteristics of the materials which the soil is made of, and the soil stratigraphy.These factors could cause amplifications and increase the seismic hazard for structures built on the nearby ground surface.
Some recent studies [1][2][3][4][5][6][7][8] have pointed out that natural or man-made underground cavities could also be a source of seismic hazards for structures.Earthquakes in a calcareous area carved by underground karstic process, like grottoes, caves, and dolines, can cause the falling in or collapse of the vaults in the former or the removal of debris in the latter.Surface depressions or subsidence of the ground can occur with the resulting danger of collapse for the building above [9].Some cases of serious damage were found in buildings over man-made cavities after the 1980-earthquake of Atella (Potenza, Italy) [1] and after the same earthquake, the most damaged buildings at Rionero in Vulture were mainly concentrated in areas where there was notable cave density [2].Gizzi [3] highlight the relationship between the presence of underground cavities and the observed building damage after the 1930-Irpinia earthquake.Similarly, in [4] it was noted that cavities can significantly increase the damage to nearby buildings in the presence of low strength rock.It also seems that the caves can reduce the seismic loading transferred to buildings.Indeed, a small reduction of the peak ground acceleration has been surveyed in areas in which the spatial distribution of the cavities is of some hundred meters [5].Man-made caves and their relationships with historical buildings were considered in defining their risk level [6].In the city of Catania seismic response analysis has been performed in proximity of several natural cavities to find that cavities represent a high risk for foundation stability of some buildings [7].Often ancient buildings are on the surface of soils where many subsurface cavities are present, natural or man-made, excavated over the centuries for housing, defensive purposes, or ritual and religious intents.Italy, for instance, is rich of underground 2 Advances in Civil Engineering cavities below ancient medieval towns, castles, monasteries, and so forth.These towns are built on essentially man-made underground cavities dating back to many centuries, which are Roman underground water systems of aqueducts and cisterns, medieval tunnels for storage and communications, and hiding places used during the Second World War [8].
To perform a correct seismic structural analyses, a crucial point is the determination of the effective ground surface acceleration due to earthquakes, and, consequently, the finding of the seismic loadings on the buildings.If caves are present underground, then their scattering effects on the seismic wave propagation must be considered for an accurate estimate of seismic accelerations at the ground level.
A first objective of the present study is to propose a numerical technique for determining the seismic frequency spectrum at the ground surface when subsurface caves are present in the soil, for a given frequency spectrum registered at the surface of a soil without caves.To implement this, the ANSYS Parametric Design Language [10] (a fortranlike programming language) was used.An ad hoc code was developed as a computational tool to assist engineers and designers.Indeed, commercial computer software dealing with wave scattering in continuum medium with voids is few or is not comprehensive.On one hand, software for structural analysis may handle sophisticated material nonlinearities and complex geometries, but often it is inadequate for problems of wave diffraction and radiation.On the other hand, codes dealing with seismic wave propagation in soils mainly concentrate on the effects of multilayered stratigraphy.One of these latter codes is the computer program SHAKE [11,12], which estimates the wave reflection and diffraction occurring at the interfaces between soil layers made of different materials.
Once the computational model was developed, a second objective of the present work is to estimate the influence of the dimensions and depth of an underground cave on the amplitude of seismic wave registered at the ground level, by performing numerical parametrical analysis.
Researches on seismic wave scattering produced by cavities mainly focus on two aspects of engineering interest.A first set of works concentrates on determining the dynamical response of subsurface caves, tunnels in particular, to seismic waves.Different approaches have been followed.In [13,14], for instance, models based on boundary elements have been developed, while, in [15], both boundary elements and finite elements have been used.In [16], the earthquake-induced strains on tunnels have been evaluated by means of a threedimensional shell theory, in [17] a multiscale finite element approach has been followed, and, in [18], simplified solutions have been found by modeling tunnels as Timoshenko beams.A second set of studies focuses on the diffracting effects of caves on the wave motion, paying a special attention to the resulting oscillations at the ground level.The present study belongs to this set.A variety of analytical solutions have been found.Only a few papers are mentioned, and the works therein quoted are referred to for more bibliographic references.A review of closed-form solutions for the scattering of pressure and shear waves by a single spherical obstacle has been proposed in [19].The diffraction problem for shear waves hitting a circular cylindrical cavity placed in an halfspace has been solved in [20,21], and the seismic amplification of surface ground motion above the cavity has been estimated.While in [20] the diffracted solutions of SV-waves have been found by means of Fourier-Bessel series, in [21] the antiplane problem related to the propagation of SH-waves has been solved by using Bessel and Hankel functions.Approximated solutions are determined in [22] by using a hybrid approach which combines the finite element method in the region around the cavity and wave eigenfunction expansions far away from the cavity, in [23] by applying the weighted residual method and in [24] by finite elements.Amplification of the seismic risk for surface structures on soils with voids has been investigated in [25].In it, the response of a rigid embedded foundation to antiplane SH-waves diffracted by a cylindrical cavity has been evaluated by using Bessel functions.In [26], the influence of SV-waves in the surface motion has been investigated numerically by means of the explicit finite difference program FLAC [27].The diffraction of waves by cavities has been studied in the real case of the town of Castelnuovo (Italy), characterized by many underground cavities, which was stricken by the violent L' Aquila earthquake in 2009.
Here we consider only SV-waves, with plane wavefront and upward direction of propagation, which are the most demanding for ancient buildings since they transmit inertial horizontal forces.The model is restricted to the twodimensional setting.A homogeneous linearly elastic halfplane with a single void placed at a certain depth from the ground is considered.Since the problem is assumed linear, the strategy we follow to solve it is based on the superposition of two solutions.
(i) The first one is the so-called free-field problem, that is, the problem of wave propagation in a homogenous half-plane.For this, the solution is available analytically.
(ii) The second one, named diffracted problem, accounts for the perturbation brought by the presence of caves to the wave propagation.This problem consists in determining the motion due to a certain distribution of forces applied on the cave boundaries, determined as function of the free-field solution.Since this problem cannot be solved analytically, it is approached numerically, by means of finite elements.A semicircular computational domain is considered, and absorbing conditions are assigned at the artificial boundary to avoid wave reflection.The two main techniques for absorbing outgoing waves are the high-order absorbing boundary conditions and the perfectly matched layer (see [28,29] for a comparison of these methods and the references therein quoted for detailed descriptions of them).Here, the first-order absorbing boundary conditions are considered [30], avoiding more complex higher-order conditions.First-order conditions perfectly filter waves which hit the boundary with null incident angle.Since in the diffracted problem waves generated at the cave boundary propagate with circular wave fronts, hitting the artificial boundary with very small incident angles, then firstorder conditions guarantee satisfactory accuracy, providing that a sufficiently large radius of the semicircular domain is considered.Once the diffracted problem is numerically solved, the solution of the global problem is obtained as sum of the diffracted and free-field solutions.
In this paper, the above described strategy is applied to the problem in the frequency domain.The frequency spectrum of seismic shear oscillations is determined at the ground level, as required for engineering applications where the frequency content of shear oscillations is a fundamental input data for seismic analyses of existing or new buildings.Several simulations are performed by considering different dimensions and depth of an arch-shaped cave.For each simulation, the seismic oscillation registered during the Umbria-Marche earthquake of 1997 is assumed as the solution of the free-field problem.The amplitude of its frequency spectrum is evaluated, and then the corresponding diffracted problem is solved by using finite elements.Finally, the global solution is estimated, and the amplitudes estimated at the ground level, right above the cavity, are compared with those of the free-field case.Depending on both dimensions and depth of the cave, frequency ranges are found in which oscillation amplitudes are amplified and frequency intervals are found in which they are reduced.A detailed description of the found spectra is made, and the differences with the free-field spectrum are discussed.
The paper is organized as follows.In Section 2, the problem is formulated.The material characteristics of soil and the cave geometry are defined, and the free-field solution is determined.In Section 3, the strategy used to solve the problem is described first in the time domain and then in the frequency domain.Section 4 deals with the absorbing boundary conditions assigned at the artificial boundaries of the computational domain.The radius of the semicircular computational domain is determined by means of numerical analyses.In Section 5 the numerical results are presented and discussed.Concluding remarks are drawn in Section 6.

Problem Statement
The scattering problem of seismic waves caused by underground cavities in soils is formulated in the two-dimensional setting.The soil is represented by a half-plane, homogeneous and linearly elastic, with mass density  = 2000 kg/m 2 , Young modulus  = 21.333⋅ 10 7 N/m 2 , and Poisson ratio ] = 0.33.These values characterize a real soil where caves are present [31].In this medium, dilatational and shear waves propagate with velocities  V = 400 m/s and   = 200 m/s, respectively.They are related to  and ] by It is assumed that a cave is placed in the soil at a certain depth ℎ from the ground level and has the shape drawn in Figure 1.
Then, it is suppose that seismic shear waves propagate in the soil from the bottom to the top, with wave velocity   .Figure 1 sketches a geometrical scheme of the problem.To represent points and vector by components, the orthogonal Cartesian frame (; e 1 , e 2 ), depicted in Figure 1, is used.The aim is to investigate the scattering produced by the cave on the propagation of shear seismic waves.In particular, the seismic oscillations registered at the ground level are evaluated, right above the cave (point  in Figure 1), and the differences with the oscillations found when cavities are absent are registered.Of interest is the analysis of the influence of the cave dimensions and of its depth on the seismic oscillation at the ground level.For a cave of the shape drawn in Figure 1, consider six different dimensions and four different depths.The smaller cave has dimensions  = 2 m and  = 2.5 m (cave width and height, resp., as shown in Figure 1), and caves of dimensions 2, 3, 4, 6, and 8 times larger are considered.The four different depths are ℎ = 1, 2.5, 5, 10 m.All the possible combinations listed in Table 1 are considered.Numbers are used to label the different cave dimensions and capital letters to distinguish the depths.Sizes 1 and 2 are representative of caves found under the typical medieval centers of Marche.Dimensions 3, 4, 5, and 6 are introduced to carry out a parametric analysis and are typical of anthropogenic cavities and natural caves found in other sites of Italy [1,4,6,32].As the seismic excitation is concerned, the north-south component of the accelerogram Due to the random nature of the seismic oscillation, its spectrum is determined by using stochastic concepts, as described in [34].The whole seismic time interval has been divided into subintervals, each one overlapping the next one for 1/3 of the subinterval length.In this way a sample of the random process is generated.The subinterval length has been set to get the necessary resolution of 0.5 Hz in generating the FFT of the acceleration.At each subinterval the FFTs of the signal are determined.Once all the resulting amplitude spectra are obtained, the data sample is generated, and its statistical characteristic minimum, maximum, and mean curves are estimated (dotted curves in Figure 2(b)).Then, these curves were approximated by means of the lognormal probability density function, as described in [35].The lognormal function has the expression

Advances in Civil Engineering
with E being the total energy of the signal,  being the mean frequency of the signal, and  being the root mean square of the signal.For , , and  evaluated from the minimum, the maximum, and the mean Fourier transform curves, the resulting lognormal curves are those plotted with solid line in Figure 2(b).We consider the average curve.From it, the frequency spectrum |û|() of the horizontal displacement amplitude is obtained by simply dividing â() by the square of ; that is, |û|() = |â|()/ 2 .The graph of |û| is plotted in Figure 3.The problem is stated as follows.
Let a linearly elastic, isotropic half-space, representing the soil, with a subsurface cavity of the shape of Figure 1 be given.Let the soil be subject to a seismic shear oscillation, characterized by the frequency spectrum of Figure 3, registered at the ground level in the case of an uniform, free-of-cavities soil.Determine the wave scattering produced by the cave, and, in particular, the oscillation frequency spectrum registered at ground level right above the cavity (point  in Figure 1).

Solving Strategy
The strategy proposed in [36], Section 6.1, for the problem of seismic waves on a soil with an excavated foundation is adapted to the problem at hand.
where the first equation is the governing equation of linear elastodynamics, and the second and third equations are the stress-free boundary conditions for the ground and the cave surface, respectively.The unit vector n is normal to Ω  and points toward the inside of the cave.The soil viscosity is neglected.The dynamical solution u of the above problem would sum to the solution of the static problem, where all the static loads are considered.For this reason volume loads, like the soil weight, are not taken into account in (3) 1 .We estimate u as sum of two terms: where u  and u  are the solutions of the free-field and scattered problems, respectively, described in the following.
Free-Field Problem.The free-field problem is the problem of seismic wave propagation in the soil without caves.It is governed by the following equation of motion and stress-free condition at the ground level: where T is the stress tensor (tr ∇u) I) , (6) with I being the identity tensor, ∇ being the gradient operator, and sym(⋅) and tr(⋅) being the symmetric part and trace operators, respectively.The field equation (5) 1 admits both volumetric and shear waves.If we restrict to waves propagating in the direction e 2 , volumetric waves have the form ℎ( V  − e 2 ⋅ x)e 2 , with ℎ being an arbitrary wave-shape and  V being the wave velocity (see (1) 1 ).Shear waves have the form ℎ(  −e 2 ⋅x)e 1 , with   being the wave velocity (see (1) 2 ) and e 1 being the direction of oscillation.Volumetric and shear waves are also called longitudinal and transverse waves, since their propagation and oscillation directions are mutually parallel and orthogonal, respectively.As justified above, only shear waves are considered.
To satisfy the stress-free conditions at the ground level, a reflected wave must be accounted for.Thus the free-field solution is The corresponding stress state is of pure shear, with T  =   (e 1 ⊗ e 2 + e 2 ⊗ e 1 ), and In the expression above, ℎ  () = ℎ()/.
Scattered Problem.The scattered solution represents the perturbation to the free-field solution due to the presence of the cave.The scattered problem is obtained by applying the decomposition (4) to the global problem (3) and using ( 5) of the free-field problem.It reads as div The scattered problem may be solved once the free-field solution is found, since, in it, the inverses of the forces induced by the free-field oscillations are applied on the boundary Ω  of the cave, representing the force source for the scattered problem.
It is assumed that the boundary of the cave is stress-free.This assumption corresponds to a cavity with no structural covering on its surface.An opposite situation could be that of a cavity covered by a rigid structure.In this case the cave would be modeled as a rigid inclusion.

Formulation in the Frequency Domain.
In order to determine the frequency spectrum at ground level, the problem is formulated in the frequency domain.
For the free-field problem, use the Fourier transform of the transverse oscillation Apply rescaling and, then, the translation rules to obtain where ĥ() = ∫ ∞ −∞ ℎ() exp(−2) is the Fourier transform of ℎ.With the use of the Euler formula, (11)  From this relation, the Fourier transforms at two different vertical coordinates  2 and  2 are related as follows: If the frequency spectrum û (0, ) at the ground level is known (its modulus is reported in Figure 3), the displacement spectrum at a certain coordinate  2 is and, from (8), the shear stress spectrum is Now the scattered problem (9) can be reformulated in the frequency domain as follows.Find û = û (x, ) which satisfies the equations div T (û  ) +  2 û = 0, in Ω, where T has expression (6) and τ is given by (15).Since the free-field solution is known from formula ( 14), only the scattered problem ( 16) must be solved.The solution is obtained numerically by means of the finite element method.A semicircular computational domain Ω  of radius  is considered, as shown in Figure 4(a), and absorbing boundary conditions are assigned on the semicircle Ω  .The assignment of the filtering conditions and the choice of the radius  are discussed in the next section.
In general, the scattered problem ( 16) admits a complexvalue solution of the form where Re{⋅} and Im{⋅} stand for real and imaginary parts.Thus the final solution û(x, ) is the sum of the free-field and scattered solutions: while the free-field solution has only the horizontal component, and it does not depend on  1 ; the solution accounts for a vertical oscillation as well and depends also on the coordinate  1 .

Absorbing Boundary Conditions
To solve the scattered problem by finite elements, a computational finite domain must be considered in place of the semi-infinite half-plane, and absorbing boundary conditions must be applied on the artificial boundaries of the computational domain.They should filter outgoing waves, without reflecting them, thus permitting a reprodution of the wave motion observable in the semi-infinite domain.The firstorder absorbing conditions proposed by Lysmer and Kuhlemeyer [30] are implemented.They consist of a distribution of normal and tangential dashpots with viscous coefficients equal to respectively, on the artificial boundary.These conditions perfectly filter volumetric and shear waves propagating in the direction normal to the boundary.Indeed, a wave of the form u(x, ) = ℎ V (x ⋅ n −  V )n + ℎ  (x ⋅ n −   )t, sum of volumetric and shear waves, which hits the artificial boundary at a point with unit normal vector n and unit tangential vector t, satisfies the equation  where the stress tensor T is given by ( 6).Thus special care must be devoted to the choice of the artificial boundary, to maximize the absorption.It should be such that outgoing waves hit it with an incident angle as close to zero.In this case, the wave source is represented by the cave boundary, where forces are applied according to (16) 3 , and we expect that, at a certain distance from the cave, the wave fronts of both shear and volumetric waves are practically semicircles centered at the cave, which propagate in the radial directions.Thus we assign to the computational domain the shape of a semicircle, as depicted in Figure 4(a).We call it Ω  and the artificial boundary Ω  .
Given the cave dimension  = max(, ), the cave depth ℎ, and the radius  of the semicircular domain, the geometrical conditions  ≫  and ℎ ≪  must be satisfied in such a way that outgoing waves arrive at Ω  with almost semicircular fronts.They require that the cave must be very distant from Ω  .It is required that the cave dimension and its depth must be small compared to the wavelength ; that is,  ≪  (hypothesis of compact source) and ℎ ≪  (hypothesis of shallow depth of the cave).These latter conditions allow the use of coarse meshes, in comparison to the geometrical dimensions  and ℎ, since the mesh size must be at least six-eight times smaller than the wavelength for an accurate numerical approximation.
Choice of the Radius .To fix the value of  the following numerical tests are performed.Consider the semicircular domain represented in Figure 4(b) and apply a harmonic horizontal force of amplitude equal to 1 N at the point of the vertical radius placed at depth ℎ.In the frequency domain, numerical solutions are determined for different radii  and frequencies .The radius  is selected in such a way that the corresponding solution does not differ significantly from those obtained with larger radii, whatever the frequency  is.
For the two depths ℎ = 1 m and ℎ = 30 m of the force, consider the radii  and the mesh sizes ℎ  reported in Table 2.For each case, four simulations are performed, by assigning the frequencies  = 2, 3, 4, 5 Hz.They belong to the frequency range (1,6) Hz in which the frequency spectra of common seismic waves attain the largest amplitudes.
The solutions of all these simulations have been analyzed by comparing the radial and tangential displacements û and û and the radial and tangential stresses σ and σ along a radius inclined of an angle of /4 (see Figure 4(b)).It was found that, for the simulations with ℎ = 1 m, a reasonable radius choice is  = 50 m, since the corresponding solutions are very close to those with the larger radius  = 100 m.For larger values of ℎ, such as ℎ = 30 m, the radius length must be increased to  = 150 m to obtain accurate results.All the comparative analysis cannot be reported here to limit the length of the paper.In Figure 5, the real parts of û , û , σ , and σ are reported, in the case of ℎ = 1 m and  = 4 Hz and for different values of the radial coordinate (see Table 2(a)).Analogous curves for ℎ = 30 m are plotted in Figure 6.
From this analysis, one concludes that satisfactory accuracy is guaranteed if a value such as  = 150 m is chosen.Furthermore, the mesh size ℎ  = 1 m is assigned, which is enough fine to capture also high frequency oscillations (as an example, for  = 20 Hz, shear waves have wavelength equal to 10 m).The mesh in the case of cave 3C of Table 1 is plotted in Figure 7.
Remark 1 (energy on the artificial boundary).The accuracy of the numerical solutions for different values of  can be evaluated also on the basis of energetic estimates, as described in the following.The energy dissipated by the dashpots on the artificial boundary Ω  in a period  = 1/ is determined for different values of .Then a solution obtained for a certain radius  is considered accurate if the dissipated energy on Ω  practically does not change for larger radii.In other words, the right radius is that at which the dissipated energy is practically stabilized on a constant value.
The energy radiated from the computational domain in a period is given by the expression where are the radial and tangential displacement components, respectively, with û and û being numerically determined.For simulations of type (a) in Table 2, the values of  as function of the four considered radii are plotted in Figure 8 for different frequencies .It can be seen that  changes very slightly passing from  = 50 m to  = 100.For simulations of type (b) in Table 2, the curves are analogous to those of Figure 8 and thus are not reported.The only difference is that, in this case, curves stabilize (i.e., become practically constant) for  > 100 m.In conclusion, this energy analysis confirms the prescriptions on the choice of  given by the above analysis on displacements and stresses.

Simulation Setting.
The semicircular domain of radius  = 150 m, depicted in Figure 4(a), is discretized with plane quadrilateral finite elements of mesh size ℎ  = 1 m.All cave geometries and positions listed in Table 1 are considered.Natural conditions are applied at the cave boundary, according to formulas (16) 3 and (15).In (15), the amplitude spectrum |û| of Figure 3 is assigned to û (0, ).Thus, at each frequency, we determine only the amplitude of the oscillation, losing information on the phase.It follows that the found frequency spectrum represents an upper-bound of the real frequency spectrum, which would be obtained by considering in (15) the complex-valued spectrum û of the oscillation registered at the ground surface.For each geometry listed in Table 1, frequencies are selected within the interval (0, 20) Hz, by using a frequency step Δ = 0.2 Hz.The absorbing conditions described in Section 4 are applied at the artificial boundary Ω  , while stress-free conditions are assigned on the ground line.
The scattered problems considered here are antisymmetric with respect to the  2 axis (Figure 4(a)).Thus they could be solved by considering one-half of the domain and imposing null vertical displacements at the symmetry boundary, with the advantage of reducing the computational effort.Since the objective of this work is to propose a computational tool available for nonsymmetric problems, the whole computational domain has been considered.

Scattered Solution.
In Figures 9 and 10, the contour plots of Re{û  1 }, the real part of the horizontal scattered displacement, and |û  1 |, the amplitude of the horizontal scattered displacement, respectively, are represented for case 3A of Table 1 and frequencies  = 4, 6, 8 Hz.Analogous trends are found in simulations with different frequencies, dimensions, and depths of the cave, not reported for the sake of conciseness.Figure 9 shows that surface shear waves, propagating at the ground level, are generated.Their wavelengths are inversely proportional to the frequency, consistent with the formula  =   /. Figure 10 points out that their amplitudes decrease   Concerning the oscillation of point , the frequency spectra of the horizontal oscillation amplitudes for all the cave dimensions and depths of Table 1 are plotted in Figure 11.As expected, the amplitudes increase as the cave dimensions increase and the cave depth reduces.It is noticed that frequencies at which the curves reach the maximum values reduce when the cave dimension increases.For caves placed at depth ℎ = 1 m, labeled by letter A, the maximum amplitude is reached at  = 5.4 Hz for the smallest cave 1A and at  = 2.8 Hz for the largest cave 6A.Since the frequency is inversely proportional to the wavelength, a relation of direct proportionality establishes between cave dimensions and wavelength of the amplified waves.It is related to two phenomena which depend on the ratio between wavelength and cave dimension: from the one side, waves with large wavelengths in comparison to the cave dimensions (small frequencies) propagate practically undisturbed, since the cave represents a small perturbation; from the other side, for waves with small wavelengths compared to the cave dimensions (large frequencies), the cave represents an obstacle which obstructs the wave propagation.
For caves of large dimensions (labels 5 and 6), secondary high frequency peaks are observed.For instance, a secondary peak presents at  = 9.0 Hz for cave 5A and at  = 8.8 Hz for cave 6A.They could be due to reflected waves trapped in between the cave and the ground surface.
The approximation of the unbounded half-space by a discretized finite mesh leads to the presence of few poles close to the real axis of the complex plane; in some circumstances this can be highlighted by the presence in the frequency response of oscillations, as can be seen in Figure 11.
A second effect of the presence of the cave is the generation of reflected waves which propagates in the downward direction, with circular wavefronts.This wave motion localizes on a cone centered in the cave, which is indicated by dashed lines in Figures 9 and 10.
The scattered problem also accounts for vertical displacements, as shown in Figure 12, where the contours of Re{û Summing up, the ground surface is subject to two types of progressive waves, a longitudinal wave characterized by horizontal oscillations (Figure 9) and a transverse wave with vertical oscillations (Figure 12).In both cases the wave speed is   and the wavelength is  =   /.
A posteriori, we can affirm that the semicircular computational domain is well-chosen.Indeed the wavefronts of the scattered solutions are semicircular and overlap to the artificial boundary, guaranteeing a very good absorption by the dashpots.
The accuracy of the proposed analysis has been checked by comparing the results with those given by simulations with different mesh sizes.It was found that the solution does not change significantly for smaller mesh sizes, and this confirms the accuracy of the proposed results.In Figure 13, the results of one convergence test are shown among the several checks which have been done.In the case of cave 3C, the four mesh sizes ℎ  = 3, 1.5, 1, 0.5 m are considered, and the mean value ) is estimated for each mesh size in the frequency range ( 1 = 0,  2 = 10) Hz, in which the largest amplitudes |û  1 | are attained.The resulting graph of Figure 13 shows that the solution for ℎ  = 1 m is practically equal to those for ℎ  = 0.5 m, and, as a result, a satisfactory level of solution convergence is achieved for ℎ  = 1 m.

Global Solution.
The frequency spectra of the global solution at point , evaluated by means of formula (18), is plotted in Figure 14.For cave 1 the spectrum is practically similar to that of the free-field case.Indeed the cave size is too small to modify the wave propagation significantly.This is in agreement with theoretical [37] and numerical [26] results, according to which the seismic wave propagation is modified only if the cave dimensions are larger than a quarter of the wavelength of the incident wave.Since caves of type 1 have dimensions  = 2 m and  = 2.5 m, they can only affect the propagation of shear waves with frequency larger than 20 Hz.For larger dimensions, the spectra deviate from the  free-field spectrum.It exhibits larger values than the free-field spectrum in certain frequency intervals and smaller values in other frequency ranges.These frequency intervals and the corresponding oscillation amplitudes vary with the cave size and depth as described in the following section.First consider the caves of type A (depth ℎ = 1 m).Their spectra are the blue curves in Figure 14.A frequency region is indicated by   , in which the curves exhibit larger values than the free-field curves.As the size of the cave increases, the region   shifts toward lower frequencies and reduces in amplitude.The frequency intervals   and the corresponding wavelength intervals   are listed in Table 3, for different dimensions of the cave.For increasing cave dimensions, a shift of   toward lower frequencies corresponds to a shift of the wavelength interval   toward larger wavelength.As expected, a relation of proportionality is established between the dimensions of the scattering object and the wavelength of the scattered waves.For frequencies  smaller than those of   , the cave is so small in comparison to the wavelength of the incident wave that it is unable to scatter the wave.For frequencies larger than those of   , the cave is so large compared with the wavelength of the incident wave that it behaves as an obstacle, deviating the wave to left and right sides.As a result, at point , placed right above the cave, the spectrum amplitude reduces.However, for sufficiently large caves, such as caves 4A, 5A, and 6A, when the frequency exceeds a certain frequency value, indicated with  in Figures 14(d)-14(f), the spectrum curve became again larger than the free-field curve.This behaviour is likely due to waves reflection which generates in the region between the cave and the ground, once the vertical direction of wave propagation is deviated by the cave.For increasing cave dimensions, the threshold value  reduces.
The maximum amplitude is registered for the "medium size" cave 4A.Indeed, small caves are ineffective, while large waves represent obstacles for wave propagation, producing a reduction of the spectrum.
Up to now, only the spectra corresponding to caves of type A have been described.However analogous considerations can be made for the spectra related to caves of types B (red line), C (violet line), and D (green line).The only difference is that as the cave depth grows, the spectrum diagrams and, in particular, the peaks reduce in amplitude.It follows that the spectra of deep large caves, caves 5D and 6D, are always below the free-field spectrum.Thus, in these cases, the cave produces a beneficial effect, whatever the frequency value is, as reported in [4,5].
From the previous results some clear trends emerge, and so some remarks can be made for engineering purposes.Figure 15 synthetically shows the peak displacement of the control point  for the different analyzed cases.It underlines that underground cavities could play a role in modifying the surface motion.In particular, when cavities are not far from the ground, the displacement could also increase of about 35% when cavities are present.In the most unfavorable case, a similar increase can be estimated for the displacements of an upper structure, if present.This surely has to be considered for the design of new structures or the seismic safety evaluation of existing structures, to resist future earthquakes.

Figure 3 :
Figure 3: Frequency spectrum of the shear oscillation at the ground level.

Figure 4 :
Figure 4: (a) Computational domain.(b) Problem solved to fix the radius .

10Figure 9 :
Figure 9: Real part of the horizontal scattered displacement, Re{û  1 }, for cave 3A at different values of the frequency .

Figure 10 :
Figure 10: Absolute value of the horizontal scattered displacement, |û  1 |, for cave 3A at different values of the frequency .

Table 1 :
Different dimensions [m] and depths [m] of the cave, labeled by numbers and capital letters, respectively.

Table 2 :
Geometry values used in the simulations performed to fix .

Table 3 :
Frequency intervals   and corresponding wavelength intervals   for different dimensions of the cave.