Abstract

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 two-dimensional 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.

1. 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 [18] 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 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 fortran-like 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 three-dimensional 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 half-space 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 two-dimensional setting. A homogeneous linearly elastic half-plane 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 first-order 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.

2. 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  kg/m2, Young modulus  N/m2, and Poisson ratio . These values characterize a real soil where caves are present [31]. In this medium, dilatational and shear waves propagate with velocities m/s and 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 , 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  m and  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  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 registered in the station of Nocera Umbra during the Umbria-Marche earthquake of 1997 is considered. The seismic acceleration is graphed in Figure 2(a), and it is taken from ITACA (Italian Accelerometric Archive), freely available in the website [33]. To study the wave scattering in the frequency domain, the frequency spectrum of the seismic acceleration is determined. From now on, a hat is used to indicate the Fourier transforms of a given time-depending function. Furthermore, indicate with the frequency and with the angular frequency.

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 with 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, . 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).

3. 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.

3.1. Time-Dependent Formulation

Let be the displacement of point of at the instant . The displacement field , with , must solve the problem 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 is normal to and points toward the inside of the cave. The soil viscosity is neglected. The dynamical solution 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). We estimate as sum of two terms: where and 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 is the stress tensor with being the identity tensor, being the gradient operator, and and being the symmetric part and trace operators, respectively. The field equation (5) admits both volumetric and shear waves. If we restrict to waves propagating in the direction , volumetric waves have the form , with being an arbitrary wave-shape and being the wave velocity (see (1)). Shear waves have the form , with being the wave velocity (see (1)) and 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 , 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 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.

3.2. 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 is the Fourier transform of . With the use of the Euler formula, (11) turns into From this relation, the Fourier transforms at two different vertical coordinates and are related as follows: If the frequency spectrum at the ground level is known (its modulus is reported in Figure 3), the displacement spectrum at a certain coordinate is and, from (8), the shear stress spectrum is Now the scattered problem (9) can be reformulated in the frequency domain as follows. Find which satisfies the equations where 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 complex-value solution of the form where and stand for real and imaginary parts. Thus the final solution 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 ; the solution accounts for a vertical oscillation as well and depends also on the coordinate .

4. 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 first-order 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 , sum of volumetric and shear waves, which hits the artificial boundary at a point with unit normal vector and unit tangential vector , satisfies the equation where the stress tensor 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), 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 , 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 m and 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 Hz. They belong to the frequency range 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 (see Figure 4(b)). It was found that, for the simulations with m, a reasonable radius choice is m, since the corresponding solutions are very close to those with the larger radius m. For larger values of , such as m, the radius length must be increased to 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 m and Hz and for different values of the radial coordinate (see Table 2(a)). Analogous curves for m are plotted in Figure 6.

