Partially embedded metabarrier to suppress surface waves in granular media^a)

Controlling surface wave propagation is essential to protect critical civil infrastructure and urban areas from the disastrous impact of some earthquakes. Locally resonant metamaterials, owing to their unconventional material properties of negative effective mass density and stiffness attained from the hybridization of the incident wave with the local resonances of the resonators, offer unique capabilities for surface wave control by generating frequency bandgaps (Boechler , 2013; Colombi , 2014; Giraldo Guzman , 2024). Particularly, Rayleigh surface waves in a homogeneous elastic half-space primarily delocalize into shear waves that propagate into the half-space within these bandgap frequencies (Colombi , 2016; Colquitt , 2017). There has been much research on exploiting local resonant metamaterials to control Rayleigh waves, both on understanding the interaction of Rayleigh waves with local resonances (Colquitt , 2017; Pillarisetti , 2022a,b; Pu , 2021) and on proposing a wide range of metabarrier configurations for low frequency applications (Colombi , 2016; Muhammad , 2021; Liu , 2020; Miniaci , 2016; Palermo , 2016; Wu , 2024). The incorporation of soil stratification, involving multi-layered soils with varying homogeneous and elastic properties, has also been extensively explored in the design of locally resonant metamaterials for surface wave suppression (Carneiro , 2024; Pu , 2020; Zeng , 2022). However, natural soils are granular and inherently heterogeneous, necessitating the consideration of material heterogeneity in the design of metabarriers (Muhammad , 2019; Pu and Shi, 2018; Zeng , 2022).

In granular media, the elastic constants exhibit power-law dependency on the ambient pressure, which linearly scales with depth. At small strains, this dependency arises from the Hertzian contacts between the grains. The Hertzian contacts exhibit a nonlinear stress-strain relation, resulting in a power-law dependence of the elastic modulus μ on the ambient pressure p, $μ∼p1/3$⁠, with the speed of sound (c) varying as $c∼p1/6$ (Goddard, 1990; Gusev , 2006; Liu and Nagel, 1992). Unlike the non-dispersive Rayleigh waves, which are the only surface wave solutions for a homogeneous and elastic half-space, considerations of the material heterogeneity result in an infinite number of surface wave solutions called the guided surface acoustic modes (GSAMs) (Aleshin , 2007; Gusev , 2006). Gusev (2006) demonstrated the existence of a discrete set of GSAMs: PSV modes, which are fundamentally a combination of longitudinal (P) and shear-vertical (SV) modes with wave polarization in the sagittal plane and shear-horizontal (SH) modes with wave polarization perpendicular to the sagittal plane. Although there are many noteworthy studies on quantifying the material heterogeneity and attenuation behavior in granular media (e.g., Bodet , 2014, 2010; Zaccherini , 2022), only a few studies have exploited these material models in designing locally resonant metabarriers for surface wave control in granular half-space.

The hybridization between the local resonance of the resonators and GSAMs was demonstrated for the first time by Palermo (2018) in an unconsolidated granular media using a metabarrier comprising buried spring-mass resonators that exhibit longitudinal resonances. Only the first order PSV mode (i.e., PSV₁) is observed to hybridize with the resonators localizing the energy at the resonator base, in contrast to Rayleigh waves, which delocalize into shear waves propagating into the half-space (e.g., Colombi , 2016). The delocalization of the surface waves to shear waves is inhibited due to the depth dependency of the longitudinal and shear wave speeds (“mirage effect”). The higher-order modes (including PSV₂) are observed to propagate almost undisturbed, with their displacement profiles exhibiting downconversion to lower-order modes for frequencies above the anti-resonant frequencies. Zaccherini (2020b) demonstrated the applicability of the multi-layered resonant metabarriers to attain higher attenuation compared to a single layer by confining surface wave energy at the base of the lowest layer of resonators, which attenuates significantly before returning to the surface. They also introduce graded metabarriers with increasing and decreasing resonant frequencies to control SH GSAMs (Zaccherini , 2020a), building on earlier studies of meta-wedges used to control Rayleigh waves in homogeneous media (Colombi , 2016; Colquitt , 2017). In a very recent study, Mora (2021) demonstrated the use of a non-trivially shaped resonator array, which is found to hybridize with the two lower-order (PSV) surface acoustic modes by considering the consolidation of granular media underneath the resonators. The hybridization is attributed to the hardening or consolidation of the medium under the resonators and the local resonances.

Extending the previously investigated spring-mass resonators (Palermo , 2018; Zaccherini , 2020b) to low frequencies in granular soil is impractical for suppressing PSV surface modes due to additional interactions that compromise the metabarrier efficiency, despite their success in suppressing Rayleigh waves in homogeneous soil (Muhammad , 2021; Palermo , 2016), as illustrated using a transient finite element analysis in the supplementary material. In this paper, we explore metabarriers comprising partially embedded rod-like resonators, owing to their simplistic geometry and the existence of only two types of resonator hybridizations: longitudinal and flexural. This is a significant consideration because while flexural resonance hybridizations with Rayleigh surface waves in homogeneous media are well understood, such interactions with PSV modes have been overlooked in previous metabarrier studies in heterogeneous granular media (Palermo , 2018; Zaccherini , 2020a). Understanding the resonance hybridizations of these rod-like resonators holds promise for the development of more realistic metabarriers.

