1 Introduction

The frequency and spatial variations introduced by local topography on seismic ground motions are known to produce concentrated damage during earthquakes (e.g., Northridge, California, 1994; Kobe, Japan, 1997; Kocaeli, Turkey, 1999). In the theory of elastodynamics the quantification of these effects implies the solution of a wave scattering problem. Modern numerical techniques and the available computer power, has led to the possibility of realistic simulations of these effects for large scale regional domains, e.g., (Bao et al. 1998; Okamoto et al. 2011; Taborda and Bielak 2011; Cupillard et al. 2012; Restrepo et al. 2012). Solution methods with improved versions of the boundary element technique have also been proposed to find the scattered field by irregular topographies (Zhou and Chen 2006, 2008; Wu and Maupin 2007; Liang and Liu 2009). At the same time, actual field records of large urban areas have become increasingly available during the last years providing additional evidence of the effects of topography. Despite these strong advances on those fields, a detailed understanding of the underlying physics of the problem is still required. In this work we study the effects of topography on seismic ground motions, using a partition of the field into physically meaningful terms based upon the concept of diffracted motions.

In earthquake engineering the terms scattering and diffraction have been used loosely as synonyms, however they are not (Sánchez-Sesma and Iturrarán-Viveros 2001; Mow and Pao 1971). A strict definition of scattered field has been first attributed to Rayleigh: "A scattered wave is the difference of the total wave field observed in the presence of an obstacle and the incident wave". In the case of an obstacle supported by a half-space, it is the difference of the total wave and the free field, Pao and Varatharajulu (1976). On the other hand, according to Keller (1962), diffracted waves are produced by incident waves which hit edges, corners or vertices of boundary surfaces or which graze such surfaces. This connection between the diffracted field and the geometric entities in the scatterer is expected to be reflected in the modifications to the ground motions imparted by the topographic effects.

The theory of diffraction has a long history. Its physical aspects have been studied in great detail, mainly in the context of electromagnetic waves. A landmark contribution and converted thereby in one of the building blocks for advancing the theory, is identified in the work of Sommerfeld (1896). He found the complete solution for the diffraction of electromagnetic plane waves by a semi-infinite crack in a homogeneous medium. Shortly after Sommerfeld’s work, MacDonald (1902) delivered the solution for the total field on a wedge, under plane and cylindrical SH waves. He wrote the total field as series expansions in terms of Bessel functions. A detailed study of Macdonald’s solution, highlighting various aspects of the diffracted field can be found in Sanchez-Sesma (1985).

A major push to the theory was also imparted by Keller (1956, 1957), who conducted studies on the diffraction of electromagnetic waves by a convex cylinder and an aperture. Later Keller (1962), initiated the development of the now well recognized geometrical theory of diffraction (GTD). In one of his main contributions he studied the canonical problem of a generalized wedge composed of an edge enclosed by a convex or a concave surface. The most salient feature of the GTD is the introduction of diffraction coefficients, which upon application on the incident rays hitting the geometric singularity would deliver diffracted rays (e.g., just like reflection coefficients are used upon the incident field in a half-space). The diffraction coefficients in Keller’s work failed on the transition regions adjacent to shadow and reflection boundaries. This was later improved by Kouyoumjian and Pathak (1974) who completed Keller’s GTD producing a workable expression to predict the diffraction by a generalized wedge.

Within the specific context of earthquake engineering, the particular contribution from the diffracted field to the total solution has not received much attention. In that area the problem has been traditionally solved in terms of the concept of a scattered field. A possible explanation may be due to practical reasons for both engineers and mathematicians and also, because of the lack of a sound theory of diffraction of mechanical waves. A particular reference to the role played by the diffracted field in the site effects problem can be found in the work of Sanchez-Sesma and its co-workers. For example in Sanchez-Sesma (1985) the author revisited Macdonald’s solution directly in the frequency domain. In that work it was pointed out the fact that large differential motions were introduced by the diffracted field in the presence of shadow zones.

