Introduction

Mathematicians and physicists have long studied the physics of waves at structured interfaces: Back in 1898, Lamb wrote a seminal paper on reflection and transmission through metallic gratings1 that then inspired numerous studies on the control of surface electromagnetic waves. The concept that periodic surface variations could create and guide surface waves has emerged in a variety of guises: Rayleigh-Bloch waves for diffraction gratings2,3, Yagi-Uda antenna theory4,5 and even edge waves localised to coastlines6 for analogous water wave systems. Most notably in the last decade, the discovery of spoof plasmon polaritons (SPPs)7 has motivated research not only in plasmonics8,9,10 but also in the neighbouring, younger, field of platonics11, devoted to the control of flexural (Lamb) waves in structured thin elastic plates. The extension of ideas such as cloaking12,13,14 to plasmonics15,16,17,18,19,20 is highly motivational as it suggests that one can take the concepts of wave physics from bulk media, as reviewed in21, and apply them to surfaces. In platonics, implementations of ideas from transformation optics22,23,24 show how, all be it for the limiting cases of thin elastic plates, similar concepts can take place in another arena of physics, in this case mechanics and elasticity. Unfortunately in much of elastic wave theory, and for many applications, one actually has the opposite physical situation to that of a thin plate, that is, one has an elastic material with infinite depth; on the surface of such a half-space elastic surface waves, Rayleigh waves, exist25,26 that exponentially decay with depth and have much in common conceptually with surface plasmons in electromagnetism. It is therefore attractive to investigate whether concepts proven in plasmonics can be translated across despite the underlying governing equations having many fundamental differences. Some experimental work has taken place in broadly related scenarios such as the attenuation of Rayleigh waves at high frequencies in a marble quarry with cylindrical holes27 and in piezo-electric substrates with pillars28, but with differing aims and not at frequencies relevant for seismic waves. Some work on structured elastic media, going beyond elastic plates, to try and create seismic metamaterials29,30,31 is underway with some success either with subwavelength resonators31 or with periodic structuring within the soil30 and trying to still utilise flexural wave modelling. Our aim here is complementary in that we want to implement ideas from transformation optics into this elastic surface wave area and investigate what can be achieved.

The desire, and need, to control the flow of waves is common across wave physics in electromagnetic, acoustic and elastic wave systems32. Acoustics is, mathematically, a simplified version of full elasticity with only compressional waves present: full elasticity has both shear and compression, with different wave speeds, and full coupling between them - leading to a formidable vector system. The quest for a perfect flat lens33,34 with unlimited resolution and an invisibility cloak12,13 that could conceal an object from incident light via transformation optics35 are exemplars for the level of control that can, at least, in theory be achieved for light36 or sound37,38,39. An invisibility cloak for mechanical waves is also envisioned40,41 for bulk waves, and experiments using latest developments in nanotechnology support these theoretical concepts42. However, some form of negative refractive index appears to be necessary to achieve an almost perfect cloak. Unfortunately, a negative index material emerges from a very complex microstructure that is not feasible for the physical length-scales, that for seismic waves at frequency lower than 10 Hz means wavelengths of the order of a hundred meters, in geophysical applications. Notably one can design locally resonant metamaterials that feature deeply subwavelength bandgaps31 and so could be used in seismic and acoustic contexts, but these effects are observed over a limited frequency band; this is being improved to allow for simultaneous protection and cloaking24 and has potential. However, here we look at protection and cloaking, applied to low frequency seismic waves, using the concept of lensing.

Elastic waves, in the same way as light, are subject to Snell’s law of refraction, and an appropriate choice of material properties leads to spectacular effects like those created by gradient index (GRIN) lenses. Compared to a classic lens where ray-paths are bent through discontinuous interfaces causing losses and aberrations, GRIN lenses are obtained with a smooth refractive index transition. Rayleigh and Maxwell themselves studied GRIN lenses, notably Maxwell’s fisheye whose spatially varying refractive index was later shown to be associated with a stereographic projection of a sphere on a plane43. As noted in43 GRIN lenses have been mainly studied and implemented for optical applications, or wave systems governed by the Helmholtz equation. The ideas behind transformation optics are not limited to metamaterial-like applications and have contributed to recent advances in plasmonic light confinement by touching metallic nanostructures44,45.