However, the presence of multiple dispersive surface modes (PSV₁ and PSV₂) poses challenges in realizing surface wave bandgaps even with such a simple metabarrier design, as the bandgaps are not directly evident from performing transient analysis, as discussed in the supplementary material. This limitation arises because the plotted spectral amplitude curves, obtained through transient analysis akin to prior literature (Palermo , 2018; Zaccherini , 2020b), only show the cumulative transmission of PSV₁ and PSV₂ modes, obscuring individual PSV₁ and PSV₂ mode bandgaps due to simultaneous excitation by the point-load excitation. Only by fully understanding the distinct hybridizations between PSV₁ and PSV₂ modes with longitudinal and flexural resonances of rod-like resonators can these PSV₁ and PSV₂ mode bandgaps be revealed, which is the primary focus of this paper. This objective is addressed by relying on frequency-domain analysis to study individual interactions of surface modes (PSV₁ or PSV₂) with resonators and identify associated PSV₁ and PSV₂ bandgaps. The formation of these bandgaps is also verified experimentally on a granular testbed comprising glass beads using three-dimensional (3D)-printed rod-like resonators at higher frequencies. We also explore the impact of resonator geometry on critical resonance-hybridized surface modes that generate PSV mode bandgaps. This parametric study aids in the rational design of metabarrier configurations and identifies geometry limitations for very low frequencies. Furthermore, we investigate graded metabarriers for the first time in the context of PSV mode suppression, featuring arrays of rod-like resonators with increasing and decreasing height to facilitate wider PSV₁ and PSV₂ bandgaps through unique mode conversions between PSV modes.

The remainder of the paper is arranged as follows. Section II identifies true resonance hybridizations of the PSV₁ and PSV₂ modes with the metabarrier comprising partially embedded rod-like resonators to identify the PSV₁ and PSV₂ bandgaps. Section III presents surface wave experiments on a granular testbed of glass beads to realize these bandgaps using 3D-printed rod-like resonators. Section IV presents a parametric analysis to study the influence of the unit-cell geometry on the PSV₁ and PSV₂ bandgaps. Sections V A and V B demonstrate the implications of graded metabarriers on the propagation of PSV₁ and PSV₂ modes. Moreover, a compound metabarrier configuration comprising a combination of graded metabarriers is demonstrated in Sec. V C to exhibit superior PSV mode suppression capabilities. Finally, Sec. VI summarizes key findings from the numerical and experimental studies on the interaction of PSV modes with rod-like resonators to generate bandgaps in granular media.

Studying the interaction of the surface modes (PSV₁ and PSV₂) with longitudinal and flexural resonances of rod-like resonators could aid in designing any other non-intuitive resonator designs, as most of the earlier proposed metabarriers for Rayleigh waves were based on tailoring their longitudinal and flexural-like resonances (Muhammad , 2021; Liu , 2020; Mora , 2021). Figure 1(a) depicts the unit-cell schematic of the metabarrier comprising rod-like resonators made from concrete, partially embedded in a granular medium composed of unconsolidated glass beads. This medium conforms to the depth-dependency of wave speeds, $cL,S=γL,S(ρgz)αL,S$⁠, where c_L and c_S represent the depth-dependent longitudinal and shear wave speeds, respectively. The density (ρ = 1600 kg/m³) and material constants (γ_L = 14.4, γ_S = 6.42, α_L = 0.31, and α_S = 0.31) of the granular media are taken from previous literature (Aleshin , 2007; Palermo , 2018). The volume of granular media beneath the resonator is considered partially consolidated, as expressed by (⁠ $cL,S=γL,S(ρg(z+M/ρS))αL,S$⁠), where (M) represents the mass of resonators applied over the contact area (S), following Mora (2021). Perfectly matched layers with fixed boundary conditions (BCs) are considered at the bottom of the granular half-space to filter unfeasible wave solutions that traverse along the bottom surface. Finally, the hybridized dispersion curves are obtained by imposing Bloch–Floquet BCs on the surfaces of the granular half-space to enforce spatial periodicity and sweeping the wavenumber in the wave propagation direction (x-direction) within the irreducible Brillouin zone (Palermo , 2018).

Following a parametric study on the half-space depth (not shown here), the depth is set to 150 m, as further increasing the depth has no discernible effect on the dispersion curves. The hybridized dispersion curves for a metabarrier configuration with W = 1 m, A = 1.75 m, H_R = 10 m, and H_B = 4 m are shown for reference in Fig. 1(b). Note that the Young's modulus at the resonator footing is around 170 MPa, and the pressure exerted on the granular media due to the self-weight of resonators (∼80 kPa) is well below the allowable bearing capacity (>250 kPa) of typical granular soils (NYC Construction Codes, 2014). The color scale (⁠ $U¯z=maxHBUz/maxHHSUz$⁠) in the dispersion curves helps identify the surface modes propagating closer to the surface (⁠ $U¯z∼1$⁠) from those that propagate under the resonators (⁠ $U¯z∼0$⁠). These dispersion curves demonstrate hybridization with longitudinal and flexural resonances of the rod-like resonators, with wave structures corresponding to the first flexural resonance (F₁), first longitudinal resonance (L₁), and second flexural resonance (F₂) shown in Figs. 1(c), 1(d), and 1(e), respectively. It is important to note that the cutoff frequencies for the longitudinal and flexural resonance modes in the hybridized dispersion curves differ from the fixed-free resonance frequencies of the resonator, as the resonators are partially buried. Interestingly, we can also observe unhybridized surface modes localized below the resonator, as evident from the wave structure shown in Fig. 1(f). While these dispersion curves encompass all possible surface mode-resonator hybridizations without displaying evident bandgaps, it is crucial to identify the specific hybridized modes (longitudinal/flexural resonance) to which the incident PSV₁ and PSV₂ modes preferentially hybridize to reveal the PSV₁ and PSV₂ bandgaps.