A partition of the field into incident, reflected and diffracted waves is also explicitly displayed by Sánchez-Sesma and Iturrarán-Viveros (2001) in the study of diffraction of plane SH waves by a finite crack in a full space. They obtained a near field solution by superposition of two semi-infinite cracks of the Sommerfeld type. These authors explicitly isolated the diffracted field. Later Iturrarán-Viveros et al. (2010) followed the same technique to find the solution for the diffraction by a cylindrical wave.

In contrast to the antiplane problem, in the case of in-plane waves the analytical solutions related to the diffracted field are just a few and the problem has been studied mainly from a numerical point of view. In this case complexities arise because of the coupling of boundary conditions and mode conversions at interfaces. One of the few particular solutions directly addressing diffraction is the one due to Achenbach (1973), who treated the problem of a semi-infinite slit under longitudinal waves obtained via integral transforms together with the Wiener-Hopf technique and the Cagniard-de-Hoop method. In the slit problem, the total field was shown to be composed of the incident P wave, reflected and diffracted P and SV waves and head waves connecting the diffracted P and SV fronts.

In a recent contribution Jaramillo et al. (2013), using the solution from Kouyoumjian and Pathak (1974) as a building block, proposed a superposition based diffraction technique (SBD), to study site effects due to topographic formations of arbitrary shape. In that method the surface topography is partitioned into several generalized wedges, each one of which contributes with a primary source of diffraction to the total field. Although the resulting method is cumbersome to apply if one wishes to find the complete solution everywhere, it is a useful tool that can be used in the interpretation of results from complex scenarios. It is in that study where we find the motivation for this work.

In our treatment of topographic effects we obtain the diffraction field from the complete solution computed numerically. Our interest lies in establishing a connection between the topographic effect and the diffracted motions. For that purpose we study the two simple problems of cylindrical canyons of semi-circular and rectangular cross sectional shapes. Those two problems differ in the number of diffraction sources. Although from the very concept of diffraction it is clear that this part of the response is related to a geometric effect, the idea has seldom been explored in the study of the effects of topography on earthquake induced ground motions. The physically based partition provided by the diffracted field, is expected to reveal conceptual aspects of the response not evident in the traditional scattering approach. In this paper we study the effects of topography on the incident ground motions using the diffracted field as a physical characterization of the site effects, directly linking the geometry of the scatterer to the motion at a given receiver. The aim of the work is to contribute with the understanding of the topographic or site effects problem by proposing an alternative methodology to isolate the geometric effects and to suggest a partition of the field that can be used to isolate also the mechanical effects.

2 Prediction of the diffracted field

In the classical approach of the scattering problem, the solution is written as the contribution from the free field and the scattered motions, Pao and Varatharajulu (1976). However, while the free-field is well understood, the scattered motions are not. In problems involving arbitrary scattering surfaces, like canyons and ridges, an alternative physically sound way of studying the topographic effect is through the use of the concept of diffraction understood in the sense of optics. As described before, the diffracted motions correspond to that part of the response that cannot be predicted by geometrical methods. For example, if in a given problem the involved scattering surfaces are at least C2-continuous and fully illuminated, a complete, continuous solution, can be constructed based on geometrical methods without formally solving the elastodynamics wave equation. However, if there are shadow zones or geometric singularities, the diffracted field is required in order to smooth out any discontinuities existing in the geometric field. Our proposed analysis method is based on this idea, where the total response is separated into the geometric and diffracted parts. In this way, a direct physical connection between the solution and the topographic features of the scatterer can be established in a source-receiver basis. In this section we describe these alternative partitions of the field and establish different relations between the involved terms. In particular we show how to obtain the diffracted field from the total displacement.

The problem domain is defined in Fig. 1 showing a homogeneous isotropic elastic-half space V0 with a scatterer V1 which are assumed to be perfectly welded along the internal contact surface S. Similarly, SF = traction-free surface of the half-space, S = external remote surface of the half-space where radiation conditions are prescribed. Whenever S is rendered finite in a computational representation, this surface becomes the internal truncation surface SD and SF becomes a finite surface denoted thereby like S ^ F . The scatterer and half-space are mechanically defined by the following parameters: ρ i  = mass density, μ i  = shear modulus, λ i  = Lame constant and where i = 1, and 0 for the scatterer and half-space respectively. The whole system is then subjected to incident harmonic plane waves with a time dependence e i ^ ω t (which is omitted here an hereafter) and whereω = circular frequency and i ^ = - 1 . The purpose of the analysis is to determine the motions in the scatterer and inside the half-space due to the incident wave.

