Abstract
A method to study the response of surface topographies submitted to incident waves is presented. The method is based on the superposition of diffracted sources described in Jaramillo et al. (2013). Since the technique proceeds in the frequency domain in terms of the superposition of incident, reflected, and diffracted waves, it has been termed like a superposition based diffraction approach. The final solution resulting from the superposition approach takes the form of a series of infinite terms, where each term corresponds to diffractions of increasing order and of decreasing amplitude generated by the interactions between the geometric singularities of the scatterer. A detailed, step-by-step algorithm to apply the method is presented with regard to the simple problem of scattering by a V-shaped canyon. In order to show the accuracy of the method we compare our time and frequency domain results with those obtained from a direct Green’s function approach. We show that fast solutions with an error of the order of 6.0% are obtained.
1. Introduction
The presence of surface or subsurface topography has been recognized for many years as an incident factor in earthquake induced ground motions resultant at a site. Classical and well-documented examples can be identified in the large ground accelerations registered at the Pacoima Dam during the 1981 San Fernando, California, earthquake and in the records of Tarzana hill during the 1994 Northridge earthquake [1, 2]. Mathematically, the quantification of topographic effects involves the solution of the elastodynamic wave scattering problem, which in the case of geometries of arbitrary shape gives rise to highly complex responses requiring a numerical solution. Despite the fact that currently existing numerical techniques together with the continuous increase in computational power provide the analyst with an actual capability to conduct large-scale simulations for highly realistic scenarios, closed-form solutions still remain important as validation frameworks, mainly since they derive into important conceptual understanding of the problem of site effects. In this work we introduce a method to determine analytic approximations to the scattering of waves induced by topographic surface irregularities, which in contrast to other analytical treatments available in the literature has the advantage that the solution is progressively built based upon physical arguments giving it great flexibility as a validation tool for numerical implementations or allowing its use as first-hand interpretation of results associated with very complex scenarios.
In the case of scattering of in-plane waves, the number of analytic solutions is limited. For instance, the scattering of incident and waves by a semicircular canyon was only recently found by Lee and Liu [3] using a wave function expansion approach, while in the case of horizontally polarized shear waves there is only a handful of contributions: an important work was conducted by Trifunac [4], who found the frequency domain solution to the scattering of plane waves by a semicircular canyon and a semicircular valley using a separation of variables approach. Similar solutions for the scattering of waves by irregularities of various shapes constructed in terms of wave function expansions have also been developed by a region matching technique in Tsaur and Chang [5]; Tsaur et al. [6]; Tsaur [7]; Liu et al. [8]; Tsaur [9]; Han et al. [10]; Gao et al. [11]; Zhang et al. [12]; Gao and Zhang [13]; Chang et al. [14]; and Tsaur and Hsu [15]. In all of these works the solutions take the form of infinite series or expansions with an infinite number of terms and despite the fact that they correspond to the simplest cases of scattering of scalar waves by strongly idealized geometries, these solutions shed light on conceptual understanding of the problem of site-specific response. One particular limitation of solutions in terms of infinite series is the fact that practical applications are only possible after a thorough numerical analysis, which may be a difficult task due to the nonuniform convergency properties of the expansion functions. These limitations may become important factors when the size of the spatial domain or the frequency content of the analysis is changed.
In a recent contribution, Jaramillo et al. [16] laid down the basis of a superposition analysis technique, where, in contrast to the standard approach behind series solutions, these authors used a partition of the field based upon a physically derived incoming motion representing the geometrical plus the diffracted field. Since the geometrical field can be obtained exactly, using well-known reflection laws, the only approximation to the solution is related to the diffracted field. In the resulting superposition based diffraction (SBD) approach these diffracted terms are derived after representing the topographic irregularity as a superposition of overlapped wedges: each one contributing with a source of diffraction emanating cylindrical waves. The resulting solution is also an infinite series, but by contrast with wave function expansion techniques, here each term corresponds to a diffracted wave of decreasing amplitude whereby truncation is decided based upon accuracy. Depending on the ratio between the wavelength of the incident motion and the characteristic dimensions of the scatterer, the series can be truncated after considering just a few terms which results in a highly economical approach.
In this paper we present a step-by-step description of the proposed SBD algorithm. In the first part we explain our idea of a physically based incoming motion, corresponding to the geometrical (or optical) field, and establish its connection with the classical free-field motion commonly used in earthquake engineering applications. We then review the fundamental solution for the scattering of waves by a generalized wedge as presented in Jaramillo et al. [16]. The paper subsequently describes how to conduct the partition of the computational domain into subdomains according to the regions of existence (or absence) of incident and reflected rays and depending upon the number of diffraction sources introduced by the topographic irregularity. The method is then validated against the solution for a symmetrical V-shaped canyon reported by Tsaur and Chang [5]. We present approximations considering diffracted waves and their interactions with adjacent wedges up to third order. This yields accuracy comparable to the analytic solution. Moreover, in order to show the applicability of the method, we determine the response of a 25° V-shaped canyon using different number of diffraction terms. The SBD results were compared with those obtained by a boundary element method (BEM) computer package. The comparison was conducted in the frequency and in the time domain where we measure the error associated with the different sources of diffraction. Economic solutions, with errors of the order of with respect to the numerical results, are obtained after including just a few terms in the series.
2. Alternative Representation of the Total Response
The basis of the current superposition based diffraction (SBD) approach is the linear character of the problem. This allows the total solution to be written in terms of the addition of different and arbitrary superpositions. One such partition is based upon the usual earthquake engineering definition of free-field motion, where the total response is constructed by the addition of an incident field ; its reflections in a half-space after removing the scatterer ; and a relative or scattered field . Recalling the definition of free-field motion, commonly used in earthquake engineering and given by , allows us to write the total field like
By reasons that will become apparent later, we refer to the free-field term like the artificial incoming motion.
An alternative partition of the field, introduced in Gomez et al. [17] and constituting the basis of our method, is now explained with reference to Figure 1 which schematically describes a scattering problem where the domain (a) has been partitioned into the scatterer (c) and its supporting half-space (b), respectively.
Based upon the above partition of the total problem into subdomains we can subsequently write for the total field:where the first three terms correspond, respectively, to the incident field; its reflections over the free surface (notice that becomes exposed after having removed the scatterer from the half-space); and the diffracted field. This last term in the solution is required in order to restore continuity of the field since the contribution from is incomplete due to the geometric singularities existing along . If we identify the first two terms in (2) like the optical field we can rewrite the total solution like
In (3) is an additional field introduced by the scatterer . The name artificial incoming motion, previously coined to the engineering free-field definition , is now evident since such definition, according to (1), leads to the concept of the scattered field which has mainly a mathematical meaning. If we now consider the term like an alternative scattered motion or relative displacement between the total solution and the optical field, then (3) can be written likewhich is analogous to the classical partition given by (1) in terms of the free-field and the mathematical scattered motion. The alternative and classical scattered fields are easily shown to be connected by
Moreover, in problems involving only a topographic irregularity where and the total field can be written like
Since the optical field (OF) or physically based incoming motion can be obtained analytically, the construction of the total solution becomes feasible as long as we find a way to compute the contribution from the diffracted field. This term can be obtained following the work from Jaramillo et al. [16] in terms of the diffracted field for a generalized infinite wedge, and after representing the topographic irregularity as a superposition of wedges of different inclinations and perceiving the incident wave at different angles. Although the formulation of the problem in the standard form of (1) is more suitable for numerical treatments of the problem, the superposition given by (6) has a stronger physical basis and turns out to be highly convenient for an analytical approach.
2.1. Fundamental Solution
In the SBD method the geometry of the topographic irregularity is constructed via superposition of rectilinear segments using as a fundamental entity the simple wedge shown in Figure 2.
The corresponding diffraction field generated by the interaction of a plane or cylindrical wave with the corner singularity in the wedge is given in
This solution was proposed by Kouyoumjian and Pathak [18] in the context of propagation of electromagnetic waves after using the well-established geometrical theory of diffraction (GTD), originally proposed by Keller [19], and particularized to the case of horizontally polarized shear waves incident against surface topographies by Jaramillo et al. [16]. In (7) is radial coordinate of the field point measured from the vertex of the wedge, is angular coordinate measured with respect to the reflection boundary, is incidence angle measured with respect to the reflection boundary, is wedge angle (with being a factor between and ), is radius of the incident cylindrical wave (for the diffraction of a cylindrical front), is wave number, and is velocity of wave propagation. The remaining terms appearing in (7) are defined as follows:
3. The Superposition Based Diffraction Method
The algorithm to construct the solution is summarized next with reference to the schematic simple V-shaped canyon shown in Figure 3. First, the surface geometry has to be approximated by overlapping wedges, each one contributing with 2 surfaces of reflection and 1 source of diffraction. In a second step, zones illuminated by the incident and reflected waves are identified and the domain is divided according to the number of rays existing in these different zones. The third step involves the identification of boundaries between illuminated and shadow zones and corresponding to zones of discontinuity to be filled out by the diffraction field and generating a second subdomain partition. The fourth step requires the computation of the optical and diffracted fields in each one of their subdomains of existence: in this last step the interaction between sources of diffraction must also be considered. The inclusion of each new source amounts to an application of the expression given in (7) with the values of the parameters properly adjusted to the particular case and with each one of them contributing with an additional term to the series solution. The details involved in the partition of the total domain into subdomains are elaborated next.
3.1. Representation of the Surface Geometry
In the first analysis step the surface topography is represented by a superposition of several wedges of the type described in Figure 2. A wedge is defined by the angular parameter and by the angle of incidence of the plane wave as given in (7). Here we denote the -th wedge by and each newly introduced -wedge contributes at the same time with a diffractor (or source of diffraction) denoted by and with two free boundaries. The free surfaces are reflectors and therefore regarded as continuous sources. The wedges required to represent the V-shaped canyon are shown in the left part of Figure 4. On the other hand, the arrows show the source of diffraction and the two reflection surfaces contributed by each wedge.
3.2. Partition of the Domain
The computational domain is now partitioned into convenient subdomains according to the regions of existence of the fields contributed by the different sources. In a first step we consider the main incident wave and its reflections at the free boundaries. In what follows we will denote this partition like the free-field-partition which is clearly dependent on the angle of incidence. Here we consider the regions of existence of different plane waves in terms of rays. In Figure 5 we show the free-field-partition for an asymmetric V-shaped canyon subjected to a vertically incident wave. The first region is determined by the incident wave (denoted in the figure by the ray ). In this problem the incident ray exists throughout the complete domain. The second region is the reflection zone defined by the existence of the downward vertically propagating waves (denoted in the figure by the rays ). These appear after reflection at the horizontal free boundaries of the half-space. The last and final subdomain is the region enclosed by the diagonal lines defining the region of existence of rays reflected at the inclined free faces of the canyon (denoted in Figure 5 by the rays and ).
In the final partition step regarding the continuous sources, the above three partitions are added (or superimposed) showing the contribution from the incoming field in terms of rays at all the different points inside the problem domain. After considering interceptions of regions a total of 6 different zones are identified: these are enclosed by circles in Figure 5. It should be noticed that only the rays originated at the inclined boundaries must be phase-corrected considering the used system of reference. Also, notice that the rays and are reflected downward as long as the angles . In a more general case, the rays reflected by the inclined free surfaces of the canyon may experience subsequent reflections at the horizontal part of the half-space leading to a different free-field-partition.
The second superposition step consists in the consideration of the diffracted field. For this purpose we perform a second partition of the computational domain, analogous to the one used for the free-field, but now dictated by regions of existence of diffraction terms. As already pointed out each wedge contributes with cylindrical diffracted waves emanated from a source at its apex producing a subdomain related to the region of existence of the originating wedge. This partition is described in Figure 6 and in what follows we refer to such specific consideration of the domain as the diffracted-field-partition.
The corresponding diffracted rays have been labeled , , and . For instance, the diffracted rays exist throughout the complete domain, while the rays labeled and exist only inside the domains delimited by the lines and , respectively. The final superposition of the 3 subdomains is represented by crossed circles with their related diffracted rays (see Figure 6).
3.3. Higher-Order Diffraction Contribution
In the above superposition process, the geometric singularity associated with each subdomain may diffract either the main plane front corresponding to the incoming field or the cylindrical waves originated at adjacent wedges. Hereafter we refer to the first effect like primary diffraction while subsequent diffraction events are referred to like higher-order diffraction. Both the first- and higher-order diffraction contributions represent terms in the series solution. However, since each diffracted wave has an amplitude smaller than its originating source, only a few terms need to be retained in the series. To clarify, let us consider the diffraction field generated by the central wedge with apex at point in Figure 6. The diffracted field generated by this source can be written likewhere the first subscript indicates the originating source, the last one indicates the diffractor, and the number of subscripts indicates the order of the diffraction. Accordingly, the term describes the diffraction of the main front at point ; the term is the second-order diffraction (indicated by two subscripts) introduced by the diffractor (indicated by the last subscript) of the cylindrical wave with source at (indicated by the first subscript). Similarly, the term is the third-order diffraction produced by the diffractor which diffracts a cylindrical wave with source at , diffracted at and diffracted again at . Figure 7 shows a schematic representation of snapshots of the different diffracted fronts generated by the wedges of the V-shaped canyon. At the time instant the cylindrical front identified at the lower apex corresponds to the first-order diffraction. Similarly, at time , there are two sets of cylindrical fronts generated at the upper wedges. The large front is once again the main or first-order diffraction, while the small ones correspond to second-order diffraction experienced by the cylindrical wave generated at . These fronts subsequently experience third-order diffraction at .
The solution process now reduces to the addition of the free and diffracted fields as imposed by their domains of existence (previously defined like free-field-partition and diffracted-field-partition). The diffracted field is obtained after successive applications of the fundamental solution given in (7) with arguments corresponding to a plane or to a cylindrical wave according to each case. In order to clarify this last step, we show in Table 1 order of the arguments for the evaluation of the second-order diffraction resulting after the cylindrical wave generated at the source point interacts with the wedges with vertices at points and . According to the notation introduced previously these terms are denoted like and , respectively.
3.4. Total Solution
The complete solution now involves the addition of the field existing within each subdomain. In general we can writeDespite what has been implied by (10), the solution takes the form of an infinite series where each term with an increasing number of subscripts corresponding to higher-order diffraction has a decreasing amplitude. Our solution has the advantage over those found using eigenfunction expansions that the approach does not require the problem geometry to correspond to a separable system of coordinates and the fact that here each term of the series has a physical meaning associated with the diffraction contributions of various orders. As such, the approximation is conducted on physical grounds rather than on convergency analysis based on mathematical arguments. Moreover, since the truncated terms are related to the diffracted field, which is known to be frequency dependent, the stopping criteria depend on the dimensionless frequency of the particular problem under study.
4. Results
4.1. Evaluation of the Symmetrical V-Shaped Canyon Studied by Tsaur and Chang [5]
As a validation of the proposed SBD approach we solved the V-shaped canyon previously studied by a wave function expansion method in Tsaur and Chang [5]. The canyon response was computed for values of the depth-to-width aspect ratio corresponding to covering both shallow and deep canyons. The analyses were conducted for a dimensionless frequency and for unitary values of the wave propagation velocity and mass density.
Figure 8 compares the spatial distribution of the amplitude function over the free surface. The first row compares the solution by Tsaur and Chang [5] and the SBD results for 1, 2, and 3 orders of diffraction, respectively, while in rows 2, 3, and 4 we show these same results independently. Similarly each column in Figure 8 compares the wave function solution with the SBD results for the different values of the aspect ratio. It is observed how, independent of the number of orders of diffraction considered in the SBD solution, the results degrade in the direction of increasing aspect ratio. However, it is evident that 3 orders of diffraction suffice to reach a solution comparable to the one reported by Tsaur and Chang [5] even for the deep canyon (), while 2 orders of diffraction are enough in the case of a shallow valley. If one takes advantage of the symmetry of the domain the solution with 3 orders of diffraction amounts to a total of 8 terms in the SBD series while the same result requires the evaluation of 24 terms if one follows the region matching technique formulated by Tsaur and Chang [5].
4.2. Response of a 25∘ V-Shaped Canyon
In order to test the accuracy of the SBD technique, we conducted frequency and time domain analysis of the simple V-shaped canyon topography using our proposed method and a numerical boundary element (BEM) based algorithm. We first found the response of the system to harmonic plane waves of unit amplitude described in terms of the dimensionless frequency , where is the canyon height and is the corresponding wavelength of the incident wave. In all our analyses we considered a canyon while the depth parameter was varied according to the dimensionless frequency . We obtained the response at values of between and . The frequency domain results were used later in order to find the time domain response using a Fourier transformation approach. For that purpose we used a Ricker pulse defined according to , where with being central frequency and being initial time for the intense phase. In this work we used a pulse with a central frequency and the initial time was adjusted to generate a time window conveniently fitting the selected computational domain.
Figure 9 displays the frequency domain results in terms of the spatial distribution of the amplitudes of the transfer function over the canyon surface. Column 2 shows contour maps of equal amplitude over the free surface of the canyon, for a range of dimensionless frequencies obtained with the SBD approach and the BEM algorithm while column 1 displays results corresponding to Fourier spectral amplitudes at 4 constant values of the receiver coordinate along the canyon surface. Similarly, columns 3 and 4 display the spatial distribution of the amplitude function at selected values of the dimensionless frequency.
(a)
(b) BEM
(c)
(d)
(e)
(f) 1st-order diffraction
(g)
(h)
(i)
(j) 2nd-order diffraction
(k)
(l)
(m)
(n) 3rd-order diffraction
(o)
(p)
Over a specific observation point the response is obtained considering the contribution from the optical field (previously found in closed form) plus the diffracted field, which in the SBD method is approximated by a truncated series. It must be emphasized that when the exact solution is represented in the frequency domain, it must contain an infinite number of diffraction terms. However, depending on the value of the dimensionless frequency parameter in the SBD method, a limited number of terms may be enough in order to achieve a predefined accuracy. The difference between our approach and other series solutions is the fact that in the proposed SBD method each term in the series corresponds to a diffracted wave with an amplitude that decays with , while in alternative solution techniques each term on the series represents a part of the total field.
The results displayed in columns 1, 3, and 4 in Figure 9 reveal how the SBD solution approaches the response obtained with the BEM numerical algorithm as we move away from the scatterer and as the dimensionless parameter increases. This trend can be explained from the dependency of the diffracted field on as identified in (7). On the other hand, the amplitudes reported in column 2 show that the larger deviation of the SBD results with respect to those obtained with the BEM algorithm appear near the apex of the conforming wedges. In particular, the results obtained with SBD technique with diffraction order near the apex exhibit a discontinuity in the total field which rapidly vanishes as the number of orders of diffraction is increased.
When evaluated in the frequency domain the performance of the SBD technique can be summarized as follows. When a large number of diffraction terms (3 to 5 terms) are required in the total solution in order to recover the half-space response expected for small scatterers. This is understood after realizing that the frequency independent incoming motion introduces a discontinuity that, in the low frequency regime, must be smoothed out almost immediately. As we have repeatedly highlighted, the amplitude of the diffracted field depends on ; therefore a low frequency implies a large amplitude in this component. At the same time, this fact implies that the higher-order diffraction terms retain large amplitudes and the scatterer contains a high level of internal interaction. By contrast as the total response requires only a few diffraction terms since now this component of the total field decays very fast.
As an additional validation Figure 10 displays frequency-space contour maps for the ratio between the amplitude of the transfer function obtained with the BEM algorithm and with the SBD method, for receivers over the free surface. Column 1 shows the BEM to SBD ratios, while column 2 displays ratios between the SBD results obtained with different orders of diffraction. It can be observed that the solution is very accurate in the high frequency regime and far away from the scatterer and it exhibits inaccuracies at low frequencies near the conforming wedges. These regions correspond to cases in which the parameter is small and there is an important contribution from the higher-order diffraction terms. The last contour map from column 1, corresponding to the solution with up to orders of diffraction, shows values close to almost everywhere. This is also evident in the last contour map from column 2 where the improvement from the solution considering orders of diffraction is evident.
(a) BEM/SBD-first-OD
(b) SBD-first-OD/SBD-third-OD
(c) BEM/SBD-second-OD
(d) SBD-first-OD/SBD-second-OD
(e) BEM/SBD-third-OD
(f) SBD-second-OD/SBD-third-OD
The computations with the SBD approach, which were progressively obtained with , , and orders of diffraction and error measures defined like , were computed for each solution. Figure 11 displays the error for receivers over the free surface and measured with respect to the BEM response considering up to 3 orders of diffraction. In computing the error we considered a range of frequencies from to . A large difference between the SBD solution with only 1 diffraction order and those with 2 and 3 orders is clearly observed. This is due to the secondary reflections experienced by the first diffracted waves over the free surface of the half-space where its amplitude is doubled. The above error measure is also shown in Figure 12 for the real and imaginary component of the response.
(a) Error for the real component of the transfer function
(b) Error for the imaginary component of the transfer function
Figure 13 displays synthetic seismograms for receivers over the canyon surface computed with the BEM algorithm and with the SBD technique using 1, 2, and 3 orders of diffraction. A qualitative comparison of these results shows how the only difference between the two sets appears in a diffracted front emanating from the inferior wedge. This front is missing from the synthetics corresponding to a single diffraction term. Such discontinuity is not present in the synthetics for 2 and 3 orders of diffraction since it corresponds to a second-order diffraction that originates when the main diffracted front from the upper wedge interacts with the inferior wedges from the canyon. In this case the second-order diffracted field experiences an additional diffraction event producing a third-order front that travels almost in phase with the second-order front. These fronts separate from each other as the distance from the scatterer increases.
(a) BEM
(b) SBD-first-OD
(c) SBD-second-OD
(d) SBD-third-OD
Figure 14 shows synthetic seismograms computed with the BEM algorithm and with the SBD technique using 1, 2, and 3 orders of diffraction at 4 different locations. It is observed, once again, how the time domain response is recovered almost completely by considering only diffraction effects up to third order. The contribution from higher-order terms has an amplitude which is already very small with respect to the amplitude from the incident field. Moreover, these higher-order diffraction effects are highly delayed with respect to the incident field. This is also verified by the error measures shown in Figure 15 corresponding to receivers over the free surface and calculated with respect to the numerical BEM solution. This field is compared with the results obtained from the SBD technique for 1, 2, and 3 orders of diffraction. Clearly, the smaller error occurs for the solution with third-order diffraction effects. It must be realized that in the time domain the magnitude of the incident field is highly important since large amplitudes imply that diffracted fields required large distances to vanish.
(a) History at
(b) History at
(c) History at
(d) History at
5. Conclusions and Further Work
We have presented a superposition based diffraction technique (SBD) to study scattering of waves by surface topographies lying on elastic half-spaces. In the method the surface geometry is approximately represented by a series of superimposed wedges. The solution is then constructed after partitioning the total domain into subdomains according to the contribution from each wedge in terms of incident, reflected, and diffracted rays. The diffracted waves are classified into first-order diffracted waves, corresponding to the diffraction of the main front with each one of the corner singularities contributed by each wedge and higher-order diffracted waves corresponding to the interaction between adjacent wedges. The final solution is then easily built by adding incident, reflected plus first- and higher-order diffracted waves. Although the exact solution involves an infinite number of diffracted waves, we have shown by comparing our results with those from a numerical boundary element algorithm that the method reaches engineering accuracy with just a few terms. For instance, a simple approximation with only one diffraction term leads to an error of the order of 5% while it is reduced to values as low as 3% when we consider up to third-order diffraction events. In order to illustrate the application of the technique we have studied the simple case of a V-shaped canyon, which is widely documented in the engineering literature. The method is applicable not only as a simple validation tool for numerical implementation but also as analysis tool to interpret results related to highly complex scenarios obtained with robust numerical algorithms.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This project was conducted with financial support from “Departamento Administrativo de Ciencia, Tecnología e Innovación, COLCIENCIAS” and from Universidad EAFIT through Research Grant no. 1216-403-20661.