Unlike the time-domain analysis (as detailed in the supplementary material), frequency-domain finite element analysis allows for the preferential excitation of specific guided modes of interest (Chillara , 2016), facilitating a detailed understanding of the resonator interactions with PSV₁ and PSV₂ modes independently. The schematic of the 3D finite element model used for the frequency-domain analysis is shown in Fig. 2(a). The model is divided into two “buffer” regions of length 0.5λ, a “loading” region of length 4λ, “incident” and “transmission” regions of length 6λ, and a metabarrier region comprising an array of 40 partially embedded rod-like resonators with the same dimensions as those considered earlier in the dispersion analysis, where λ represents the wavelength of the incident PSV₁ or PSV₂ mode at any desired frequency. Perfectly matched layers of thickness λ are employed on all the exterior boundaries of the model to minimize boundary reflections, with periodic BCs imposed on the lateral surfaces (in the y-direction) of the model to impose periodicity in the resonator arrangement along the lateral direction. Three-dimensional brick and tetrahedral elements with a maximum element size of $λ/10$ are used in the analysis. To ensure the excitation of pure PSV₁ or PSV₂ modes at a given frequency, we apply the corresponding stress distribution of the surface mode with depth as a body load to the “loading” region. Simulations are repeated at several selected frequencies within the frequency range (8–40 Hz) for both PSV₁ and PSV₂ mode excitations. To demonstrate the efficacy of the body load excitation approach, the wavenumber spectra computed from the complex-valued displacement data (⁠ $Ux+Uz$⁠) on the surface of the metabarrier region, in the absence of the metabarrier, are overlaid on the PSV₁ and PSV₂ mode dispersion curves for all excitation frequencies. As expected, Figs. 2(b) and 2(c) demonstrate preferential excitation of pure PSV₁ and PSV₂ modes, respectively, as the peaks in all the wavenumber spectra align with the corresponding modes across the intended frequency range.

To identify the true resonance hybridizations of PSV₁ and PSV₂ modes, the wavenumber spectra computed in the presence of the metabarrier for PSV₁ and PSV₂ excitations are overlaid on the hybridized dispersion curves within the desired frequency range (8–40 Hz), as shown in Figs. 3(a) and 3(d), respectively. The wavenumber spectrum at each excitation frequency is normalized with the maximum spectral amplitude in the wavenumber spectrum obtained without the metabarrier to eliminate discrepancies caused by variable excitation amplitudes during body load excitation. Furthermore, the transmitted PSV₁ and PSV₂ modes are analyzed in the transmission spectrum panels in Figs. 3(a) and 3(d) by normalizing the spectral amplitude of the PSV₁ and PSV₂ modes in the wavenumber spectra computed over the transmission region to the spectral amplitude of the incident surface mode (PSV₁ or PSV₂) in the wavenumber spectra computed over the incident region.

The wavenumber spectra for PSV₁ mode excitation demonstrate a strong hybridization of the incident PSV₁ mode with the longitudinal resonance (L₁) of the resonators and a weaker hybridization with the second flexural resonance (F₂) present in the same frequency range. This is evident from the pronounced peaks along the longitudinal-resonance hybridized mode for all excitation frequencies. This strong longitudinal resonance hybridization aligns with previous studies on PSV₁ mode interaction with spring-mass resonators (Palermo , 2018; Zaccherini , 2020b), though these studies did not account for flexural-like resonances. Despite potential flexural resonance hybridizations, the presence of a cutoff frequency (f_L) for the longitudinal-resonance hybridized mode initiates a bandgap for the PSV₁ mode above this cutoff frequency, as highlighted in red in Fig. 3(a) and evident in the transmission spectra. Below the cutoff frequency, the transmission spectra reveal some mode conversion of the incident PSV₁ mode to the PSV₂ mode throughout the frequency range. The out-of-plane displacement fields [Real(U_z)] before [Fig. 3(b)] and after [Fig. 3(c)] the longitudinal-resonance cutoff frequency indicate suppression of the PSV₁ mode after the cutoff frequency, primarily due to reflection of the incident PSV₁ mode.

In contrast, for PSV₂ mode excitation, we observe a strong hybridization of the incident PSV₂ mode with the second flexural resonance (F₂) and a notably weaker hybridization with the longitudinal resonance (L₁) in the same frequency range. This is evident from the wavenumber spectra showing strong peaks along the flexural resonance hybridized mode for all excitation frequencies. This unique hybridization initiates the PSV₂ mode bandgap in the absence of the hybridized flexural resonance mode, highlighted in the transmission spectra [Fig. 3(d)]. Outside the bandgap, there's a slight mode conversion from the incident PSV₂ mode to PSV₁ mode, as observed in the transmission spectra. Comparing the out-of-plane displacement fields [Real(U_z)] for frequencies outside [Fig. 3(e)] and within [Fig. 3(f)] the PSV₂ mode bandgap reveals PSV₂ mode suppression by localizing the incident PSV₂ mode into the half-space, similar to the mechanism behind bandgaps for Rayleigh waves in a homogeneous half-space (Colombi , 2016; Colquitt , 2017). Therefore, achieving a strong PSV₂ mode bandgap requires more resonators compared to the PSV₁ bandgap, where mode suppression is mainly achieved through reflection. To further validate the presented numerical findings, we demonstrate the presence of these PSV₁ and PSV₂ bandgaps in a laboratory-scale setup in Sec. III.

Surface wave experiments are conducted in a steel tank measuring 122 cm × 61 cm × 45 cm, filled with unconsolidated glass beads of 150 μm diameter up to a height of 40 cm, achieving a packing density of 1600 kg/m³, as depicted in Fig. 4(a). To minimize boundary reflections, all tank boundaries are covered with 5 mm corrugated cardboard sheets. Glass beads are poured into the tank using a mesh sieve, with the sides of the tank gently shaken at multiple fill levels, following the filling protocol described in previous studies (Bodet , 2014). Cardboard baffles are inserted 20 cm deep on the side boundaries to scatter sidewall reflections, and a flat cardboard piece is used to manually level the surface without applying excessive pressure that could otherwise cause consolidation. A low-frequency magnetic shaker (Vibration Test Systems, Aurora, OH) excites a broadband Ricker waveform centered at 700 Hz into the granular medium, aided by a 150 mm long, 4.76 mm diameter stainless steel connecting rod buried 2 mm within the granular medium, positioned approximately 2 cm from one tank edge at a 35° inclination from vertical, as illustrated in Fig. 4(a). Particle velocity on the granular media surface is measured using a laser Doppler vibrometer (Polytec OFV-525; Polytec, Baden-Württemberg, Germany) mounted on a 25 cm scan length micro-precision scanning stage (Newport ILS250PP; Newport, Irvine, CA). The Newport XPS-RL motion controller and Polytec OFV-500 vibrometer controller control the scanning stage and vibrometer, respectively. Two amplifiers are utilized—one for transmission (Techron 5530; AE Techron, Elkhart, IN) and the other for reception (Olympus 5077PR; Olympus, Shinjuku, Japan). Finally, a National Instruments data acquisition system (chassis, PXIe-1071; National Instruments, Austin, TX) with arbitrary waveform generator (PXIe-5423) and oscilloscope (PXIe-5172) cards, controlled via python commands, simultaneously manages excitation, reception, and scanning motion.