Fig. 1
figure 1

Definition of the problem domain

Fig. 2
figure 2

Classical partition of the total solution into free field plus scattered motions

In the classical formulation we write the total solution inside the half-space like;

u T = u IN + u A R + u S u A 0 + u S
(1)

where uIN = incident field, u R A  = field reflected over the free surface of the half-space in the absence of the scatterer: the field defined from the superposition uIN + u R A u 0A is the free field evaluated along the (fictitious) contact surface S. Due to its mathematical nature we also refer to this component of the total motion as the artificial incoming field. In Eq. (1), it is evident that the scattered contribution uS, is the difference between the total motions and the artificial incoming field u 0A . This construction process of the final complete solution is schematically represented in Fig. 2.

Alternatively, the total solution could be written generalizing the concept of scattering as a relative displacement as we elaborate next: With the aid of Fig. 3, in a first analysis step we solve the problem with the scatterer being removed but leaving the contact surface S as a free surface where traction vanishes. The solution to this intermediate problem can be constructed superimposing the diffracted field uD to the geometric solution u 0 P  = uIN + u R P , where u R P  = physical reflection of the incident rays over the common interface S. The resulting partition is hence written like;

u T = u IN + u P R + u D u P 0 + u D
(2)

where in analogy with the free-field motion u 0 A defined in Eq. (1), we refer to the superposition uIN + u R P u 0 P as the physical incoming field.

Fig. 3
figure 3

Partition of the solution into incident, physically reflected and diffracted fields

If we now add the scatterer domain V1 to the exterior half-space domain with free surface SFS, it will introduce an additional field uM so the total final solution can be written like;

u T = u IN + u P R + u D + u M u P 0 + u D + u M
(3)

In contrast to Eq. (1), the superposition expressed by Eq. (3) is physical, as it naturally separates the geometric contribution—represented by the term u R P  + uD—from that associated to the additional scattered motions uM.

Assuming now that the total field uT has been obtained (for instance, by a numerical method), a comparison between Eqs. (1) and (3) yields the scattered displacements

u M = u T - u P 0 - u D .
(4)

Although we use the superscript ()M to indicate in some sense a mechanical effect in the last term of Eq. (4), it is clear that this part of the field is also affected by the geometry of the coupling surface S. More precisely, the term uM is composed of the refracted field associated to the incident rays and to the diffracted field itself, plus their further diffractions as they interact with the geometric singularities inside the scatterer.

Now, if there is no scatterer (e.g., no impedance contrast), then uM = 0 and there are only geometric modifications captured by the diffraction term uD as;

u D = u T - u P 0 .
(5)

Notice that the general solution expressed by Eq. (3) captures the unmodified half-space solution in which case we have that uM = u 0A  − u 0P  − uD which results after writing

u T = u P 0 + u D + u M u A 0 .
(6)

In order to obtain the diffraction term as indicated by Eq. (5), the problem must first be solved numerically, while the term u P 0 can be obtained by geometrical means.

3 Boundary integral formulations of the scattering problem

In this work we solve the scattering problem in terms of two direct boundary integral statements, derived from the elastodynamics representation theorem (Eringen and Şuhubi 1975; Aki and Richards 2002). The resulting integral equations are used to formulate direct boundary element methods (BEM), through a classical collocation scheme, Banerjee (1994). In each case, the discrete equations are written in the form of force-displacements relationships in terms of a stiffness matrix as in the finite element method: this allows the half-space and the corresponding plane wave excitation to be coupled into existing finite element codes in the form of a half-space-super-element (HSSE). Details of the half-space super element and its implementation in general finite element codes can be found in Zhang et al. (1998). Since in the current work we are interested in the scattered part of the response, only the integral representation for the half-space is considered. The domain V1, corresponding to the scatterer, may also be described by an integral equation or even by a variational statement like in the finite element technique.

In the integral representations we consider formulations with a half-space and a full-space Green tensor (i.e. Lamb and Stokes tensors respectively). The representation in terms of the half-space Green function—simultaneously satisfying radiation and free-surface boundary conditions—results in a BEM discretization involving only boundary elements along the coupling surface S (Fig. 1): that approach is regarded herein as the exact solution and is given in Eq. (7) as follows;