A classical example of a GRIN lens is the circular Luneburg lens43,46, once a plane wave enters the lens, its radially varying refraction index steers the ray-path towards a focal point located on the opposite side of the lens47. The eponymous fisheye lens of Maxwell48 is another well documented and fascinating GRIN lens and it has been proposed as a non-Euclidean cloak. To date, applications have been mainly limited to scalar wave systems governed by Helmholtz type operators for transversely polarized electromagnetic waves49, for pressure waves or platonics where composite and thickness modulated plates have recently been proposed50,51,52,53. In full elasticity, where the propagation is described by the vector Navier equation that is neither simply a Helmholtz equation or scalar, a proof of concept is still missing. Experimentally, for seismic waves, the realization of any proposed lensing arrangement becomes a real challenge because wavelengths range from 20 m to 500 m. For real applications, wave control must be achieved over a broad frequency band, to cover the various wavelengths of interest; unlike other metamaterial cloak designs such as those based around subwavelength resonators24, the lens proposed here is effective across a very broad spectrum of frequencies.

In civil engineering the structuring or reinforcement of soils is commonplace with geotechnical solutions aimed at improving soil seismic performances54 (e.g. soil improvement by dynamic compaction, deep mixing and jet grouting), typically implemented prior to construction of structures, aimed to rigidify or decouple the building response and not at rerouting the seismic input. In the conceptual lens we design, the structuring of the soil is feasible and we use material parameters typical of poorly compacted sediments. Using the ideas of transformation optics, and detailed 3D numerical simulations, we show that a square arrangement of four Luneburg lenses can completely reroute waves around an area for seismic waves coming with perpendicular incident directions (e.g. x and y in Fig. 1b). Not only are the waves rerouted leaving the inner region protected, but the wave-front is reconstructed coherently after leaving the arrangement of lenses.

Figure 1: The vertical component of the displacement field is depicted at the onset time along with the meshed model of the halfspace that is used as input for the SPECFEM3D simulations.
figure 1

The force generating the plane surface waves is represented only in its vertical component. The building is made of stiffer material and it is also meshed with the halfspace. Dimensions are not to scale. (b) Same as (a) but here we see the pillars forming the lens that enclose and protect the building. (c) The velocity profile as function of the radius is depicted in blue for equation (1) and in gray for the effective velocity of each cell. The inset shows a zoomed view of the cylindrical pillar and the cell.

Full size image

Luneburg Lens for Seismic Surface Waves

In his seminal work, Luneburg46 derived a spherical optical lens with radially varying refractive index that focused a beam of parallel rays to a point at the opposite side of the lens; a two dimensional variant is straightforward to deduce. Of course this relies upon the governing equation being the Helmholtz equation, which the full elastic equations certainly are not. However, it was realised in elasticity many years ago55 that one could deduce governing equations specifically for surface waves and that asymptotically, for Rayleigh waves, a governing wave equation upon the surface can be deduced whose wavespeed is the Rayleigh wave56. Using this insight we model the Rayleigh wave solution as obeying a surface Helmholtz equation and then utilise the usual transformation optics approach. Thus the concepts advanced by Luneburg translate unchanged to a plane wavefront of seismic waves travelling on the Earth surface through a lens to a point located on the boundary of the lens. One can also think of this using a ray-theory approach and note that an asymptotic surface wave ray theory is well-developed and used in seismology57.