Laser Doppler vibrometry (LDV) scanning is performed in steps of 0.5 cm along the direction of wave propagation over a distance of 1 m (201 points), starting from 18 cm from the source. At each scanning point, the out-of-plane particle velocity is recorded for 40 ms at a 100 kHz sampling frequency and averaged 32 times to enhance the signal-to-noise ratio. The normalized waterfall plot without any metabarrier (baseline), after applying a bandpass filter (100–1000 Hz), demonstrates the propagation of PSV surface modes [Fig. 4(b)]. Fast-moving quasi-longitudinal P-waves, along with low-amplitude scattered sidewall reflections and strong backwall reflections near the end of the scanning region, are also observable in the waterfall plot. To mitigate the influence of sidewall reflections and quasi-P-waves in visualizing surface modes, a wavenumber-frequency dispersion spectrum is constructed by performing a two-dimensional (2D) fast Fourier transform (FFT) over the waterfall data between two border lines (marked in white dashed lines) where surface waves are expected to propagate [Fig. 4(c)]. The dispersion spectrum reveals the propagation of PSV₁ and PSV₂ surface modes, with the numerical dispersion curves derived from the measured properties of glass beads (γ_L = 21, γ_S = 8.8, α_L = 0.3, and α_S = 0.33) through curve-fitting plotted for reference.

Two metabarrier configurations comprising a 2D array of 3D-printed rod-like resonators made of polylactic acid (PLA) material are investigated, as shown in Fig. 4(a). Configuration #1 features a 4 × 6 arrangement with dimensions W = 2 cm, H_R = 6 cm, H_B = 2 cm, and A = 3.75 cm to illustrate the formation of the PSV₁ mode bandgap. Configuration #2 features a 20 × 12 arrangement with dimensions W = 0.5 cm, H_R = 8 cm, H_B = 5 cm, and A = 1.25 cm [Fig. 1(a)] to illustrate the formation of the PSV₂ mode bandgap. In both configurations, the resonators are buried 33 cm from the source, allowing the PSV₁ and PSV₂ modes to separate before interacting with the resonators, leaving transmission regions of about 55 cm and 45 cm for configurations #1 and #2, respectively. The resonator dimensions, buried depth, and spacing are carefully chosen through a parametric study to yield bandgaps within the 600–1000 Hz frequency range, as highlighted in the corresponding hybridized dispersion curves [Figs. 5(a) and 5(d)]. The analysis frequency range in experiments is limited by the transmission region length (less than two wavelengths of the PSV₂ mode below 500 Hz) and the magnetic shaker's inability to excite above 1000 Hz. For configuration #1, with the cutoff frequency of the hybridized longitudinal resonance mode at 750 Hz, the PSV₁ mode bandgap is expected to initiate above 770 Hz [Fig. 5(a)]. For configuration #2, the absence of the hybridized flexural resonance mode between 710 and 870 Hz is expected to initiate the PSV₂ mode bandgap [Fig. 5(d)].

The waterfall plots obtained in the presence of metabarrier configurations #1 and #2 are shown in Figs. 5(b) and 5(e), with the metabarrier region and transmission region (TR) marked for reference. Surface mode attenuation is evident in both configurations in the transmission region. Configuration #1, which targets the PSV₁ mode bandgap, exhibits more pronounced wave reflection, consistent with our frequency-domain analysis [Fig. 3(b)]. The normalized dispersion spectra for configurations #1 and #2, computed from the waterfall data over the transmission region, are shown in Figs. 5(c) and 5(f), respectively. The dispersion spectrum of configuration #1 clearly shows PSV₁ mode attenuation above 750 Hz, while configuration #2 shows a notable PSV₂ mode attenuation from 650 to 800 Hz [marked in Figs. 5(c) and 5(f)].

Finally, the transmission loss for PSV₁ and PSV₂ modes within the 600–1000 Hz frequency range in the presence of a metabarrier is shown in Figs. 5(g) and 5(h) for configurations #1 and #2, respectively. The transmission loss is assessed by normalizing the spectral amplitudes (averaged over two experiments) obtained along the PSV₁ and PSV₂ dispersion curves in the dispersion spectrum evaluated over the transmission region when a metabarrier is present, with corresponding spectral amplitudes (averaged over three experiments) extracted without a metabarrier. For effective normalization, the time series and spectral content across the 600–1000 Hz frequency range are ensured to be similar at the initial scan point (before the metabarrier region) for all experiments. To ensure reliability, the steel tank is gently shaken before each experiment, and the magnetic shaker is remounted: If shaking causes consolidation of the glass beads, the tank is refilled using a mesh sieve to restore the material properties. To account for minor discrepancies in the numerical and experimental dispersion curves, which may arise from variations between experiments due to the consolidation of glass beads as the tank is shaken between trials, we consider a bound of ±7.5 rad/m for each frequency around the numerical dispersion curves of PSV₁ and PSV₂ to compute the maximum spectral amplitude between the bounds and obtain the transmission spectra for the PSV₁ and PSV₂ modes. The shaded regions in the transmission plots depict fluctuations in the estimated PSV₁ and PSV₂ spectral amplitudes across various baseline experiments, underscoring the measurement sensitivity, as observed similarly in prior studies (Mora , 2021). Additionally, to corroborate the experimental findings, the transmission loss derived from full-scale 3D time-domain simulations performed over a finite element model replicating the experimental setup for both metabarrier configurations is plotted in Figs. 5(g) and 5(h). Details concerning the transient analysis are presented in the supplementary material.