u i S ( ξ , i ^ ω ) = S G i j H S ( x , i ^ ω ; ξ ) t j S ( x , i ^ ω ; n ^ ) d S ( x ) - S H i j H S ( x , i ^ ω , n ^ ; ξ ) u j S ( x , i ^ ω ) d S ( x ) for ξ V 0 .
(7)

The above representation is equivalent to the one derived by (Pao and Varatharajulu 1976) with the only difference being the fact that here we use scattered instead of total fields. In Eq. (7) G i j H S ( x , i ^ ω ; ξ ) and H i j H S ( x , i ^ ω , n ^ ; ξ ) are the displacement and tractions Green’s tensors for a half-space. On the other hand, the terms u S j and t S j refer to scattered displacements and tractions respectively. The resulting BEM algorithm now only involves the mesh along the coupling surface S as shown in Fig. 1. This approach may result computationally expensive since the free surface boundary condition is difficult to satisfy. In this work that condition is enforced using the discrete wavenumber boundary element method (DWBEM), (Kawase 1988; Kim and Papageorgiou 1993).

In contrast, if one uses the full-space Green’s function, the discretization of the half-space must be extended laterally beyond the coupling surface resulting in the following integral representation;

u i S ( ξ , i ^ ω ) S G i j ( x , i ^ ω ; ξ ) t j S ( x , i ^ ω ; n ^ ) - H i j ( x , i ^ ω , n ^ ; ξ ) u j S ( x , i ^ ω ) d S ( x ) - S ^ F H i j ( x , i ^ ω , n ^ ; ξ ) u j S ( x , i ^ ω ) d S ( x ) for ξ V 0 .
(8)

and where G i j ( x , i ^ ω ; ξ ) and H i j ( x , i ^ ω , n ^ ; ξ ) are the full-space displacement and tractions Green’s tensors respectively. The relevant surfaces and parts of the domain in this integral equation were already described in Fig. 1. Since the traction-free surface S F has to be rendered finite in the computational model the resulting representation is an approximation where S F has been denoted like S ^ F .

4 Diffraction effects in the rectangular and semi-circular canyons

In order to study the connection between the diffracted field and the topographic effect, we selected two simple problems. The first case corresponds to the semi-circular cylindrical canyon previously studied by Wong (1982) and Kawase (1988). That problem has two singular sources of diffraction located at the top rims of the canyon. As a second study case, we selected a canyon with a rectangular cross-section and four concentrated sources of diffraction, located at the top and bottoms rims. In the case of oblique incidence both problems will also exhibit diffraction in the corresponding shadow zones.

The canyons are geometrically defined by their depth H = 1.0 km , width B = 2.0 km (for the rectangular case) and radius a = 1.0 km (for the semi-circular case). Note that B = 2a. The half-space is mechanically described by a mass density ρ = 1000 kg/m 3 , a velocity of shear wave propagation β = 1.0 km/s and a Poisson’s ratio ν = 1/3. The excitation is a Ricker pulse of central frequency fc = 1.0 Hz, maximum frequency fmax = 4.0 Hz and a time duration Tmax = 8.0 s. We considered incidence angles (defined with respect to the vertical) corresponding toθ = [0°, 30°]. The analysis was conducted in the frequency domain with a frequency step of Δ f = f max / 32 . The model responses were studied for a normalized frequencyη = ωa/π β.

The amplitude of the frequency domain transfer function associated to the diffracted field, was obtained from Eq. (5), after computing the total field uT numerically, while the physically-based incoming field u 0 P , was obtained analytically. All the analysis were conducted with the formulations given by Eqs. (7) and (8). To integrate the half-space Green’s functions we used the discrete wavenumber boundary element method reported in Kawase (1988) and Kim and Papageorgiou (1993). That technique however, was found to be computationally expensive due to the computations involved in enforcing the traction’s free surface boundary condition by the half-space Green’s function. By contrast, the results obtained with the full-space Green’s function, Eq. (8), exhibited a comparable accuracy at low computational cost. The accuracy of the implementations was tested by benchmarking our results for the semi-circular canyon with those reported in Kawase (1988) and Wong (1982) showing excellent agreement.