From this analysis, one concludes that satisfactory accuracy is guaranteed if a value such as m is chosen. Furthermore, the mesh size m is assigned, which is enough fine to capture also high frequency oscillations (as an example, for Hz, shear waves have wavelength equal to 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 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 m to . 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 m. In conclusion, this energy analysis confirms the prescriptions on the choice of given by the above analysis on displacements and stresses.

5. Numerical Results

5.1. Simulation Setting

The semicircular domain of radius m, depicted in Figure 4(a), is discretized with plane quadrilateral finite elements of mesh size m. All cave geometries and positions listed in Table 1 are considered. Natural conditions are applied at the cave boundary, according to formulas (16) and (15). In (15), the amplitude spectrum of Figure 3 is assigned to . 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  Hz, by using a frequency step 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 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.

5.2. Scattered Solution

In Figures 9 and 10, the contour plots of , the real part of the horizontal scattered displacement, and , the amplitude of the horizontal scattered displacement, respectively, are represented for case 3A of Table 1 and frequencies 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 going away from the cave. Looking at the amplitudes at the ground level (Figures 10(b), 10(d), and 10(f)), it is observed that the values registered at point , right above the cave, are quite small in all the cases. This is due to the fact that the cave represents an obstacle which breaks the horizontal wavefront of the propagating seismic shear waves. The largest scattered displacements are registered at the points labeled with letters and in Figure 10, which are placed on the right and on the left of point , at a certain distance from it. These maximum amplitudes are , , and  m, for Hz, respectively. These reduce as the frequency increases.

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 m, labeled by letter A, the maximum amplitude is reached at Hz for the smallest cave 1A and at 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 Hz for cave 5A and at 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 (Figures 12(b), 12(d), and 12(f)) and (Figures 12(a), 12(c), and 12(e)) are plotted for frequencies Hz.

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 m are considered, and the mean value is estimated for each mesh size in the frequency range  Hz, in which the largest amplitudes are attained. The resulting graph of Figure 13 shows that the solution for m is practically equal to those for m, and, as a result, a satisfactory level of solution convergence is achieved for m.

5.3. 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 m and 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 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.

6. Concluding Remarks

In this paper, a numerical model to analyze the effects of underground cavities on the propagation of seismic shear waves has been proposed. It is based on the decomposition of the wave propagation problem into a free-field problem, whose solution is available analytically, and a scattered problem, which depends on the free-field solution, and can be solved by means of standard finite elements.

A first distinctive feature of the proposed model is its simplicity, which makes it easily implementable on commercial codes. For this reason, it is addressed to engineers and represents a tool for the design of new structures or for the seismic assessment of existing buildings placed on the surface of soils characterized by underground cavities. Indeed, it provides the frequency spectrum at the ground level, a crucial design datum in seismic structural engineering.

A second characteristic is its versatility. Here it has been applied to the very simple case of a single arch-shaped cave in a homogeneous soil. However, multiple caves of different shapes can be considered, and, nonhomogeneous stratified soils can be supposed, providing that a more complex solution of the free-field problem is considered. These further cases are left to future works.

For the single arch-shaped cave considered in this paper, an in-depth parametric analysis has been performed to determine the influence of the cave dimensions and depth on the frequency spectrum of seismic shear oscillations, in particular at the ground level. First, the scattered solution has been studied. It is characterized by transverse and longitudinal surface waves, which propagate on the ground surface from above the cave and by a reflected shear wave which propagates from the cave in the downward direction. Concerning the oscillation amplitude at the ground level we notice that the maximum values are not attained at point right above the cave, but in nearby points placed on the left and on the right of it, in agreement with analogous results found in [20, 26] by following different approaches. The positions of these latter points are not defined because they depend on the frequency of the incident wave and on the cave dimensions and depth. Thus we focus attention on the oscillation of point and on its relations with the cave dimensions and depth. We found that the cave amplifies the response frequency spectrum only within certain frequency range. Outside these ranges, the spectra exhibit smaller values than the corresponding free-field spectra, and thus the cave produces a shadow effect. The frequency intervals, where amplifications are registered and the extent of the amplifications is strictly related to both the dimensions and depth of the cave, are summarized in the following points.(1)For fixed cave depth, caves of “medium size” produces the maximum amplifications (cave 4 in the proposed simulations; see Figure 14). Indeed small caves have little effect, and large caves stand as an obstacle to the wave propagation.(2)For fixed cave depths, small caves amplify high frequency waves and large caves amplify low frequency waves. Indeed a relation of proportionality may be established between the wavelength of the amplified incident waves and the cave dimensions.(3)For fixed cave dimensions, the largest spectrum peaks are registered for shallow caves. As the cave depth increases, the amplified part of the spectrum reduces. This result was also found in [20].(4)For large shallow caves, amplification peaks are registered also at high frequencies because of reflected waves trapped in the soil portion between the cave and the ground surface (see Figure 14(f), cave 6A).

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This research was partially supported by the Italian Ministry of Education, Universities and Research (MIUR) by the PRIN funded Program 2010/2011 no. 2010MBJK5B.