In the model configuration presented here the elastic energy is primarily carried by Rayleigh surface waves; they are a particular solution of Navier’s equation for elastodynamics for a half-space bounded by a traction-free surface, e.g. the Earth’s surface. Well known in seismology, for the idealised situation of isotropic and homogeneous media Rayleigh waves are non-dispersive, elliptically polarized and in practical terms25 they have a velocity very close to that of shear waves: where μ is the shear modulus and ρ the density58 so for simplicity we will simply use the shear wave speed in our analysis. Shear horizontally polarized waves (SH) are also present in our numerical model, and they also propagate with wavespeed vs; notably SH waves are governed by a Helmholtz equation without any approximation. We do not consider Love waves here, which can also be important is seismology, as they only exist for stratified layered media and we assume that our elastic half space is vertically homogeneous, that is, the material parameters do not vary with depth. In Cartesian coordinates we take z to be the depth coordinate and x, y to be in the plane of the surface, then the Rayleigh waves can be represented using a Helmholtz equation on the surface and we consider a circular lens on the x-y plane as in Fig. 1c, is characterized by a radially varying refraction profile47. This lens, and the associated material variation, then extends downwards and the material is considered vertically homogeneous; we distinguish the material outside the lens to have parameters with a subscript 0 and that inside to have subscript 1.

The refraction index n between two media, say, material 0 and material 1 can be formulated in terms of the ratio of velocity contrast . For a Luneburg lens we require the refractive index, n(r), to be:

where is the radial coordinate and R the outer radius of the lens (Fig. 1c). We tune the material velocity within the lens to reproduce the index given in 1 so

Taking a continual material variation is perfect for theory, but from a practical perspective it is not possible to realize a circular structure 10’s of meters in depth and radius, whose soil properties change smoothly (e.g. on the scale of Fig. 1c). Instead we create a composite soil made of bimaterial cells such that their effective material properties have the variation we desire, this provides a realistic lens using actual soil parameters that could be created using conventional geotechnical techniques54,59.

In Fig. 1c the circular surface of the lens is discretized using equally spaced cells on a periodic square lattice. Each cell contains an inclusion of softer material that, in our illustration, is represented by a pillar extending down into the soil; the exponential decay of the Rayleigh wave amplitude with depth means that for the computational model we can truncate this and a depth of 30 m is more than sufficient. The diameter of each pillar is determined using the effective velocity prescribed for each cell based upon its radial position (r) from the center of the lens. Assuming a square section cell of width l on the x-y plane the filling fraction is defined using the surface area occupied by the pillar in the cell. For cylindrical pillars with diameter d (Fig. 1c) we have a geometrical filling fraction, f, with . The Maxwell-Garnett formula60,61, which we take for two-dimensional planar composites with cylindrical inclusions thereby consistent with our viewpoint of the surface being governed by a surface wave equation, relates the filling fraction with the corresponding effective property:

where vse is the effective shear velocity in the cell and vi is the shear velocity of the inclusion (the pillar). We combine the geometrical definition of f with (3) to obtain the effective velocity as a function of inclusion size. Hence, by tuning the pillar diameter we obtain the required velocity variation desired in Eq. 2 and use this to define the structure and variation for each of the Luneburg lenses one of which is shown in (Fig. 1c).

We now place four Luneburg lenses as shown in Fig. 1b and use these to protect an object placed in the space between them. The idea is simply that a plane wave incident along either the x or y axes will be focussed by the lens to a single point, the point at which the cylinder touches its neighbour, which will then act as source into the next Luneburg lens and the plane wave will then re-emerge unscathed; the building to be protected should, in this perfect scheme, be untouched. We are aiming to demonstrate the concept not in a perfect scenario, but using realistic parameters and a setting in which the effective medium approach provide a discrete velocity profile, yet the protection achieved is considerable.

To construct the Luneburg lenses, to reach the minimum vs prescribed in Eq. 1, vi needs be lower than 350 m/s. We choose a vi of 200 m/s which is a value that is realistic for poorly consolidated soil (sand or water filled sediments)62,63. In the lens configuration depicted in Fig. 1b,c for each lens there are 26 elementary cells (~6 × 6 m) along the radial axis of the lens and the diameter of the pillars increases towards the center of the lens as discussed earlier. In the frequency range we investigate (3–8 Hz), the inclusion is deeply subwavelength and non-resonant. The only parameter of interest for the lens design using composite soil is the filling fraction, so there is no bound on the size of the elementary cell so long as it remains subwavelength to avoid Bragg scattering phenomena; in our simulations the minimum size of each cell is bounded for numerical reasons (explained in the next section). For an actual implementation of the lens the cell could be chosen to be even smaller than our choice here with corresponding decrease in the pillar diameters. A maximum diameter of approximately 2 m would permit the pillars to be realised with existing geotechnical machineries54,59.