The experimentally estimated transmission loss for both PSV₁ and PSV₂ modes closely aligns with those obtained from transient finite element simulations for both configurations, thereby demonstrating the reliability of experiments in identifying the PSV₁ and PSV₂ bandgaps [Figs. 5(g) and 5(h)]. A notable drop in the PSV₁ mode transmission (−15 dB) occurs after 770 Hz in both experiments and the time-domain simulation for configuration #1, effectively validating the presence of the PSV₁ mode bandgap; meanwhile, the PSV₂ mode transmission oscillates around a 3 dB loss [Fig. 5(g)]. Conversely, a notable transmission dip (∼10 dB) is observed in the PSV₂ mode between 650–875 Hz for configuration #2, slightly less pronounced than observed in simulations (15 dB), effectively illustrating the presence of the PSV₂ mode bandgap [Fig. 5(h)]. The slight compromise in the PSV₂ mode transmission dip may be attributed to the slightly varying buried depths across the metabarrier, expected due to manual resonator placement. Additionally, the obtained PSV₂ mode transmission dip, computed both numerically and experimentally, is slightly wider than the expected dip from the dispersion analysis (710–870 Hz), as the wave speed decreases with the gradual decrease in the slope of the flexural resonance mode in the dispersion curves closer to the cutoff frequencies [Fig. 5(d)]. The experimental and numerical PSV₁ mode transmission loss for configuration #2 oscillates around –4 dB throughout the frequency range without exhibiting any notable wide bandgap, except for a very narrow dip observed around 900 Hz in experiments. The presented experimental investigations confirm the formation of PSV₁ and PSV₂ mode bandgaps, which arise due to their preferential hybridization with longitudinal and flexural resonances, respectively, as demonstrated earlier through frequency-domain analysis (Fig. 3). However, it is worth noting that these bandgaps are distinct, requiring a comprehensive understanding of the influence of geometric parameters of the metabarrier to motivate a rational design aimed at suppressing both PSV₁ and PSV₂ modes, as discussed in the Sec. IV.

The influence of geometric parameters on two critical factors governing the PSV₁ and PSV₂ bandgaps is investigated: (i) the longitudinal-resonance mode cutoff frequency, determining the PSV₁ mode bandgap; and (ii) the absence of the flexural resonance mode in the hybridized dispersion curves (referred to as the “flexural resonance bandgap”), governing the PSV₂ mode bandgap. Inspired by the analysis of partially buried prismatic resonators for Rayleigh wave suppression in the homogeneous media (Liu , 2020), we define a new parameter β as the ratio of the buried depth (H_B) to the resonator height (H_R). To understand the combined influence of the geometric parameters on the longitudinal-resonance cutoff frequency, we vary the parameters β from 0.25 to 0.75, resonator height (H_R) from 4 m to 10 m, and resonator width (W) from 1 m to 4 m while maintaining the filling fractions (⁠ $FF=W2/A2$⁠) at 0.11 [Fig. 6(a)], 0.25 [Fig. 6(b)], and 0.44 [Fig. 6(c)]. The lower limit for β (0.25) ensures the resonator is buried at least one-fourth of its length for stability, while the upper limit (0.75) provides a substantial range for parametric analysis. This analysis reveals β as a useful tuning parameter for the longitudinal-resonance mode cutoff frequency; lower β values correspond to lower longitudinal-resonance cutoff frequencies, enabling effective suppression of the PSV₁ mode at lower frequencies.

It is crucial to note that while the initiation of the PSV₁ mode bandgap is ensured above the longitudinal-resonance mode cutoff frequency [as in Fig. 3(a)], the full extent of the PSV₁ mode bandgap is not accurately reflected in these results. This is because complex higher-order modes tend to initiate at higher frequencies, which can enable PSV₁ mode transmission, as documented in the supplementary material in the section on frequency-domain analysis to identify the upper limit of the PSV₁ mode bandgap. Additionally, the influence of β on the longitudinal resonance mode cutoff frequency diminishes with increasing resonator width, particularly noticeable for lower filling fractions. Further, the resonator height is demonstrated to have little influence on the variation of the cutoff frequency with β, especially for a larger resonator width. In conclusion, the analysis suggests that the cutoff frequency can be reduced to initiate PSV₁ bandgaps at frequencies as low as 7 Hz [Fig. 6(a)]. However, this comes at the expense of larger resonator width (W) and spacing (lower filling fraction), as well as reduced buried depth (lower β) that enable weaker constraints on the cantilever-like beam imposed by the granular media (Liu , 2019, 2020).

Moving forward with identifying the geometric parameters contributing to the PSV₂ mode bandgap, we observe that increased resonator width (not shown here) counteracts the formation of the flexural resonance bandgap necessary to suppress the PSV₂ mode. Therefore, we restrict the analysis to metabarriers with a resonator width (W) of 1 m and study the influence of the buried depth (H_B), resonator height (H_R), and filling fraction (FF) on the flexural resonance bandgaps. We examine buried depths (H_B) of 2, 4, and 6 m, adjusting the resonator heights (H_R) for each depth to vary β (⁠ $HB/HR$⁠) from 0.75 to 0.25, and present results for the filling fractions of 0.25 and 0.44, similar to the previous analysis. No evident bandgaps are observed for the filling fraction of 0.11 and are therefore ignored. For reference, the lower (β = 0.25) and upper (β = 0.75) β limits are represented as vertical lines in Figs. 7(a) and 7(b) for each buried depth.

