Large scale mechanical metamaterials as seismic shields

Earthquakes represent one of the most catastrophic natural events affecting mankind. At present, a universally accepted risk mitigation strategy for seismic events remains to be proposed. Most approaches are based on vibration isolation of structures rather than on the remote shielding of incoming waves. In this work, we propose a novel approach to the problem and discuss the feasibility of a passive isolation strategy for seismic waves based on large-scale mechanical metamaterials, including for the first time numerical analysis of both surface and guided waves, soil dissipation effects, and adopting a full 3D simulations. The study focuses on realistic structures that can be effective in frequency ranges of interest for seismic waves, and optimal design criteria are provided, exploring different metamaterial configurations, combining phononic crystals and locally resonant structures and different ranges of mechanical properties. Dispersion analysis and full-scale 3D transient wave transmission simulations are carried out on finite size systems to assess the seismic wave amplitude attenuation in realistic conditions. Results reveal that both surface and bulk seismic waves can be considerably attenuated, making this strategy viable for the protection of civil structures against seismic risk. The proposed remote shielding approach could open up new perspectives in the field of seismology and in related areas of low-frequency vibration damping or blast protection.


Introduction
Of all the possible natural hazards, earthquakes are among the most catastrophic in terms of human, socioeconomic and environmental impacts. Every year more than a million earthquakes (roughly two earthquakes per minute) occur worldwide, accounting for nearly 60% of all disaster-related mortality [1,2]. Traditional seismic isolation systems aim at extending the lifetime of protected structures by means of various passive, active and hybrid control techniques [3,4]. In general, these systems are inefficient for large earthquakes and cannot be adapted to structural changes [5]. In addition, they produce dangerously large horizontal displacements [6] and ignore soil-foundation interactions that play a key role in the overall earthquake response of buildings [7]. The dynamic behaviour of a structure embedded in the ground is included in a recently-proposed type of foundation with incorporated vibrating inclusions, which can significantly reduce seismic wave energy in certain frequency ranges [5]. However, newly built foundations cannot be used to protect existing buildings of civil, historical, cultural or economic importance. The attenuation of seismic waves before they reach critical targets would be a largely preferable strategy, additionally providing the means to protect distributed areas rather than individual structures. This approach can be implemented via seismic wave barriers made of mechanical metamaterials (phononic crystals and locally resonant metamaterials), which provide wave manipulation possibilities, exploiting their unconventional properties such as negative refraction [8,9], cloaking [10], frequency band gaps (BGs) [11], etc. From a historical point of view, the concept of metamaterials was derived from electromagnetism [12] in which properties associated to negative permittivity and permeability were first observed. Following this approach, Wu et al [13] developed an effective medium theory for elastic metamaterials showing that effective parameters (bulk modulus, shear modulus and mass density) can become negative near resonances by choosing appropriate resonant scatterers, leading to eight possible types of wave propagation. In the following years, drawing inspiration from double-negative electromagnetic materials, an elastic metamaterial comprising fluid-solid composite inclusions was proposed [14], with both negative shear modulus and negative mass density over a large frequency region, and the unique property of only allowing transverse waves to propagate with negative dispersion, while forbidding longitudinal waves. Finally, an innovative metamaterial was proposed (a so-called 'hybrid elastic solid'), exhibiting multiple resonances in its building blocks and allowing band structures to display two negative dispersion bands [15]. This material exhibits a region supporting only compressional waves and another displaying 'super anisotropy' in which compressional waves and shear waves can propagate only along different directions. Among the aforementioned properties of metamaterials, in the present work we exploit their ability to exhibit frequency BGs, i.e. frequency ranges within which wave propagation is inhibited regardless of the incidence angle of an incoming wave.
Seismic waves are generally divided into body and surface waves, both of which are a superposition of longitudinal and shear bulk waves [16]. Data relative to different earthquakes show that the main frequency components of these waves span the range of 1-20 Hz [5,17]. Surface waves travel slower than body waves, with exponentially decaying amplitudes into the depth, and cause surface ground motion that can be destructive for buildings [18]. The first attempts to apply mechanical metamaterials for seismic surface waves were performed experimentally, showing the presence of BGs in scaled marble quarries with cylindrical holes at kHz frequencies [19,20] and numerically for scarcely realistic structures with kilometre-sized cylindrical holes [21]. Later, it was proposed to attenuate the shear component of surface waves by using gigantic chiral locally resonant metamaterials around isolated buildings [22].
Seismic body waves are usually considered less destructive than surface waves. However, for specific geological configurations, such as a deep layer with a smaller shear-wave velocity compared to that in overlying layers, a substantial part of the seismic energy can be channelled in the stiffer layer and propagate in the form of guided or Lamb waves [23]. Recent experiments on real-size and scaled metamaterials [24][25][26] consisting of arrays of cylindrical bore holes or tubes with local resonators, have demonstrated the presence of BGs for surface waves around 50-100 Hz, which remain above the most destructive frequencies of any earthquake. The theoretical and numerical modelling in all of these studies was performed separately for surface and bulk waves and only for simplified two-dimensional-models, whereas seismic waves arise from the superposition of these waves and are essentially three-dimensional (3D) . In particular, 2D modelling may lead to BG width overestimation in the case of locally resonant metamaterials [26] or phononic crystals [27], since coupling with out-of-plane modes is ignored. Furthermore, wave dissipation in the ground is usually completely ignored, while it may play a key role in the earthquake response of structures and should be accounted for in the design of real seismic shields. In addition, the results obtained for scaled mechanical metamaterials are not fully transferable to real earthquake scales [26], while the thus far proposed real-size structures appear to be unfeasible and difficult to manufacture by means of existing technologies due to the huge sizes involved. Finally, no attempt has been made thus far to compare the shielding performance of phononic crystals with that of locally resonant metamaterials.
In this paper, we propose and numerically analyse 3D Large-Scale Mechanical MetaMaterials (LSM 3 ) for the shielding of seismic waves propagating in dissipative soils. We perform a detailed investigation of the influence of geometric and mechanical parameters on the attenuation potential of feasible phononic crystal and locally resonant metamaterial configurations in typical frequency and intensity ranges for seismic waves. To do this, finite-element modal analysis and dynamic transient simulations are performed for both surface and guided seismic waves. The efficiency of the proposed metamaterial structures to protect sensitive sites are evaluated through full-scale numerical experiments and practical guidelines are provided.