Numerical Implementation of The Luneburg Lens

Three dimensional numerical simulations of seismic elastic surface waves are implemented using SPECFEM3D: a parallelized, time domain, and spectral element solver for 3D elastodynamic problems widely used in the seismology community64.

The reference computational domain, depicted in Fig. 1a, is a 40 m depth halfspace 350 × 500 m wide of homogeneous sedimentary material with background velocity vs set to 400 m/s. The computational elastic region, apart from the surface itself, is surrounded by perfectly matched layers (PMLs)65 that mimic infinity and prevent unwanted reflections from the computational boundaries and these are standard in elastic wave simulation66.

A small building with a flexural fundamental mode of approximately 4 Hz is located at the center of the model. We place a line of broadband Ricker source time functions centered at 5 Hz to generate an almost perfectly plane wavefront (Fig. 1a,b). The driving force F has equal components in all 3 orthogonal directions and hence the resulting excitation is not only made of Rayleigh but also of SH waves. Body waves leave the computational domain and pass into the PML at the bottom and side boundaries of the computational domain; their interaction with the structure and the lens is negligible. This configuration could represent the upper layer of a deep sedimentary basin with some strategic structure located at the surface (power plants, data centers, hospitals) desired to be shielded from a seismic energy source. In Fig. 1b representing the building surrounded by the square array of lenses, the pillars are inserted as softer inclusions in the background velocity model. This approach is commonly used67,68 to simplify the discretization of complex models. We use very fine meshing (down to 2 m) and this, combined with the 5th order accuracy in space of the spectral element method, allows us to accurately model pillars down to a smallest diameter of 0.3 m. This approach was validated against a model where the pillars were meshed explicitly; the only notable difference was a factor of 5 increase in the runtime for the explicit case vis-a-vis the regular mesh. SPECFEM3D is a parallel code and simulations are run on 64 cores for a total of approx. 30 corehours for a simulated time of 1.5 seconds with the 3D wavefield saved for post-processing; SPECFEM3D is the standard geophysics code used in academic and industry applications with a long history of development and application64.

Results

This simulation is shown for wave-field snapshots for different propagation distances in Fig. 2 both with, and without the lenses. The sources generating the plane wavefront in Fig. 1 are located at the surface and so most of the seismic energy propagates as Rayleigh and SH waves. The vertical component of the displacement u shown in Fig. 2, is dominated by the elliptically polarized motion of the Rayleigh waves. Although not visible, SH waves behave very similarly to Rayleigh waves for the model here discussed, body waves have far lower amplitude and are not relevant to our analysis. Figure 2a shows, as one would expect, the strong scattering by the building and its highly oscillatory motion. When the Luneburg lenses are inserted in the model (Fig. 2b) the simulation shows that the Rayleigh wave front splits and then progressively converges to the focal points of lenses L1 and L2. Given the square layout, the focal points lie exactly at the touching points of the L1–L3 and L2–L4 lenses. This lens does not support any subwavelength phenomena (the evanescent part is lost) hence, the size of the focal spot is diffraction limited at λ/2n and some energy is backscattered during the focusing process creating some low amplitude reflections. The second lenses (L3 and L4) behave in a reciprocal manner converting the two point (secondary) source-like wavefields back into a plane wavefront. During the entire process, the inner region where the building is placed has experienced a very low seismic excitation as compared to the reference unprotected case. Figure 3 presents the motion of the roof of the reference building on a dB scale and it shows the vibration is drastically reduced. The snapshots in the bottom row of Fig. 2 showing the wavefront as it emerges from the lenses shows that despite the strong alteration of the ray-path, the reconstruction of the wavefront after the lenses is surprisingly good. To validate the wavefront reconstruction, Fig. 2 shows equal-propagation-distance snapshots rather than equal-time to compensate the different propagation speed in the reference and in the protected case. Hence this device combines the some cloaking behaviour with the seismic protection. Considering the broad frequency bandwidth of the input signal this is an interesting result as most cloaks so far proposed have problems with broadband excitation. The velocity structure of the lenses is such that the propagation of the wavefront in Fig. 2b is slightly slower than the reference configuration of Fig. 2a. Thus, we observe cloaking functionality to be valid for the wave envelope but not for the phase. This is not particularly relevant in the present seismic context where the only application of cloaking is to avoid very directive scattering, it would be unfortunate and undesirable to scatter or refocus the signal to a neighbouring building, while simultaneously realising seismic protection.