We observe that for a lower filling fraction (FF = 0.25) and buried depth (H_B = 2 m), flexural resonance bandgaps are absent for smaller resonator heights [Fig. 7(a)]. With increasing buried depth, though these bandgaps initiate for all the considered resonator heights, they occur at higher frequencies. Obtaining the flexural resonance bandgap at lower frequencies requires longer resonators, as evident in Figs. 7(a) and 7(b), making it challenging to suppress the PSV₂ mode. However, a higher filling fraction of the metabarrier (FF = 0.44), though it increases the design complexity, is observed to initiate bandgaps even for shorter resonators buried at shallow depths and generate wider bandgaps compared to metabarriers with a lower filling fraction (FF = 0.25). While the intrinsic mechanism behind generating flexural resonance bandgaps for PSV₂ modes is complex, the presented parametric analysis indicates that wider flexural resonance bandgaps at lower frequencies necessitate smaller resonator width, longer resonators, and higher filling fraction. This requirement contrasts with the need for a large resonator width in metabarriers to achieve lower longitudinal-resonance mode cutoff frequencies necessary to suppress the PSV₁ mode at lower frequencies [Fig. 6(a)]. Consequently, designing metabarriers capable of simultaneously suppressing both PSV₁ and PSV₂ modes over a wide frequency range becomes challenging and requires new innovative approaches, such as exploiting graded metabarriers as studied in the subsequent section.