Figure 4 displays the results corresponding to the case of an SV wave, incident atθ = 0° andθ = 30° over the rectangular canyon (columns 1 and 2) and semi-circular canyon (columns 3 and 4). Rows 1 to 3 correspond to the vertical displacements, while rows 4 trough 6 show the horizontal components. We show in each case, the spatial distribution of the transfer function at the normalized frequencyη = 1.0 associated to the diffracted field and to the total field. We also depict in the bottom part of the figure the related synthetic seismograms along these same surfaces (rows 3 and 6). For the rectangular canyon the range of x [ - 1.0 , 1.0 ] corresponds to the bottom of the canyon, while x > 1 corresponds to the canyon walls.

Fig. 4
figure 4

Response of the rectangular and semi-circular canyons to SV waves incident atθ = 0° andθ = 30° obtained with the exact BEM algorithm. The results shown inrows 1 and2 correspond to the spatial distribution over the free surface of the frequency domain transfer function at the characteristic frequency fc = 1.0 Hz associated with the diffracted (row 1) and total (row 2) vertical displacement component. The results inrow 3 are the synthetic seismograms over the canyon surface in each case. Thefirst two columns correspond to the rectangular canyon andcolumns 3 and4 to the circular canyon

For the rectangular shape under vertical incidence, the diffraction effect is revealed by the large amplification near the rims of the canyon, where the amplitude of the transfer function corresponding to the vertical displacement component at η = 1.0 reaches a value close to 3.0. In this case, as the main incident front encounters the bottom corners of the canyon, these geometric singularities become sources of diffracted waves. From analytical solutions of the diffraction of a P wave by a slit (Achenbach 1973), we know that such a singularity, generates head waves and P and SV diffracted waves. These last two exhibiting cylindrical fronts. In the current problem all these fronts propagate vertically over the canyon walls and horizontally over the bottom surface.

To reinforce the analysis we also show in Fig. 5 snapshots of the propagation patterns at four different time instants. At t = t1, the incident SV wave has already been diffracted by the corner singularities at the bottom of the canyon, producing cylindrical diffracted P and SV waves. Over the horizontal bottom surface, the leading diffracted P wave travels in phase with a head wave (required to match the traction free boundary condition). The diffracted SV wave follows along the same bottom surface and both (the P and the SV component) are tied by the head wave extending from the surface to the inside of the computational domain. Similarly, over the vertical walls and propagating upwards, one encounters the incident grazing SV wave, the cylindrical diffracted SV wave, the cylindrical diffracted P wave and a joining head wave. It is interesting to observe how the diffracted SV wave traveling alone over the bottom surface is enough to match the traction free boundary condition, while it is coupled to the incident front over the vertical wall.

Fig. 5
figure 5

Response of the rectangular and semicircular canyon to a vertically and a 30° incident SV wave. The snapshots correspond to full particle motions and time progresses fromleft toright

The incident and diffracted cylindrical waves, generated after the first interaction with the bottom singularities, will also experience further diffraction as they encounter the opposite corners of the canyon. To clarify this aspect of the response, one could use Huygen’s principle ideas to think of the bottom corners like sources of diffracted cylindrical and head waves. If we refer to this diffracted field as the primary diffraction, then it is clear that the primary diffracted field will experience further diffraction as it encounters the opposite corners located over the horizontal and vertical surfaces. This field should, by obvious reasons be termed the second order diffraction effect.

Eventually, the process of activation of Huygen’s sources of diffracted waves continues indefinitely forcing energy to become trapped along the canyon surface in the form of higher order diffraction. This effect is observed in the synthetic seismograms and also from the snapshot at t = t2, as the increase in duration of the signals over the bottom surface. The larger amplification at the canyon rims observed in the frequency domain response (see Fig. 4) is now evident since at these locations the incident, reflected and primary diffraction fields will converge.