Figure 2
figure 2

(a) Snapshot of the vertical component of the wavefield for the homogeneous case for different propagation distances from the source (approximated values). The source position, from which distance is measured, is shown in the top snapshots. The colorscale is in normalised unit and saturated. (b) Same as (a) but with the lenses present. Full video available as supplementary material.

Full size image
Figure 3
figure 3

(a) Maps of the elastic energy distribution at the surface (z = 0) for the homogenous case. Same as (b) but for the lenses case. (c) Spectral density of the rooftop motion of the building in dB. In the upper panel, the blue trace is calculated for the homogeneous case while the black is obtained when the lenses enclose the building. The lower panel is identical, but now for the vertical component.

Full size image

A quantitative analysis of the wavefield is presented in Fig. 3a,b which show the energy maps for the reference and protected cases. The energy is calculated at the surface (z = 0) taking the L2 norm of the three components of the displacement field u(x, y, 0). In the homogeneous case (Fig. 3a) of the unprotected building the energy is almost uniform across the whole computational domain making the building resonate. In the protected case, the energy is focused towards the axes of symmetry of the lenses, leaving a central region relatively undisturbed; in the two stripes shown in Fig. 3b the energy (and equally the amplitude) is much higher than elsewhere. The motion of the Rayleigh and SH waves is maximal at the surface (z = 0) and decays exponentially with depth. Given their long wavelength (50–100 m) relative to the pillar length (30 m), the energy distribution will also decay exponentially with depth but the energy flow pattern of Fig. 3b (which is normalised in each depth cross-section) will remain identical for the whole pillar length.

Strong resonances in buildings are typically due to waves amplified by the underlying soft sedimentary soil (characterised by vs ≤ 500 m/s) in a frequency range overlapping the horizontal resonant modes of the structure. In the extreme case of an earthquake this phenomenon may lead to structural failures with well known dramatic consequences. Damping of strong ground motion and soil-structure decoupling are therefore two key concepts in civil engineering. The rooftop horizontal displacement (roof drift)69 is a diagnostic of the amplification phenomena due to the resonance frequency of the structure. Since this device is aimed at reducing anthropic noise in vibration-sensitive infrastructures rather than offering earthquake protection, also the vertical component of the motion (uz) must be considered. Vibrations produced by human activity are characterised by both in-plane and out-of-plane waves (SH, Rayleigh or higher modes) e.g.70. Figure 3c shows the frequency response function for the horizontal and vertical components of the motion recorded at the top of the building with, and without, the lenses. Over the whole spectrum an average amplitude reduction of 6 dB is achieved which is reduction of almost an order of magnitude in the vibration. While the frequency response in the vertical component is almost flat, the horizontal direction depicts more complex dynamic characterised by the fundamental horizontal mode of the building occurring at ~4.5 Hz. The vertical mode occurs at much higher frequency and it is not excited by the source. Complete cancellation of the wavefield is not achieved primarily because the evanescent field slightly couples with the building and as we focus the wavefield to a point source we introduce some back scattering that also interacts with the building (Fig. 2b). Nonetheless the concept is demonstrated and shows that one can successfully translate ideas from optics into the field of elastic seismic waves; Figs 2 and 3 should inspire the use of these concepts to both reroute surface waves and to reduce their impact on surface structures.