Despite their demonstrated potential to enhance bandgaps for Rayleigh waves in homogeneous media (Colombi , 2016), graded metabarriers (GMs) with varying resonance frequencies have not yet been explored for suppressing PSV modes in granular media. The steady decrease in the flexural resonance bandgap frequencies with increasing resonator height (Fig. 7) motivates the application of GMs with varying resonator heights to widen the PSV₂ mode bandgap. While it is not immediately intuitive to use GMs to suppress the PSV₁ mode, this section aims to demonstrate the ability of GMs to induce mode conversions between the PSV₁ and PSV₂ modes that can be leveraged to simultaneously suppress both modes. To illustrate this application, consider two GMs (GM #1 and GM #2), each consisting of 40 partially embedded rod-like resonators. The other tunable parameters, except the resonator height, are kept constant in these two metabarriers: W = 1 m, A = 1.75 m, and H_B = 4 m, as previously used for frequency-domain analysis in Fig. 3. In GM #1, the resonator heights (H_R) range from 7 m to 11 m, utilizing the flexural resonance bandgap between the first (F₁) and second (F₂) flexural resonance modes, which spans the desired 8–40 Hz frequency range, as illustrated in Figs. 8(a) and 8(b). In contrast, resonator heights (H_R) are considered to vary from 12 m to 20 m in GM #2 to utilize the flexural resonance bandgap between the second (F₂) and third (F₃) flexural resonance modes, also spanning the desired 8–40 Hz frequency range, as illustrated in Figs. 8(b) and 8(c).

A notable contrast between these two GMs is the existence of a mode coupling frequency in GM #2 within the desired frequency range (8–40 Hz), where the first longitudinal resonance mode (L) intersects with the second flexural resonance mode (F), evident across the hybridized dispersion curves for all the resonator heights considered in the metabarrier [Fig. 8(c)]. Due to the preferential hybridization of the PSV₁ and PSV₂ modes with the longitudinal and flexural resonances [Figs. 3(a) and 3(d)], respectively, a significant mode conversion between these modes is anticipated at this mode coupling frequency. In GM #2, the mode coupling frequency shifts up or down within the 10–30 Hz range, depending on whether the resonator heights increase or decrease in the direction of the incident wave [see the hybridized dispersion curves for the shortest and longest resonator in Fig. 8(c)]. This broad frequency range for mode conversion can be exploited to suppress both PSV₁ and PSV₂ modes simultaneously, as will be illustrated in this section. To that front, we individually explore the interaction of PSV₁ and PSV₂ modes with the GMs #1 and #2, considering both the increasing and decreasing resonator heights, to highlight the nuances in the mode conversions between the PSV₁ and PSV₂ modes.

The interaction between the incident PSV₁ and PSV₂ modes and GMs #1 and #2, featuring increasing resonator height, is depicted in Fig. 9 through frequency-domain finite element analysis, analogous to the analysis performed earlier in Sec. II. In addition to the transmission spectra, we include the reflection spectra in the analysis to comprehensively evaluate the metabarrier efficiency against the incident PSV₁ or PSV₂ modes. In GM #1, where no mode coupling frequency exists, the transmission spectra for the incident PSV₁ mode, owing to its affinity towards longitudinal-resonance mode hybridization, reveal significant PSV₁ mode transmission across the frequency range until reaching the minimum longitudinal-resonance mode cutoff frequency $fL(min)$ = ∼34 Hz [marked in Figs. 9(a) and 9(c)], attributed to the longest resonator in the metabarrier. This substantial PSV₁ mode transmission is further evidenced by the out-of-plane displacement field depicted at one of the excitation frequencies, as shown in the inset of Fig. 9(a). Although the mode conversions to transmitted (and reflected) PSV₂ modes remain relatively low, consistent with the uniform resonator height metabarrier [Fig. 3(a)], we note significant PSV₁ mode reflection across the frequency range, particularly for frequencies closer to and above $fL(min)$ [Fig. 9(a)].

On the other hand, as expected with the flexural resonance bandgap spanning the 8–40 Hz frequency range [Fig. 8(a)], we observe the PSV₂ mode suppression across all excitation frequencies, emphasizing its superiority over a uniform resonator height metabarrier [compare Fig. 9(c) and Fig. 3(b)]. The displacement field in the inset of Fig. 9(c) illustrates a notable localization of the incident PSV₂ mode into the granular half-space beneath the resonators, indicating a similar mechanism of PSV₂ mode suppression observed with the uniform resonator height metabarrier [Fig. 3(d)]. Furthermore, we note a relatively weak mode conversion (reflected and transmitted) to the PSV₁ mode across the frequency range in the transmission/reflection spectra, accompanied by a substantial reflection as a PSV₂ mode at higher frequencies. In conclusion, while this metabarrier configuration is ineffective for PSV₁ mode suppression, it is highly effective in attenuating the PSV₂ mode across a wide frequency range.

GM #2, exhibiting a mode coupling frequency across the 10–30 Hz frequency range, demonstrates unique interactions with the incident PSV₁ and PSV₂ modes, as depicted in Figs. 9(b) and 9(d). Although the longitudinal-resonance mode persists throughout the frequency range until $fL(min)$ = 25.5 Hz [marked in Figs. 9(b) and 9(d)], the incident PSV₁ mode is predominantly reflected upon encountering the resonator whose mode coupling frequency in its hybridized dispersion curves matches the excitation frequency. This behavior is illustrated in the transmission/reflection spectra in Fig. 9(b) and the inset showing the out-of-plane displacement field. We also observe an anomaly in the transmission/reflection spectra around 20 Hz with considerable transmission. While the cause of this phenomenon is unknown, we speculate that the efficacy of the mode coupling in reflecting the incident PSV₁ mode diminishes as frequencies approach the $fL(min)$ frequency. Additionally, mode conversions to transmitted and reflected PSV₂ modes are weak across the frequency range.

Remarkably, for an incident PSV₂ mode, we observe a broadband suppression of the PSV₂ mode due to its conversion to the PSV₁ mode for frequencies below $fL(min)$⁠, as illustrated in Fig. 9(d) and the inset showing the simulation snapshot. This “downconversion” in the surface mode order (PSV₂ to PSV₁) is attributed to the mode coupling exhibited by the resonators within the graded metabarrier. After surpassing the $fL(min)$⁠, most of the incident PSV₂ mode is observed to be significantly reflected back as the PSV₂ mode with a weak mode-converted PSV₁ mode, as illustrated in the transmission/reflection spectra. In conclusion, while the GM #2 with increasing resonator heights effectively suppresses the PSV₂ mode, albeit through a distinct mechanism from GM #1, and efficiently reflects the incident PSV₁ mode, the mode conversion of the incident PSV₂ mode to PSV₁ compromises the metabarrier's efficacy in suppressing the PSV₁ mode.

The interaction of the incident PSV₁ and PSV₂ modes with the GMs #1 and #2, featuring decreasing resonator heights, is illustrated in Fig. 10. As expected, the performance of GM #1, with decreasing resonator heights, towards incident PSV₁ and PSV₂ modes closely resembles that of increasing resonator heights [compare Fig. 10(a) and Fig. 9(a) and Fig. 10(c) and Fig. 9(c)], owing to the absence of mode coupling frequency among resonators within the frequency range of interest [Fig. 8(a)]. Therefore, the conclusions drawn for the metabarrier with increasing resonator heights, namely, their effectiveness in predominantly suppressing the PSV₂ mode across a broad frequency range, apply to the case of decreasing resonator heights. Hence, in the absence of mode coupling, the behavior of graded metabarriers towards surface modes in granular media remains unaffected by the choice of grading (increasing or decreasing), which contrasts with their behavior towards surface waves in a homogeneous medium (Colombi , 2016).

In GM #2 with decreasing resonator heights, PSV₁ mode suppression is attained through a notable mode conversion (“upconversion” in the surface mode order) of the incident PSV₁ mode to PSV₂ [Fig. 10(b)]. This mechanism of PSV₁ mode suppression contrasts with that exhibited by the same metabarrier configuration but with increasing resonator heights, where the incident PSV₁ mode is significantly reflected without any mode conversion [compare Fig. 10(b) and Fig. 9(b)]. However, above the $fL(min)$ [marked in Figs. 10(b) and 10(d)], the metabarriers with increasing and decreasing height behave similarly, with the incident PSV₁ mode predominantly reflected. We also observe a similar anomaly with enhanced PSV₁ mode transmissions around 20 Hz, as seen in the case of increasing resonator height [Fig. 9(b)]. On the other hand, unlike the mechanism where PSV₂ mode suppression is achieved by mode converting to the PSV₁ mode for GM #2 with increasing resonator heights, in the case of decreasing resonator heights, PSV₂ mode suppression predominantly occurs by localizing the PSV₂ mode into the half-space without significant mode conversions, similar to the GM #1 [compare Fig. 10(d) and Fig. 9(d)]. At higher frequencies beyond the $fL(min)$⁠, most of the incident PSV₂ mode is observed to be reflected, similar to that demonstrated for GM #2 with increasing resonator heights.

In conclusion, while GM #2 with decreasing resonator heights effectively suppresses the PSV₁ mode and efficiently localizes the incident PSV₂ mode within the half-space, the mode conversion of the incident PSV₁ mode to PSV₂ compromises the metabarrier's efficacy in suppressing the PSV₂ mode, which is strikingly opposite to the behavior observed with increasing resonator heights. Overall, these investigations on graded metabarriers suggest that while these structures can generate wider bandgaps for either PSV₁ or PSV₂ modes individually, achieving broadband suppression of both modes simultaneously using graded metabarriers alone remains challenging. However, the distinctive mode-conversions demonstrated by these graded metabarriers offer a pathway to simultaneously suppress the PSV₁ and PSV₂ modes through compound metabarriers, which involve appending two distinct metabarriers, as discussed in the subsequent section.

Suppressing the PSV₁ mode poses a greater challenge compared to the PSV₂ mode because graded metabarriers can achieve wider PSV₂ mode bandgaps using shorter resonators (such as the GM #1) without relying on mode coupling frequencies [see Figs. 9(a) and 9(c) and Figs. 10(a), and 10(c)]. One way to design a compound metabarrier capable of simultaneously suppressing both PSV₁ and PSV₂ modes is to combine this graded metabarrier with PSV₂ mode suppression capability with another metabarrier capable of suppressing the PSV₁ mode. There are two feasible options in metabarriers for suppressing PSV₁ modes at low frequencies. First, the uniform resonator height metabarriers achieve this by tuning the longitudinal resonance mode cutoff frequencies $(fL)$ to very low values. However, this approach requires larger resonator widths (W > 2 m) to obtain a cutoff frequency below 15 Hz [Figs. 6(a) and 6(b)]. Second, the GMs with decreasing resonator heights and mode coupling frequencies (GM#2) can suppress the PSV₁ mode by mode converting to PSV₂ mode, which requires thinner but longer resonators [Fig. 10(b)]. In this paper, we selected the latter option as a case study for a compound metabarrier. It merges GM #2, featuring decreasing resonator heights, as previously discussed for its efficacy in suppressing the PSV₁ mode by mode converting to the PSV₂ mode [Fig. 10(b)], with an appended GM #1 featuring increase resonator heights, which effectively suppresses both the incident PSV₂ mode and the mode-converted PSV₂ mode transmitted from the former metabarrier [Fig. 9(c)].

The performance of this compound metabarrier over the incident PSV₁ and PSV₂ modes is illustrated through the transmission/reflection spectra in Figs. 11(a) and 11(b), respectively. For the incident PSV₁ mode, as depicted in the displacement field inset [Fig. 11(a)], the majority of the energy undergoes mode conversion to the PSV₂ mode at the first graded metabarrier (GM #2), which is then directed into the half-space by the subsequent metabarrier (GM #1). A notable suppression (>75%) of the incident PSV₁ mode is evident in the transmission spectra across the frequency range of 8–40 Hz, except for frequencies nearer to the $fL(min)$⁠, where the mode conversion from the incident PSV₁ mode to the PSV₂ mode is compromised, as previously highlighted for GM #2 in Fig. 10(b). Most of the incident PSV₁ mode energy after the $fL(min)$ is significantly reflected following the observations from Fig. 10(b).

Similarly, we observe significant suppression of the incident PSV₂ mode (> 50%) for the incident PSV₂ mode throughout the frequency range (10–40 Hz). A substantial reflection as a PSV₂ mode at higher frequencies is observed, similar to that observed for GM #1 in Fig. 9(c). The residual transmission of the PSV₂ mode observed at lower frequencies for compound metabarrier [Fig. 11(b)] is attributed to the slight recovery of PSV₂ mode energy localized into the half-space, which resurfaces due to the mirage effect (Palermo , 2018), as similarly observed for GM #1 in Fig. 9(c). By increasing the number of resonators or employing multiple, spaced, compound metabarriers, we can further enhance the suppression of both PSV₁ and PSV₂ modes over a wider frequency range. This is because the mode-converted PSV₂ mode in the case of the incident PSV₁ mode and the transmitted PSV₂ mode in the case of incident PSV₂ mode can be gradually localized into the half-space by employing multiple metabarriers.

A comprehensive understanding of the interactions of the surface modes with a partially buried metabarrier in granular media is crucial for the rational design of metabarriers. In contrast to the hybridization of the non-dispersive Rayleigh surface waves with both the longitudinal and flexural resonances of the surface resonators in an isotropic and homogeneous half-space, which result in the formation of surface wave bandgaps, the hybridization behavior of the PSV dispersive surface modes in heterogeneous granular media with the local resonances (longitudinal/flexural) is complex and requires extensive numerical analysis. In this study, the efficacy of a metabarrier comprising partially buried rod-like resonators in suppressing the relevant surface modes (PSV₁ and PSV₂) in granular media is numerically and experimentally evaluated. The dispersion analysis of the unit cell comprising a partially buried resonator and granular half-space reveals multiple hybridizations between the PSV surface modes and the longitudinal and flexural resonances of the resonators with no discernible bandgaps contrary to the well-studied hybridization of Rayleigh waves with the longitudinal and flexural resonances of the resonator.

The frequency-domain analyses employing body load excitation reveal preferential hybridization of PSV₁ and PSV₂ modes with longitudinal and flexural resonances, respectively. The preferential hybridization of the PSV₁ mode with the longitudinal resonance facilitates PSV₁ mode suppression after the longitudinal-resonance cutoff frequency. In contrast, the preferential hybridization of the PSV₂ mode with the flexural resonances enable PSV₂ mode suppression at frequency ranges where no possible hybridized flexural resonance modes exist (flexural resonance bandgaps). Unlike the incident PSV₁ mode, which predominantly reflects from the metabarrier above the longitudinal-resonance cutoff frequency, the incident PSV₂ mode is observed to localize and propagate under the metabarrier with significant suppression at these flexural resonance bandgaps. The presence of PSV₁ and PSV₂ mode bandgaps is experimentally validated in a granular media testbed comprising glass beads at scaled frequencies (600–1000 Hz) using 3D-printed PLA rod-like resonators. We present an extensive parametric analysis to understand the influence of unit-cell geometry on the positioning of the PSV₁ and PSV₂ mode bandgaps. Exploiting the dependency of the flexural bandgap positioning on resonator length and the intriguing mode coupling behaviors, we demonstrate an innovative application of graded metabarriers as broadband filters for PSV₂ mode and as mode-conversion filters that facilitate conversion between PSV₁ and PSV₂ modes. Finally, using the mode-conversion capability of the graded metabarriers, we propose to simultaneously suppress both PSV₁ and PSV₂ modes at low frequencies (8–40 Hz) using a compound metabarrier configuration, a combination of two distinct graded metabarriers. It is apparent that these initial metabarrier designs have quite large dimensions; Nonetheless, they offer valuable insights into the hybridization of the surface modes in granular media with longitudinal/flexural resonances and establish a foundation for more optimized designs. The natural next steps include considerations of soil's non-elastic behavior, stratification, etc., for more realistic metabarrier designs.

See the supplementary material for the 3D time-domain finite element analysis demonstrating the challenges in suppressing surface waves in granular media using previously proposed spring-mass-type metabarriers and the proposed rod-like metabarriers. The supplementary material also includes a 3D full-scale time-domain analysis mimicking the experimental setup to demonstrate PSV1 and PSV2 mode bandgaps in partially embedded rod-like resonators, as well as a frequency-domain finite element analysis performed to identify the upper limit of the PSV1 mode bandgap.

The authors gratefully acknowledge the support of the National Science Foundation under Grant No. 1934527. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. Computations for this research were performed on the Pennsylvania State University's Institute for Computational and Data Sciences' Roar supercomputer. The authors gratefully thank Luke B. Beardslee from Los Alamos National Laboratory and Jeremy Keirn and Prabhav Borate from Penn State University for their support during our development of the experimental setup. The authors thank the reviewers of this manuscript for providing valuable insights, which helped improve the quality of the manuscript.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.