Methods
The proposed strategy aims at attenuating typical frequencies of seismic waves via ad-hoc designed LSM 3 s. The metamaterial configurations are realized by means of cavities, also called boreholes [24,25], and periodic inclusions in the soil in a square array surrounding the structure to be protected, as schematically shown in figures 1(a)-(c) (a representative unit cell of the corresponding reciprocal space is indicated in figure 1(d)). Three different architectures are considered for the periodic unit cells representative of LSM 3 structures: (i) a cross-like cavity (figure 1(e)), proven to be more effective in inducing large BGs [28] compared to other geometries with the same cavity volume (see supplementary material); (ii) a hollow cylinder (figure 1(f)), made of a stiffer material with respect to the surrounding material, filled with soil, chosen for its fabrication simplicity [29]; (iii) a locally resonant inclusion ( figure 1(g)), made of a soft rubber layer around a heavy core cylinder, chosen due to its ability to generate subwavelength BGs compared to characteristic seismic wavelengths [30].
The corresponding geometrical parameters are summarised in table 1. We consider a sandy-type soil matrix with (Young's modulus . The soil is at first modelled as a linear elastic isotropic material, while the effects of intrinsic viscoelastic losses are considered later. Note that though the introduced elastic parameters characterize typical sandy soils, results can be generalised for other ground materials by scaling the BG frequencies and unit cell dimensions according to the wave velocity value. The sensitivity of the results to mechanical properties of the soil is studied in detail in the next section. k k , , x y the frequencies w are then determined by solving an appropriate eigenvalue problem that allows the construction of dispersion curves w ( ) k . Performing standard finite-element discretization, each unit cell is meshed by means of 4-node Lagrange tetrahedral linear elements. Good convergence is obtained considering at least 8 elements for the shortest wavelength associated to the highest frequency mode reported (for instance, a maximum finite element size = L 1 m FE is used for the case of cross-like holes considering a maximum frequency of 6 Hz). In the case of simulations for viscoelastic soils, the dispersion relations are calculated by means of an inhouse Matlab-based code capable of determining complex-valued solutions of eigenvalue problems. Since the viscoelastic stress-strain relation is usually represented in terms of hereditary integrals [16], the wave propagation problem in a viscoelastic medium is essentially nonlinear, and the dispersion relation cannot be derived straightforwardly as in the elastic case. However, the problem can be solved in the frequency-domain by employing the classic elastic-viscoelastic correspondence principle [16,35]. According to the latter, the solution of the transformed viscoelastic equations is obtained in the same manner as for the corresponding elastic case (i.e. for a representative unit cell of the same configuration and dimensions, subject to the same initial and boundary conditions as the viscoelastic structure) by replacing the material constants with their frequencydependent counterparts. The bulk modulus ( ) K x is thus replaced by To model physically realistic scenarios, we assume the frequency ω to be a real-valued independent variable, and the wave vectors k , x k y to be complex-valued vectors representing spatial wave attenuation. Thus, the dispersion curves are calculated for specified values of ω by solving the complex-valued eigenproblem for k x or k y [36]. The boundary conditions and mesh resolution are the same as for the linear elastic structures. This in-house code is validated by calculating the dispersion spectra for zero viscosity values and comparing results for propagating modes with those for corresponding linear elastic structures evaluated by COMSOL Multiphysics.
Finally, time-transient 3D FEM analysis is performed using the commercial code Abaqus©. A mesh with approximately 4 000 000 linear hexahedral elements of type C3D8R is used with software integrated hourglass control for time integration in order to guarantee solution accuracy.

Lamb waves
The seismic wave shielding potential of the proposed LSM 3 is first evaluated through dispersion analysis of guided Lamb waves. Dispersion diagrams for the three considered unit cells are shown in figure 2 for different height to width ratios (H/a=1, 2, 3). In this case free-free boundary conditions are considered on the top and the bottom of the unit cell in order to reproduce the scenario of quasi-Lamb waves observed in some seismic events [37,38] (when for instance one portion of the soil is in contact with another portion or half space that is less stiff, or if the soil is bounded by a deep layer filled with gas, or in the presence of a horizontal fracture within a consolidated rock, in correspondence with which the stiffness virtually vanishes [38]). Guided waves usually dominate in seismic records of local and regional events and may propagate over thousands kilometres [34]. Table 1. Case studies and corresponding geometric parameters used in the calculations. In-plane parameters for the cross-like holes and circular inclusions are taken from [28] and [32], respectively. Shielding potential in specific frequency ranges according to I/I 0 =exp[−(N/ N C ) β ] is also provided for the first two geometries.

Parameter
Shielding potential (β; N C ) Case study For cross-like cavities (figure 1(e)), the height of the unit cells plays a fundamental role in the nucleation of the BGs. In general, as / H a increases, the BG size is reduced, and more modes appear in the same frequency range. The physical reason for this is the increasing appearance of higher order slab modes as the ratio / H a increases [39]. Thus, for small values of / H a, the BG is opened due to the 'folding' of the three lowest bands [39]. For larger / H a ratios, higher-order modes appear with the cut-off frequency decreasing with the increase of / H a and finally close the BGs. Moreover, the coupling of 2D in-plane with out-of-plane modes results in the splitting of a wide BG predicted by 2D simulations into two separate BGs with mid frequencies around 2.5 and 4 Hz. Indeed, the broad BG in the left diagram (blue curves) for / = H a 1 is divided into 2 parts by the 7th pass band, whose mode shape is shown in figure 2(d) for the high-symmetry X point at = f 3.28 Hz. The corresponding motion mainly originates from the out of plane bending of a square block and is similar to a flexural vibration of a clamped beam. The adjacent pair of blocks in a unit cell vibrates in anti-phase, whereas the diagonal pair vibrates in phase. By increasing / H a, the rotational component of the block motion becomes more conspicuous, since the rigidity of the narrow connectors decreases, whereas, in the mode located above the BG with the mode shape reported in figure 2 4.58 Hz , flexural motion dominates. For the inclusion-based phononic unit cell shown in figure 1(f), the unit cell height influences the inhibited frequencies less compared to the previous case, at least for the lowest BG. As the / H a ratio increases, additional slab modes appear, but all of them are located above 7 Hz. Hence, the size of the lowest BG is preserved, while the higher BG, in general, decreases in size. Note that for / = H a 2 the second BG is even larger than for / = H a 1, and its bounds are related to the modes localised in the inclusion, e.g. as shown in figure 2(g).
In the case of the locally resonant inclusion shown in figure 1(g), the chosen thickness of the rubber coating provides the best compromise between BG width and sufficiently low BG frequency (a parametric study relative to the influence of the coating thickness is provided in the supplementary material). Only one absolute BG is found in the considered frequency range, and a dependence of the BG size on the unit cell height is also observed. In this case, no additional slab modes are generated as the ratio / H a increases, and instead, the BG size decreases. This can be explained by considering the mode shape at the lower bound of this BG, the frequency of which increases with the increase of / H a. As shown in figure 2(h), this mode shape corresponds to a localised mode with mixed in-plane and out-of-plane displacements. With the increase of height, the coupling of in-plane and out-of-plane motions occurs at higher frequencies, finally leading to the closing of the BG. Vibration patterns for other modes are also characterised by localised motions within the inclusion (e.g., see figure 2(i)).

Surface waves
We next investigate the performance of the two chosen structures in attenuating surface waves by checking whether the previously found BGs for guided waves are preserved in a half-space structure. The solutions describing waves in a half-space can be found by solving an equivalent PML problem, when the PML boundary condition is applied to the bottom of the unit cell. The dispersion diagrams of the surface guided waves for both unit cell types are shown in figure 3 for different height to width ratios. The radiative region, or sound cone, is indicated by the shaded violet area [40]. The boundary of this region is formed by the slowest bulk wave propagating in the soil [40]. Therefore, propagating surface waves, the velocity of which is slower than that of bulk waves [16,41], are located outside the cone, while bulk modes, leaky surface modes as well as spurious (unphysical) PML modes can be found within the cone [42]. For the unit cell with a cross-like cavity, two or three BGs for surface waves exist in the corresponding band structure for various values of the unit cell height. These BGs are located at almost the same frequencies as those predicted for the Lamb waves. Again, as for guided waves, an increase of the LSM 3 height leads to a decrease in BG size (figure 3). The mode shapes for propagating waves are shown in figures 3(d) and (e) indicating the localisation of motion near the free surface of the unit cell, as expected for surface modes.
For phononic structures with hollow cylinders, the BGs predicted for guided waves are completely preserved in the dispersion diagrams shown in figure 3(b). However, the BGs are located inside the sound cone and can only inhibit the propagation of bulk or surface leaky waves. Hence, phononic-type metamaterials appear to be inefficient for shielding propagating surface seismic waves, as was also found in previous works [24,25].
Thus, among the three considered LSM 3 configurations, it is mainly the one with cross-like cavities that is capable of efficiently attenuating destructive surface waves. Since the cavities are dug in sandy soil, the problem of the technical feasibility and stability of the cavities may arise. However, this problem could be addressed by containing the soil in a cavity boundary made of a thin layer of stiffer material, e.g. aluminium. This additional layer results in a shift the BGs to slightly higher frequencies, but ensures sufficient soil compaction to avoid failure due to flexural-torsional deformations (see supplementary material). Additionally, the presented analysis remains valid for wave dispersion in stiffer soils, for which stability problems are less relevant, with appropriate scaling to higher frequencies. 3 containing an array of LSM 3 s. The shielding region is formed by an array of unit cells of cross-like cavities or hollow cylinders with = = H a 10 m, arranged in a square grid around the area to be protected, as shown in figure 4(a). The unit cells are arranged in four concentric square rings around an area of 10×10 m 2 . The mechanical constants for the soil are listed in the previous section. In this case, soil damping is not accounted for, to evaluate the attenuation capability of the introduced periodic structure only. Since unit cells with cross-like cavities can attenuate both bulk and surface waves, as shown in the previous section, here we present results for this geometry only (see supplementary material for the case of hollow cylinders).

Time transient analysis
To assess seismic shield performance, two different input signals are considered, both representative of a seismic event with a multi-frequency content centred in the range of interest for the protected structures. The two pulse shapes and their frequency contents are shown in figure 4 (d), respectively. This analysis shows that ground vibration is efficiently attenuated for surface waves with a frequency content falling inside the BGs, and therefore the seismic shield can provide a virtually unperturbed area, even if reflections at the free boundaries of the volume are included. Furthermore, the amplitude of the incoming wave rapidly decays on the surface and through the depth after the first 2-3 rows of the LSM 3 inclusions. Moreover, the LSM 3 deviates the wave stress field below the shielded area at an angle that is sufficient to protect the building foundations ( figure 4(d)). Results show that only four LSM 3 rows in the shield are sufficient to reliably reflect the incoming wave energy, effectively limiting the impact of seismic waves on the protected area. Time transient displacements in the z direction are recorded at the two acquisition points A 1 and A 2 for broadband and narrowband excitations (upper boxes of figures 4(e) and (f)), to allow quantitative comparison between the motion registered in an ordinary portion of the soil (point A 1 ) and inside the designed LSM 3 region (point A 2 ), respectively. After acquisition, signals are Fourier transformed and their frequency content is compared considering both actuation pulses (lower boxes of figures 4(e) and (f)). BGs predicted for Lamb waves in the infinite system are also highlighted as shaded grey rectangles. Results show that the introduction of the LSM 3 region reduces the ground motion considerably. In the case of the broadband excitation, the peak displacement is decreased by approximately one order of magnitude with only 4 rows (from 3.32 down to 0.34-displacements are normalised to the maximum positive value of the displacement registered inside the LSM 3 ). Analysis of the frequency content of the signals shows that in the BG regions the spectrum of the response at A 2 is reduced by a factor varying between 10 and 43 compared to that of signals taken at A 1 . Similar results are obtained when the depth of the source of excitation is varied, with an increasing reduction in barrier efficiency for increasing impinging angles (see supplementary material). However, the surface source hypothesis remains valid for surface waves that decay more slowly and are the most destructive among seismic waves.
The effect of the proposed seismic shields is further verified including model structures in the simulations. One example is a cultural heritage site such as the 'Castello delle Quattro Torra', near Siena, modelled schematically in the FEM simulation. Figure 5(a) shows the normalised Von Mises stress field map (blue and red colours represent minimum and maximum stress, respectively) for a shielded and non-shielded structure at time step t=5.3 s for a 5 Hz centred wave (19 modulated sine cycles). The non-shielded structure clearly exhibits very high stress and deformation levels, inevitably leading to its collapse, whilst the shielded one is subjected to sufficiently attenuated deformation levels to guarantee its structural integrity (see supplementary material).
The screening efficiency of the proposed LSM 3 s is also investigated for different numbers of unit cell rows N, and the risk mitigation capability is quantified relative to the MSK-76 Intensity scale and the associated peak values of ground displacement [43]. Numerical models of cross-like-cavity or hollow-cylinder structures with 0, 1, 2, 3 and 4 rows of LSM 3 s are considered. Figure 5 figure 5(c) for hollow-cylinder geometry LSM 3 and wave excitation centred at its corresponding BGs at 6 and 8 Hz. In all cases, the wave amplitude significantly decays as the number of unit cells increases, with an exponential-type decrease. Notably, the introduction of 4 rows of LSM 3 can lead to a seismic threat reduction between two and five orders of magnitude, from a highest devastating-event magnitude of X to a medium level earthquake of magnitude VII or a moderate seism of level V, with a decrease of the corresponding perceived shaking level. The calculated reduction as a function of LSM 3 rows can be fitted with an exponential-type phenomenological relation (also included in figure 5(a) and (b)) as N N 0 C where β and N C are fitting parameters, the latter indicating a 'critical' number of LSM 3 rows that provides a shielding factor of e. Parameter values are reported in table 1 for the corresponding LSM 3 geometry. This equation can provide an approximate rule of thumb for the design of LSM 3 s in the lowfrequency range, according to the required shielding potential and specific BG frequency range. Overall results lead to the conclusion that all considered LSM 3 geometries can usefully be exploited and adapted to the specific requirements, depending on the soil characteristics and properties (e.g. resonant frequencies) of the shielded structures.

Parametric study
Next, we analyse the effects of the elastic parameters of the unit cell constituents on the BG frequencies. This study is performed for LSM 3 comprising hollow and locally resonant inclusions. In the case of hollow cylinders, the material density r and Young's modulus E are varied, while Poisson's ratio and the cell size a remain unaltered to the reference parameters mentioned above. Simulations show that mass density does not play a significant role, while BGs are nucleated when the cylinder material is at least 250 times stiffer than the matrix. In the case of a sandy soil, this means that concrete or steel inclusions would be adequate. Figure 6 This analysis allows us to define optimal parameters for this type of geometry, i.e. 8 m diameter, 1 m thick, hollow concrete cylinders.
In the case of locally resonant inclusions, only one BG appears, which is narrower compared to the previous case. By changing the mechanical parameters, i.e. the E E rubber soil ratio, a shift of the BG central frequency occurs ( figure 6(d)). The increase of both the internal core mass density (figure 6(e)) and the R a e ratio (figure 6(f)) allows enlarging the BG size. However, the widest BG is smaller than the BGs nucleated in the case with hollow cylinders. Also, the BG is found to nucleate only above the threshold of = R a 0.47. e Thus, the design solutions using cross-like cavities or hollow cylinders prove to be preferable compared to the case of coated cylinders for seismic applications. Optimal geometric and material parameters for a square unit cell arrangement are indicated in figures 6(a)-(f) with dashed-dotted lines, considering both sufficient BG properties and ease of practical realisation of the structures. Further BG frequency enlargement may be achieved by using a triangular arrangement of the scatterers, choosing the position in accordance with the suppression of the shear potential energy of the first optical band, as demonstrated by Lai and Zhang [44] in a two-dimensional triangular lattice of aluminium cylinders embedded in an epoxy host.
In general, seismic isolation systems based on local resonance effects allow the use of much smaller unit cells (a=2 m instead of a =10 m for the other cases) to obtain BGs at approximately the same frequencies. However, the chosen configuration of the locally resonant metamaterial is characterised by very narrow inhibited frequency ranges and, thus, seems to be practically inefficient compared to two other proposed configurations. Moreover, a small unit cell size in the depth direction may be counterproductive, since the seismic waves with larger wavelengths can propagate below the locally resonant inclusions. Therefore, detailed analysis of the wave attenuation performance for LSM 3 has been limited to unit cells with cross-like cavities and hollow cylinders in this work.

Viscoelastic effects
Due to the presence of underground water, oil and gas, soil exhibits energy losses and time-dependent mechanical characteristics. These effects can be taken into account by modelling the sandy soil, at least as the first approximation, as a linear viscoelastic medium [45]. We thus analyse the influence of soil viscoelasticity on the shielding capability and the BG sizes for the proposed LSM 3 in a low-frequency range. The simulations are performed here for Lamb waves in the case of cross-like cavities with H=10 m (table 1). Material damping in the soil is presumed to increase linearly with frequency, as normally assumed when modelling wave propagation in viscoelastic metamaterials [41,46]. The shear modulus of the material G ( ve) therefore becomes complexvalued: , whilst the bulk modulus is assumed to remain unchanged from the elastic case to simplify the analysis: ]. Since viscoelastic parameters depend on specific soil structure and differ for various compositions, we choose representative values for the viscosity  G in order to perform a qualitative study of these effects (a parametric study is presented in the supplementary material).
For viscoelastic metamaterials, the difference between propagating (with real-valued wavenumbers) and evanescent (with imaginary or complex-valued wavenumbers) modes disappears due to the losses, and the concept of BGs no longer applies [47]. Therefore, the analysis is performed by solving dispersion relations between complex-valued wavenumbers = + ( ) * k k k Re i Im and real-valued frequencies f . The attenuation properties of viscoelastic metamaterials can be evaluated by analysing the so-called attenuation spectrum relating the wave attenuation level h = ( ) ( ) k k 2Im Re to the frequency [48]. The attenuation spectrum for waves propagating in the G -X direction is given in figure 7 with red lines representing the attenuation of seismic waves in the viscoelastic LSM 3 , blue and black lines showing attenuation of propagating and evanescent waves, respectively, in the corresponding elastic structure (i.e. with zero viscosity). According to the definition of h, the wave attenuation for the propagating waves in the elastic LSM 3 is zero outside the BGs and has finite values for the evanescent modes. The BGs for the linear elastic case are highlighted by shaded regions in figure 7. According to the adopted assumption, the wave attenuation increases linearly with the frequency but proves to be more pronounced for frequencies falling between the BGs, i.e. between 3.2 and 3.4 Hz (highlighted rectangle  figure 7). At these frequencies, the attenuation is comparable to that inside the lowest BG and further increases with the growth of the viscosity G″ (see supplementary material). Therefore, the viscoelastic properties of the soil improve the wave attenuation near the BG boundaries and may lead to the additional effect of 'joining' closely located BGs. Hence, we can conclude that the viscoelasticity of the ground enhances the shielding capability of the proposed LSM 3 s and results in wider BGs with respect to those predicted by the corresponding linear elastic model. The same conclusion remains valid for waves propagating in other directions and for hollow cylinder or locally resonant inclusion cases (see supplementary material). In the case of locally resonant inclusions, the viscoelastic behaviour of the matrix may lead to the decrease of the attenuation level inside the BG. However, the shielding performance of the whole structure can be improved by choosing for the soft coating a more viscous material than the matrix [46].