Discussion

We have combined concepts from transformation optics and plasmonics, composite media, and elastic wave theory to create an arrangement of seismic Luneburg lenses that can reroute and reconstruct plane seismic surface waves around a region that must remain protected. The lens is made with a composite soil obtained with columns of softer soil material with varying diameter distributed on a regular lattice. The use of softer soil inclusions emphasises that the methodology we propose is not reflection of waves, or absorption, or damping for which rigid or viscoelastic columns might be more intuitive; the softer inclusions are designed to progressively alter the material itself so that waves are “steered” and the reconstruction of the wavefronts after exiting the arrangement illustrates the mechanism. The Luneberg lens arrangement proposed here could be tested in a small scale experiment using elastic materials, such as metals, as the concept itself is very versatile or on larger scale test areas where one could then evaluate effects such as nonlinearity. Although presented in the context of seismic engineering there are everyday ground vibration topics that could benefit from this design. The damping of anthropic vibration sources (e.g. train-lines, subway, heavy industry) is very important for high precision manufacturing process, to reduce structural damage due to fatigue71 or simply to decrease domestic or commercial building vibrations.

Our aim here has been to present the concept of steering elastic surface waves completely in context and with a design that could be built, this should motivate experiments, further design and the implementation of ideas now widely appreciated in electromagnetism and acoustics to this field. One important practical point regarding our design is that we have presented normal incidence to the four Luneberg lens and this is practical, of course, if the position of the vibration (a railway line for instance) is known. Figure 4a shows the average differential energy between reference and protected case calculated underneath the building for various incidence angles of the wavefront. The protection progressively deteriorates as the incidence angle of the plane wave increases; at 45°, the focal point is in the centre of the region (Fig. 4b) and the energy is steered to this point.

Figure 4
figure 4

(a) Plot of the differential energy between reference and protected case vs. incidence angle of the wavefront with respect to the lenses. (b) Same as Fig. 3b but for a plane wave approaching the lenses with a 45° incidence angle.

Full size image

The concept of using soil mediation to steer surface elastic waves is clear and the Luneburg arrangement is a clear exemplar of the ideas. The proposed design is easily adapted to higher frequency bands and smaller regions. If one is only interested in seismic protection, and less in the wavefront reconstruction, the lenses can be structured differently with only one or two lenses. The other well-known GRIN lenses offer different extents and types of wave control (e.g. Maxwell, Eaton72) can provide an isotropic wave shielding. Other types of GRIN lenses and layout (for instance using 4 half-Maxwell lenses as shown in optics47) can be utilised; however the Luneburg lens requires the lowest velocity contrast between the lens and exterior region and we choose to use it as practically it could be implemented. The main practical difficulty is that these lenses are either singular (an Eaton lens has vs = 0 m/s in the center) or prescribe stronger velocity contrast (Maxwell) requiring difficult (or not yet available) soil engineering solutions.

Methods

The propagation of seismic waves in a 3D halfspace is a well-known problem in numerical seismology and modeled applying PMLs condition on all boundaries but for the top-surface that is traction-free. The accuracy of the method has been thoroughly tested using plate and rods as input model and it has delivered excellent results. The 3D time domain simulations are carried out using SPECFEM3D a code that solves the elastic wave equation using finite difference in time and the spectral element method in space. The parallelization is implemented through domain decomposition with MPI. The mesh is made of hexahedra elements and it is generated using the commercial software CUBIT. Simulations are then run on a parallel cluster (Curie at TGCC Paris) on 64 CPUs. 3D plots and video have been generated with Paraview and Matplotlib.

Additional Information

How to cite this article: Colombi, A. et al. Transformation seismology: composite soil lenses for steering surface elastic Rayleigh waves. Sci. Rep. 6, 25320; doi: 10.1038/srep25320 (2016).