An additional point of interest regarding the diffraction field, is observed in the response inside the computational domain in the snapshots at t = t3 and t = t4. In the first case, the cylindrical SV diffracted wave, propagates in phase with the first reflected field. It is observed how these two waves deplete each other as they move away from the canyon. In the second case, the diffracted field generated by the top corners, propagates together with the reflected wave, in such a way that it restores the front in the far field and as it moves away from the canyon. From these two opposite cases it is clear how the diffracted part of the response, has the ability to complete or deplete the field as required.

In the cases of incidence at 30° (columns 2 and 4 in Fig. 4), the propagation pattern is complicated since the incoming field is composed of three sets of rays leading to the activation of multiple diffraction sources. The first diffractor is activated when the incident SV wave, meets the corner singularity in the illuminated zone as depicted in the first snapshot from Fig. 5. Once again the cylindrical diffracted P and SV waves are clearly identified. This primary diffraction must deplete, in the far field, the reflection of the P wave in the left-most canyon vertical wall. Similarly, in the snapshot at t = t2 the main incident SV wave and its associated diffracted field, experiences second order diffraction by the right bottom corner, while the reflected P wave is undergoing its first order diffraction. In this case the cylindrical P and SV diffracted waves propagating over the horizontal surface, will diffract once again by the right bottom corner. In the snapshot at t = t3, the diffraction of the main incident wave is already filling up the shadow zone and experiencing third order diffraction by the canyon top singularity. At that instant the reflected P wave is experiencing second order diffraction. Finally at t = t4, the multiple diffracted waves are starting to complete the fronts of the free field, both in the shadow and illuminated zones.

A similar behaviour is observed in the case of the semi-circular canyon under a vertically incident field. For this problem the time domain response is depicted in the snapshots in Fig. 5. In the frame at t = t1 the incident field has just been reflected in the form of SV and P waves. Along the canyon surface these three waves travel in phase, until they experience the first diffraction event as shown at t = t2 where now the field is composed of the incoming wave plus the diffracted P , S V and head waves. At t = t3 the diffracted field is approaching the edges where it will undergo second order diffraction. In that frame and at t = t4 the reflected SV front is already being restored. Due to space limitations the complete set of results is included as supplementary material to be accessed by the reader.

5 Conclusions

We have studied the incidence of the diffracted field in the problem of seismic site effects in earthquake engineering. For that purpose we have established the difference between scattered and diffracted fields, and have also established the connection between these two parts of the motion. This has led to the idea of a physical incoming motion, which is different to the classical definition of free field motion, which we have referred to like the artificial incoming field. With these definitions at hand, it has been shown that the solution to a scattering problem can alternatively be expressed like the superposition of the physical incoming field plus purely geometric and mechanical effects.

Our study has been conducted with regards to the two simple problems of a semicircular and a rectangular canyon with two and four singular sources of diffraction respectively. We have subjected both canyons to vertically and obliquely incident in-plane P and SV waves and have found the total response in terms of frequency domain transfer functions along the canyon surfaces and in terms of synthetic seismogrmas for receivers over the canyons surfaces. The transfer function associated to the diffracted field in each case, has been obtained after subtracting from the total response, the physically based incoming motion. The responses were evaluated with two direct boundary element methods.

A comparison between the total response and that containing only the diffraction term, has shown that the true mark of the topographic effect is effectively captured by the diffracted field. Furthermore, we have shown that when the problem is studied from the point of view of the diffracted field its is possible to connect the total response with the topographical features of the scatterer even in the time domain. Overall, larger amplifications have been found in the rectangular canyon.

The results from this study indicate that a thorough method based on diffraction and intended to support the understanding of the site effects problem, requires the solution to simpler canonical problems, where a detailed knowledge of the diffracted part of the motion becomes available. In particular, it was observed how the diffraction appeared in the form of cylindrical waves and head waves, exhibiting a complex dependence of its amplitude and phase on parameters like the angle of incidence, the geometry of the singularity and the distance from the diffractor. For instance, in the case of the rectangular canyon under vertically incident SV waves, it was observed how the cylindrical front associated to the diffracted SV wave was able to match the traction free boundary conditions traveling independently over the bottom surface and in phase with the incident field over the canyon vertical walls. The function describing this cylindrical wave clearly has a complex angular dependence being able to capture these two opposite conditions. Knowing such functions will allow the analyst to anticipate the relevance of the site effects problem based on simple parameters of the scatterer.