Discussion and conclusions
In summary, we have applied the concepts of phononic crystals and mechanical metamaterials to the field of seismic protection and demonstrated their feasibility in attenuating low-frequency bulk, plate and surface waves. The proposed systems, contrary to traditional structural foundation isolation strategies, which cause a shift in the fundamental vibrating frequency of civil engineering structures, reduce the seismic wave energy by means of BG attenuation mechanisms and prevent it from reaching the protected site. While this concept has already been proposed in a limited number of recent papers, the present work is the first realistic and comprehensive study on the topic, given that the study is performed with 3D simulations, it takes into account the layered structure of the soil, includes viscoelastic effects, and is based on both phononic crystal and locally resonant structures. Various design configurations, based on both cross-like cavities and cylindrical inclusions, have been explored and all display shielding potential in the seismic frequency range of interest of a few Hz. All analysed geometries can be constructed relatively easily and cheaply with the current technological state-of-theart. The lateral dimension of the phononic-type LSM 3 s are approximately 10 m, but reduce to only 2 m in the case of locally resonant structures, making the practical realisation of these structures possible from a technical point of view and feasible from an economic point of view when compared to the costs involved with earthquake damage. Locally resonant systems have rather narrow BGs, so that alternative design solutions may be required to broaden the frequency attenuation range, e.g. arrays of locally resonant inclusions each having different BG frequencies [26], depending on the specific site and on the structures to be protected. Parametric analysis shows that BG nucleation in locally resonant structures occurs for LSM 3 depths of the order of their lateral size (∼2 m ), and that a further increase of the depth results in decreased BG size for both surface and Lamb waves. Since these depth dimensions are small compared to wavelengths of typical seismic waves, vertical arrays of inclusions can be used down to the necessary depth to improve the attenuation efficiency of these structures. In the case of both considered cylindrical geometries, optimal stiffness and thickness values can be found to maximise the BG size. Furthermore, analysis of soil viscoelastic effects shows that any level of viscosity (from small to medium value) improves the attenuation potential of the designed structures and their extension to wider frequency ranges.
These results indicate that the proposed seismic isolation strategy is theoretically effective as well as technically feasible. For example, during the El Centro earthquake, the maximum measured displacement was of approximately 21.64 cm, corresponding to an event between IX and X degrees of severity in the MSK-76 intensity scale [43]. The application of 4 rows of hollow-cylinder geometry LSM 3 would have locally reduced the event by five degrees in the 8 Hz range, in which only slight damage to a few poorly constructed buildings would be expected. The performed 3D time-transient numerical simulations confirm results obtained from unit cell analysis and prove the shielding capability of the proposed LSM 3 s. An approximate exponential-type relation is also derived for the normalised earthquake intensity as a function of LSM 3 rows in the low-frequency regime, to provide an indicative design tool for seismic shields. The presented analysis can be exploited to provide schematic practical guidelines for the design of LSM 3 seismic shields in real cases, as follows: -Estimate the main frequencies at which risk is concentrated (e.g., the resonance frequencies of the structure); -Determine the mechanical properties of the surrounding soil and choose LSM 3 unit cell dimensions based on the resulting maximal threat wavelength; -Select optimal LSM 3 type, parameters and number of rows to obtain the desired shielding effect based on results from figure 6 and table 1.
These structures can be further optimised and adapted to the specific application of interest, e.g. to distributed as well as localised sites. Nevertheless, the present study already demonstrates the feasibility of the proposed innovative seismic isolation strategy, which could potentially lead to a considerable reduction of human casualties and socio-economic impacts.