# Seismic wave interaction with buried cavity networks: Analytical modeling and resonance effects

Agnès Maurel, Stéphane Brulé, Sébastien Guenneau, Kim Pham

► To cite this version:

Agnès Maurel, Stéphane Brulé, Sébastien Guenneau, Kim Pham. Seismic wave interaction with buried cavity networks: Analytical modeling and resonance effects. Wave Motion, 2026, 141, pp.103666. 10.1016/j.wavemoti.2025.103666. hal-05361262

HAL Id: hal-05361262

https://cnrs.hal.science/hal-05361262v1

Submitted on 12 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ELSEVIER

<figure><img src="image_3.png" /></figure>

# Seismic wave interaction with buried cavity networks: Analytical modeling and resonance effects

Agnès Maurel<sup>a, b, ∗</sup>, Stéphane Brulé<sup>b</sup>, Sébastien Guenneau<sup>c, d</sup>, Kim Pham<sup>d</sup>

## ARTICLE INFO

### Keywords:

- Seismic waves

- Buried cavities

- Homogenization

- Resonances

## ABSTRACT

We study the scattering of elastic waves by a periodic array of cavities buried in an elastic half-space. This configuration is relevant in seismology, where shallow voids can locally amplify ground motion. Building on homogenized interface models developed for infinite media, we extend the approach to account for the presence of a stress-free surface. The resulting model yields an analytical solution to the 2D elastodynamic problem for incident longitudinal L and transverse T waves. A semi-analytical multimodal solution is used for validation. The analysis reveals the conditions under which resonances occur in the soil layer between the cavity tops and the surface, with particular emphasis on the low-frequency resonance that dominates in seismic contexts. The model identifies the key parameters governing resonance and provides insights into the transition from infinite to finite cavity arrays. It offers a simplified yet accurate framework for assessing site-specific seismic amplification.

## 1. Introduction

The study of imperfect interfaces in elastic media has long attracted attention, both for its theoretical interest and its practical applications in nondestructive evaluation, seismology, and layered composite analysis [1–4]. Among the most emblematic configurations, periodic arrays of thin cracks or voids embedded in a homogeneous elastic solid have served as canonical models to develop effective models in terms of transmission conditions capable of capturing complex scattering effects in a simplified form. A major step in this direction was achieved by Achenbach and co-workers, who investigated the scattering properties in many situations including equally spaced collinear cracks [5,6], inclined cracks [7] and spherical cavities [8]; a review is presented in [9]. Beyond their practical use in engineering, such models have proven essential for understanding the interaction between elastic waves and periodic microstructures.

When such periodic structures are embedded in a half-space rather than in an unbounded medium, the presence of the free surface introduces new and significant physical phenomena. In particular, buried defects can interact with incident surface or body waves in ways that generate complex scattering patterns, amplification effects, and resonances. These effects are especially relevant in seismology, where the elastic wavefield is strongly influenced by the subsurface geometry. This configuration is of both theoretical and practical importance, as it can lead to localized energy trapping and potentially strong surface amplifications. Such behavior was documented in a detailed historical study of the 1930 earthquake in southern Italy, specifically in the town of Rionero in Vulture,

<footnote><sup>a</sup> Institut Langevin, CNRS, ESPCI-ParisTech, 1 rue Jussieu, Paris, 75005, France</footnote>

<footnote><sup>b</sup> Menard Company, 1 route du Dôme, Chaponost, 69630, France</footnote>

<footnote><sup>c</sup> UMI 2004 Abraham de Moivre - CNRS, Imperial College, SW7 2AZ, London, UK</footnote>

<footnote><sup>d</sup> LMI, ENSTA Paris, 828 bd des maréchaux, Palaiseau, 91762, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: agnes.maurel@espci.fr (A. Maurel), stephane.brule@menard-mail.com (S. Brulé), s.guenneau@imperial.ac.uk (S. Guenneau), kim.pham@ensta.fr (K. Pham).</footnote>
where highly nonuniform seismic damage was observed — damage that could not be explained by building typology alone [10,11]. A systematic correlation was found between the most severely affected areas and the presence of dense networks of man-made underground cavities, with resonance effects suspected when the soil cover above the voids was only a few meters thick. These observations highlight the need for simple yet accurate predictive models capable of capturing the key parameters that control such resonances, particularly in the low-frequency range relevant to urban seismic risk. However, most studies to date consider only the effect of a single buried cavity [12–16]. A noticeable exception is the study by [17], which compares the effects of one and two cavities, and shows that amplification nearly doubles when moving from a single cavity to a pair. This suggests that the single-cavity case is unlikely to be representative of a cavity network, and that collective effects in such periodic systems must be addressed specifically. These collective effects, arising from the periodic arrangement of subsurface heterogeneities, are now well recognized in the broader field of wave physics. In particular, the past decade has seen the emergence of seismic metamaterials, where artificially structured soils or periodic arrays of boreholes are designed to control the propagation of seismic waves through Bragg or local resonance mechanisms [18–21].

The present work aims to provide a physically grounded model for such buried cavity systems by adapting the homogenized imperfect interface approach developed in [22] for infinite media to the case of a periodic array of cavities buried in an elastic half-space. We derive an analytical solution to the 2D elastodynamic scattering problem for both longitudinal and transverse incident waves, using a generalized transmission condition to represent the cavity layer. The proposed model allows us to identify the key geometric and mechanical parameters controlling the resonance frequencies and amplitudes, including not only the burial depth but also the density and cavity shape of the array. To assess the validity of this homogenized approach, we construct a reference numerical solution based on a semi-analytical multimodal formulation, detailed in Section 2 alongside the formulation of the effective model. Section 3 presents a systematic comparison between numerics and effective model and offers a broader discussion of the resonance mechanisms at play. In particular, we investigate the extent to which results obtained for an infinite periodic system can inform the understanding of more realistic configurations involving a finite number of cavities. Finally, Section 4 focuses on the low-frequency resonance mode – of primary importance for seismic applications – providing explicit formulas and scaling laws for its frequency and amplitude, and illustrating the potential of the model for practical site response assessment in the presence of shallow subsurface heterogeneities.

## 2. Configuration of interest

The configuration of interest is shown in Fig. 1. We consider a two-dimensional homogeneous elastic medium characterized by its Lamé coefficients $(\lambda, \mu)$ and mass density $\rho$, occupying the half-space $z \in (-\infty, 0)$, except in a region $z \in (-d-a, -d)$, which contains an array of cavities, or voids. This array is perfectly periodic with period $h$, and the cavities have characteristic dimensions $e$ and $a$. In the following, we denote $d_a = d + a$, and define the array density as $\phi = e/h$. The study focuses on the scattering of incident longitudinal (L) and transverse (T) waves.

### 2.1. Governing equations for the direct and effective problems

We consider the two-dimensional plane strain equations of linear elastodynamics in the $(x, z)$ plane, within the lower half-space defined by $x \in (-\infty, +\infty)$ and $z \in (-\infty, 0)$. The governing equations for the displacement vector field $\mathbf{u} = (u_x, u_z)^T$ are given by:

$$ \mathrm{div}\sigma + \rho\omega^2\mathbf{u} = 0, \quad \sigma = 2\mu\epsilon + \lambda\mathrm{tr}(\epsilon)\mathbf{I}, \quad \epsilon = \frac{1}{2}(\nabla\mathbf{u} + (\nabla\mathbf{u})^T), \quad (1) $$

where $\sigma$ is the stress tensor, $\epsilon$ is the strain tensor, $\omega$ denotes the angular frequency, $\rho$ is the mass density, and $(\lambda, \mu)$ are the Lamé parameters. The identity tensor is denoted by $\mathbf{I}$, and $\mathrm{tr}(\epsilon)$ denotes the trace of the strain tensor. The components of the stress tensor $\sigma$ and the displacement vector $\mathbf{u}$ are written explicitly as:

$$ \sigma = \begin{pmatrix} \sigma_{xx} & \sigma_{xz} \\ \sigma_{xz} & \sigma_{zz} \end{pmatrix}, \quad \mathbf{u} = \begin{pmatrix} u_x \\ u_z \end{pmatrix}. $$

These equations are subject to stress-free boundary conditions at the surface $z = 0$, namely:

$$ \sigma_{xz}(x, 0) = \sigma_{zz}(x, 0) = 0. \quad (2) $$

#### 2.1.1. The direct problem

In the absence of cavities – that is, in a homogeneous elastic half-space – the problem defined by Eqs. (1)–(2) admits an explicit solution. This is the case, for instance, when the source consists of an incident plane L- or T-wave; see, e.g., Chapter 5 in [23]. This solution is referred to as the free-field. When cavities are present, however, the direct problem must be solved in the region $z \in (-d_a, -d)$, which contains the cavities. In this case, the boundary conditions read:

$$ \sigma \cdot n = 0, \quad \text{on the boundaries of the cavities,} \quad (3) $$

where $n$ denotes the unit normal vector to the cavity boundaries. The direct problem defined by (1)–(3) does not admit closed form solution and must therefore be solved numerically. In Section 2.3, we present the multimodal method used to compute the reference, numerical, solutions.
<figure><img src="image_4.png" /><figcaption>Fig. 1. Configuration of interest: Scattering of incident longitudinal ($A_L^{\text{inc}}$) and transverse ($A_T^{\text{inc}}$) waves by a periodic array of buried cavities in an elastic half-space in two dimensions. In the effective problem, the region of the cavity array is replaced by an imperfect interface with effective jump conditions (4).</figcaption></figure>

### 2.1.2. The effective problem

The effective problem consists in solving (1) and (2), as in the direct problem, but with the key difference that the cavity layer $z \in (-d_a, -d)$ is replaced by effective, imperfect interface. Indeed, as shown by [22], the effect of such an array can be rigorously encapsulated by a set of effective jump conditions, leading to the following homogenized interface model:

$$ 
\begin{align} 
[u_x] &= hB \left\langle \frac{\partial u_x}{\partial z} \right\rangle + (hB - a) \frac{\partial \langle u_z \rangle}{\partial x}, \tag{4} \\
[\sigma_{xz}] &= (1 - \varphi)a \left\langle \frac{\partial \sigma_{xz}}{\partial z} \right\rangle - \mu h \frac{\partial}{\partial x} \left( C_x \frac{\partial \langle u_x \rangle}{\partial x} + C_z \left\langle \frac{\partial u_z}{\partial z} \right\rangle \right), \quad [\sigma_{zz}] &= (1 - \varphi)a \left\langle \frac{\partial \sigma_{zz}}{\partial z} \right\rangle - \varphi a \frac{\partial \langle \sigma_{xz} \rangle}{\partial x}. 
\end{align} 
$$

where:

$$ [\nu] = \nu(-d) - \nu(-d_a), \quad \langle \nu \rangle = \frac{1}{2} (\nu(-d) + \nu(-d_a)), $$

denotes the jump of a field across the cavity layer and the average value and ($B, B_x, B_z, C_x, C_z$) are dimensionless effective coefficients obtained by solving elementary static problems (see Appendix A). These depend only on the geometry of the cavities and the material properties of the surrounding soil, and are independent of frequency. In the limit where the cavities vanish – i.e., when $a = e = 0$ – all effective coefficients reduce to zero, and the classical continuity conditions for displacement and normal stresses are recovered.

### 2.2. Potential formulation of the solution

In the remainder of this section, we describe the numerical resolution of the direct problem using potential theory. We introduce Helmholtz potentials ($\phi, \psi$) defined such that:

$$ (\lambda + 2\mu)\Delta\phi + \rho\omega^2\phi = 0, \quad \mu, \Delta\psi + \rho\omega^2\psi = 0. \tag{5} $$

The displacement components are reconstructed from these potentials as:

$$ u_x = \frac{\partial \phi}{\partial x} - \frac{\partial \psi}{\partial z}, \quad u_z = \frac{\partial \phi}{\partial z} + \frac{\partial \psi}{\partial x}. \tag{6} $$

The associated normal stress components $\sigma_{xz} = \mu (\partial_x u_z + \partial_z u_x)$ and $\sigma_{zz} = (\lambda + 2\mu)\partial_z u_z + \lambda \partial_x u_x$ are given by:

$$ \frac{\sigma_{xz}}{\mu} = 2 \frac{\partial^2 \phi}{\partial x \partial z} + \frac{\partial^2 \psi}{\partial x^2} - \frac{\partial^2 \psi}{\partial z^2}, \quad \frac{\sigma_{zz}}{\mu} = -k_T^2 \phi - 2 \frac{\partial^2 \phi}{\partial x^2} + 2 \frac{\partial^2 \psi}{\partial x \partial z}, \tag{7} $$

Without loss of generality, we consider incident plane waves of the form:

$$ \phi^{\text{inc}}(\mathbf{r}) = A_L^{\text{inc}} e^{i\alpha_L z} e^{i\beta x}, \quad \psi^{\text{inc}}(\mathbf{r}) = A_T^{\text{inc}} e^{i\alpha_T z} e^{i\beta x}, \tag{8} $$

where $A_L^{\text{inc}}$ and $A_T^{\text{inc}}$ are the amplitudes of the incident longitudinal (L) and transverse (T) waves, respectively (Fig. 1). The wavenumbers are defined by:

$$ k_L = \sqrt{\frac{\rho}{\lambda + 2\mu}} \omega, \quad k_T = \sqrt{\frac{\rho}{\mu}} \omega, \quad \alpha_L = k_L \cos \theta_L, \quad \alpha_T = k_T \cos \theta_T, \quad \beta = k_L \sin \theta_L = k_T \sin \theta_T. \tag{9} $$

with $\theta_L$ and $\theta_T$ being the incidence angles for the longitudinal and transverse waves.
## 2.3. Resolution of the direct problem: numerical method

To compute the reference solution to the full scattering problem, we employ a multimodal expansion technique directly inspired from [24]. This method leverages the periodicity of the cavity array to reduce the computation to a single elementary cell, repeated with horizontal periodicity $h$.

Let us define the computational domain as the union of three subdomains:

$$ \mathcal{Y} = \mathcal{Y}^{-} \cup \mathcal{Y}_{L} \cup \mathcal{Y}^{+} \qquad (10) $$

where:

$$ \mathcal{Y}^{-} = \{ x \in (-\frac{h}{2}, +\frac{h}{2}), z \in (-\infty, -d_a) \} \text{ corresponds to the semi-infinite soil layer beneath the cavities,} $$

$$ \mathcal{Y}^{+} = \{ x \in (-\frac{h}{2}, +\frac{h}{2}), z \in (-d, 0) \} \text{ corresponds to the soil layer above the cavities,} $$

$$ \mathcal{Y}_{L} = \{ x \in (-\frac{h-e}{2}, +\frac{h-e}{2}), z \in (-d_a, -d) \} \text{ corresponds to the vertical solid ligaments between two neighboring cavities, (see Fig. 2(a)).} $$

In the following, the potentials are expanded in suitable modal bases in each subdomain. In $\mathcal{Y}_L$, the basis consists of the Lamb modes of an elastic ligament of width $(h - e)$, which naturally satisfy the traction-free boundary conditions at $x = \pm(h - e)/2$. In the substrate $\mathcal{Y}^{-}$, the field representation given later in Eq. (11) involves the vertical wavenumbers associated with the longitudinal and transverse bulk modes, ensuring that the radiation condition as $z \to -\infty$ is exactly satisfied for all incidence angles, since only the downward-propagating (or exponentially decaying) components are retained.

**Modal expansions in the substrate** – In both $\mathcal{Y}^{-}$ and $\mathcal{Y}^{+}$, the solution is expanded as a superposition of incident, reflected, and transmitted modes. For a point $r = (x, z)$, we write:

$$ r \in \mathcal{Y}^{-} : \quad \begin{pmatrix} \phi(\mathbf{r}) \\ \psi(\mathbf{r}) \end{pmatrix} = \begin{pmatrix} \phi^{\text{inc}}(\mathbf{r}) \\ \psi^{\text{inc}}(\mathbf{r}) \end{pmatrix} + \sum_{n=-\infty}^{+\infty} \begin{pmatrix} R_{L,n} e^{-i\alpha_{L,n}z} \\ R_{T,n} e^{-i\alpha_{T,n}z} \end{pmatrix} F_n(x), \qquad (11) $$

$$ r \in \mathcal{Y}^{+} : \quad \begin{pmatrix} \phi(\mathbf{r}) \\ \psi(\mathbf{r}) \end{pmatrix} = \sum_{n=-\infty}^{+\infty} \begin{pmatrix} b_{L,n}^{+} e^{i\alpha_{L,n}z} + b_{L,n}^{-} e^{-i\alpha_{L,n}z} \\ b_{T,n}^{+} e^{i\alpha_{T,n}z} + b_{T,n}^{-} e^{-i\alpha_{T,n}z} \end{pmatrix} F_n(x), $$

where the modal functions in $x$ are defined by:

$$ F_n(x) = e^{i\beta_n x}. \qquad (12) $$

The vertical wavenumbers are given by the dispersion relations:

$$ \left\{ \begin{array}{ll} \alpha_{L,n} = \sqrt{k_L^2 - \beta_n^2}, & \alpha_{T,n} = \sqrt{k_T^2 - \beta_n^2}, & n = -\infty, \dots, +\infty, \\ \beta_n = \beta + \frac{2n\pi}{h}. & & \end{array} \right. \qquad (13) $$

The choice of square roots ensures upward or downward propagation depending on the sign of $\text{Im}(\alpha_{L,n})$ or $\text{Im}(\alpha_{T,n})$, and enforces radiation conditions. The periodic structure implies a pseudo-periodicity constraint on the solution $F_n(-h/2) = F_n(h/2)e^{i\beta_n h}$.

**Lamb modes in the solid ligaments** – In the ligament region $\mathcal{Y}_L$ separating two cavities, the solution is expanded in terms of symmetric and antisymmetric Lamb modes — i.e., the solutions of a traction-free elastic plate of finite width $(h - e)$. These modes inherently satisfy the stress-free conditions on the vertical boundaries of the cavity.

We write:

$$ r \in \mathcal{Y}_L : \quad \begin{pmatrix} \phi(\mathbf{r}) \\ \psi(\mathbf{r}) \end{pmatrix} = \sum_{n=0}^{+\infty} \left( \begin{pmatrix} c_{Sn}^{+} e^{i\alpha_{Sn}z} + c_{Sn}^{-} e^{-i\alpha_{Sn}z} \\ \psi_{Sn}^{+} e^{i\alpha_{Sn}z} + \psi_{Sn}^{-} e^{-i\alpha_{Sn}z} \end{pmatrix} \begin{pmatrix} \phi_{Sn}(x) \\ \psi_{Sn}(x) \end{pmatrix} + \sum_{n=0}^{+\infty} \left( c_{An}^{+} e^{i\alpha_{An}z} + c_{An}^{-} e^{-i\alpha_{An}z} \right) \begin{pmatrix} \phi_{An}(x) \\ \psi_{An}(x) \end{pmatrix} \right), \qquad (14) $$

where: $(\alpha_{Sn}, \alpha_{An})$ are vertical wavenumbers associated with the symmetric and antisymmetric modes, respectively, $(\phi_{Sn}(x), \psi_{Sn}(x))$ and $(\phi_{An}(x), \psi_{An}(x))$ are the corresponding horizontal modal functions, defined by:

$$ \left\{ \begin{array}{ll} \begin{pmatrix} \phi_{Sn}(x) \\ \psi_{Sn}(x) \end{pmatrix} = \begin{pmatrix} \cosh(\beta_{Sn}x) \\ B_{Sn} \sinh(\beta_{Sn}x) \end{pmatrix}, & B_{Sn} = -\frac{2i\alpha_{Sn}\beta_{Sn}}{\hat{\beta}_{Sn}^2 + \alpha_{Sn}^2} \frac{\sinh\left(\frac{\beta_{Sn}h}{2}\right)}{\sinh\left(\frac{\beta_{Sn}h}{2}\right)}, \\ \\ \begin{pmatrix} \phi_{An}(x) \\ \psi_{An}(x) \end{pmatrix} = \begin{pmatrix} -\sinh(\beta_{An}x) \\ B_{An} \cosh(\beta_{An}x) \end{pmatrix}, & B_{An} = \frac{2i\alpha_{An}\beta_{An}}{\hat{\beta}_{An}^2 + \alpha_{An}^2} \frac{\cosh\left(\frac{\beta_{An}h}{2}\right)}{\cosh\left(\frac{\beta_{An}h}{2}\right)} \end{array} \right. \qquad (15) $$

Here, $(\beta_{Sn}, \beta_{An})$ and their hats $(\hat{\beta}_{Sn}, \hat{\beta}_{An})$ are horizontal wavenumbers associated with the respective modes. These quantities satisfy the coupled dispersion relations for symmetric and antisymmetric Lamb waves:

$$ \left\{ \begin{array}{ll} \alpha_{Xn} = \sqrt{k_L^2 + \beta_{Xn}^2}, & \alpha_{\bar{X}n} = \sqrt{k_T^2 + \beta_{\bar{X}n}^2}, & X=A,S, \quad n=0, \dots, +\infty, \\ 4\alpha_{Sn}^2 \beta_{Sn} \hat{\beta}_{Sn} \sinh\left(\frac{\beta_{Sn}h}{2}\right) \cosh\left(\frac{\beta_{Sn}h}{2}\right) = (\alpha_{Sn}^2 + \beta_{Sn}^2)^2 \cosh\left(\frac{\beta_{Sn}h}{2}\right) \sinh\left(\frac{\beta_{Sn}h}{2}\right), & (16) \\ 4\alpha_{An}^2 \beta_{An} \hat{\beta}_{An} \cosh\left(\frac{\beta_{An}h}{2}\right) \sinh\left(\frac{\beta_{An}h}{2}\right) = (\alpha_{An}^2 + \beta_{An}^2)^2 \sinh\left(\frac{\beta_{An}h}{2}\right) \cosh\left(\frac{\beta_{An}h}{2}\right). & \end{array} \right. $$
<figure><img src="image_5.png" /><figcaption>Fig. 2. Domains used for (a) the numerical (direct) solution, see Eqs. (11) and (14), and (b) the effective model, see Eq. (23). In the latter, the cavity region $z \in (-d_a, -d)$, $d_a = d + a$, is not explicitly modeled, as its effect is accounted for by the homogenized jump conditions (4).</figcaption></figure>

**Truncation and matching conditions** To render the problem tractable numerically, we truncate the infinite series in Eqs. (11) and (14). Let:

$\mathcal{N}_1 = \{-n_1, \dots, n_1\}$ with $N_1 = (2n_1 + 1)$ the number of modes in $\mathcal{Y}^{\pm}$,

$\mathcal{N}_2 = \{0, \dots, n_2\}$ with $N_2 = (n_2 + 1)$ the number of Lamb modes in $\mathcal{Y}_L$.

We impose the continuity of displacements and normal stresses at the two horizontal interfaces $z = -d_a$ and $z = -d$, and the stress-free condition at $z = 0$. To do so, we project each continuity condition onto the truncated modal bases. We define the scalar products:

$$ (f, g)_1 = \int_{-h/2}^{h/2} f(x) \overline{g(x)} dx, \quad (f, g)_2 = \int_{-e/2}^{e/2} f(x) \overline{g(x)} dx, \qquad (17) $$

where the overline denotes complex conjugation. These inner products are used to enforce the interface and boundary conditions in weak form

$$ \begin{array}{ll} \text{at } z = -d_a, \quad m \in \mathcal{N}_2, \quad \left\{ \begin{array}{l} (u_x(x, -d_a^{-}), \phi_m)_2 = (u_x(x, -d_a^{+}), \phi_m)_2, \\ (u_z(x, -d_a^{-}), \psi_m)_2 = (u_z(x, -d_a^{+}), \psi_m)_2, \end{array} \right. & \text{for } \phi_m = \phi_{Am}, \phi_{Sm}, \\ & \text{for } \psi_m = \psi_{Am}, \psi_{Sm}, \end{array} \qquad (18) $$

$$ \begin{array}{ll} m \in \mathcal{N}_1, \quad \left\{ \begin{array}{l} (\sigma_{xz}(x, -d_a^{-}), F_m)_1 = (\sigma_{xz}(x, -d_a^{+}), F_m)_2, \\ (\sigma_{zz}(x, -d_a^{-}), F_m)_1 = (\sigma_{zz}(x, -d_a^{+}), F_m)_2, \end{array} \right. & \text{for } \phi_m = \phi_{Am}, \phi_{Sm}, \\ & \text{for } \psi_m = \psi_{Am}, \psi_{Sm}, \end{array} \qquad (19) $$

$$ \begin{array}{ll} \text{at } z = -d, \quad m \in \mathcal{N}_2, \quad \left\{ \begin{array}{l} (u_x(x, -d^{-}), \phi_m)_2 = (u_x(x, -d^{+}), \phi_m)_2, \\ (u_z(x, -d^{-}), \psi_m)_2 = (u_z(x, -d^{+}), \psi_m)_2, \end{array} \right. & \text{for } \phi_m = \phi_{Am}, \phi_{Sm}, \\ & \text{for } \psi_m = \psi_{Am}, \psi_{Sm}, \end{array} \qquad (19) $$

$$ m \in \mathcal{N}_1, \quad \left\{ \begin{array}{l} (\sigma_{xz}(x, -d^{-}), F_m)_2 = (\sigma_{xz}(x, -d^{+}), F_m)_1, \\ (\sigma_{zz}(x, -d^{-}), F_m)_2 = (\sigma_{zz}(x, -d^{+}), F_m)_1, \end{array} \right. \qquad (19) $$

$$ \text{at } z = 0, \quad m \in \mathcal{N}_1, \quad \left\{ \begin{array}{l} (\sigma_{xz}(x, 0), F_m)_1 = 0, \\ (\sigma_{zz}(x, 0), F_m)_1 = 0. \end{array} \right. \qquad (20) $$

These projections provide $(6N_1 + 8N_2)$ relations on the $6N_1$ unknown coefficients $(R_{Lh}, R_{Th}, b_{Lh}^{\pm}, b_{Th}^{\pm})$ in Eqs. (11) and $4N_2$ coefficients $(c_{Sh}^{\pm}, c_{Ah}^{\pm})$ in Eqs. (14). As in [24], the system is overdetermined and solved in the least-squares sense using QR decomposition (as implemented in MATLAB's backslash operator $\backslash$). In the following, we denote

$$ R_L = R_{L0}, \quad R_T = R_{T0}, \quad b_L^{\pm} = b_{L0}^{\pm}, \quad b_T^{\pm} = b_{T0}^{\pm}, \qquad (21) $$

the scattering coefficients associated with the propagating modes.

The numerical convergence of the multimodal formulation was verified by varying the truncation orders $N_1 = N_2 = N$. Comparisons with a reference solution obtained for $N = 100$ show that, over the frequency range considered (0–100) Hz, the scattering coefficients differ by less than 1% for $N = 5$, and by less than 0.3% for $N = 20$. All numerical results presented in this work were obtained with $N_1 = N_2 = 20$, which ensures a well-converged solution.
## 2.4. Resolution of the effective problem

The solution of the effective problem is formulated in two homogeneous domains:

$$Z = Z^{-} \cup Z^{+}, \quad Z^{-} = \{x \in (-\infty, +\infty), z \in (-\infty, -d_a)\}, \quad Z^{+} = \{x \in (-\infty, +\infty), z \in (-d, 0)\}, \qquad (22)$$

(see Fig. 2(b)). The intermediate layer $z \in (-d_a, -d)$ is not explicitly modeled, as the physical effect of the periodic cavities has been encapsulated into the jump conditions introduced in Eq. (4). This effective formulation restores translational invariance in the horizontal direction $x$ as the periodic structure is no longer present. Consequently, no evanescent (non-propagating) modes are included in the solution: their influence has been entirely accounted for through the effective interface model. Thus, the solution reads

$$r \in Z^{-}: \quad \begin{pmatrix} \phi(\mathbf{r}) \\ \psi(\mathbf{r}) \end{pmatrix} = \begin{pmatrix} \phi^{\text{inc}}(\mathbf{r}) \\ \psi^{\text{inc}}(\mathbf{r}) \end{pmatrix} + \begin{pmatrix} R_L e^{-i\alpha_L z} \\ R_T e^{-i\alpha_T z} \end{pmatrix} e^{i\beta x}, \qquad (23)$$

$$r \in Z^{+}: \quad \begin{pmatrix} \phi(\mathbf{r}) \\ \psi(\mathbf{r}) \end{pmatrix} = \begin{pmatrix} b_L^+ e^{i\alpha_L z} + b_L^- e^{-i\alpha_L z} \\ b_T^+ e^{i\alpha_T z} + b_T^- e^{-i\alpha_T z} \end{pmatrix} e^{i\beta x}, \qquad (24)$$

where $(\alpha_L, \alpha_T, \beta)$ are defined in (9). The six scattering coefficients, being the reflection coefficients $(R_L, R_T)$ in $Z^{-}$ and amplitudes of the upward $(b_L^+, b_T^+)$ and downward $(b_L^-, b_T^-)$ propagating waves in $Z^{+}$ are determined by enforcing the two stress-free conditions at the free surface (2) and the four effective jump conditions (4). Inserting the expressions from (23) into these interface and boundary conditions yields a linear system of six equations for the six unknowns. This system can be written in compact form as:

$$mC = S, \qquad (25)$$

where $m$ is a 6 × 6 matrix, $C$ is the vector of the unknowns coefficients and $S$ contains the source terms arising from the incident wave. Specifically, denoting $e_{L,T} = e^{i\alpha_{L,T}d}$, $e_{L,T}^* = e^{i\alpha_{L,T}d_a}$ and introducing

$$\xi = \frac{\sin \theta_L}{\sin \theta_T} = \frac{k_T}{k_L} = \sqrt{\frac{2(1-\nu)}{(1-2\nu)}}, \qquad (25)$$

they read

$$m = \begin{pmatrix} m_{11} & -\overline{m_{11}} & m_{11} & m_{14} & \overline{m_{14}} & -m_{14} \\ m_{21} & \overline{m_{21}} & -m_{21} & m_{24} & -\overline{m_{24}} & m_{24} \\ m_{31} & \overline{m_{31}} & -m_{31} & m_{34} & -\overline{m_{34}} & m_{34} \\ m_{41} & -\overline{m_{41}} & m_{41} & m_{44} & \overline{m_{44}} & -m_{44} \\ m_{51} & m_{52} & 0 & m_{54} & m_{55} & 0 \\ m_{61} & m_{62} & 0 & m_{64} & m_{65} & 0 \end{pmatrix}, \quad C = \begin{pmatrix} b_L^+/e_L \\ b_L^-/e_L \\ R_L e_L^* \\ b_T^+/e_T^{-1} \\ b_T^-/e_T \\ R_T e_T^* \end{pmatrix}, \quad S = \begin{pmatrix} \overline{m_{11}} \\ \overline{m_{21}} \\ \overline{m_{31}} \\ m_{41} \\ 0 \\ 0 \end{pmatrix} A_L^{\text{inc}}/e_L^* + \begin{pmatrix} \overline{m_{14}} \\ \overline{m_{24}} \\ \overline{m_{34}} \\ \overline{m_{44}} \\ 0 \\ 0 \end{pmatrix} A_T^{\text{inc}}/e_T^*, \qquad (26)$$

where

$$m_{11} = \cos \theta_L - \frac{i k_L h}{2} (B_z \cos^2 \theta_L + B_x \sin^2 \theta_L),$$

$$m_{21} = \sin \theta_T \left(1 - \frac{i \alpha_L}{2}(2hB - a)\right),$$

$$m_{31} = \cos 2\theta_T \left(1 - \frac{i \alpha_L a}{2}\right) + \frac{i \omega \alpha_L a}{2},$$

$$m_{41} = -\sin 2\theta_L \left(1 - \frac{i(1-\varphi)\alpha_L a}{2}\right) - \frac{i\beta h}{2} (C_z \cos^2 \theta_L + C_x \sin^2 \theta_L),$$

$$m_{51} = -\sin 2\theta_L e_L, \quad m_{52} = \sin 2\theta_L / e_L, \quad m_{54} = \xi^2 \cos 2\theta_L e_T,$$

$$m_{61} = \cos 2\theta_L e_L, \quad m_{62} = \cos 2\theta_L / e_L, \quad m_{64} = \sin 2\theta_L e_T, \quad m_{65} = -\sin 2\theta_L / e_T.$$

$$m_{14} = \sin \theta_L \left(1 - \frac{i \alpha_T h}{2} (B_z - B_x)\right),$$

$$m_{24} = -\cos \theta_T + \frac{i k_T}{2} (hB \cos 2\theta_T + a \sin^2 \theta_T)$$

$$m_{34} = \sin 2\theta_T \left(1 - \frac{i \alpha_T a}{2}\right) + \frac{i \varphi \beta a}{2},$$

$$m_{44} = \xi^2 \cos 2\theta_T \left(1 - \frac{i(1-\varphi)\alpha_T a}{2}\right) - \frac{i \alpha_T h}{2} (C_z - C_x) \sin^2 \theta_L,$$

$$m_{55} = \xi^2 \cos 2\theta_T / e_T,$$

$$m_{65} = -\sin 2\theta_L / e_T.$$

This system is fully explicit, providing a closed-form solution for the effective problem.

## 2.5. The free-field solution in the absence of cavities

As previously mentioned, the free-field solution refers to the configuration without cavities — that is, the solution of the elastic wave problem posed in the homogeneous half-space $Z = \{x \in (-\infty, +\infty), z \in (-\infty, 0)\}$, governed by (1) and (2). In this case, an explicit analytical solution is available, which describes the incident wave and its reflection from the free surface. The reflected field is expressed in terms of reflection coefficients that depend on the incidence angle and the type of incident wave (L or T), namely

$$r \in Z: \quad \begin{pmatrix} \phi^{\text{free}}(\mathbf{r}) \\ \psi^{\text{free}}(\mathbf{r}) \end{pmatrix} = \begin{pmatrix} \phi^{\text{inc}}(\mathbf{r}) \\ \psi^{\text{inc}}(\mathbf{r}) \end{pmatrix} + \begin{pmatrix} R_L^{\text{free}} e^{-i\alpha_L z} \\ R_T^{\text{free}} e^{-i\alpha_T z} \end{pmatrix} e^{i\beta x}, \qquad (27)$$

with

$$R_L^{\text{free}} = r_{LL}^{\text{free}} A_L^{\text{inc}} + r_{LT}^{\text{free}} A_T^{\text{inc}}, \quad R_T^{\text{free}} = r_{TT}^{\text{free}} A_L^{\text{inc}} + r_{TT}^{\text{free}} A_T^{\text{inc}}, \qquad (28)$$

and

$$r_{LL}^{\text{free}} = r_{TT}^{\text{free}} = \frac{\sin 2\theta_L \sin 2\theta_T - \xi^2 \cos^2 2\theta_T}{\sin 2\theta_L \sin 2\theta_T + \xi^2 \cos^2 2\theta_T}, \quad r_{LT}^{\text{free}} = \frac{2 \sin 2\theta_L \cos 2\theta_T}{\sin 2\theta_L \sin 2\theta_T + \xi^2 \cos^2 2\theta_T}, \quad r_{TT}^{\text{free}} = -\frac{2\xi^2 \sin 2\theta_T \cos 2\theta_T}{\sin 2\theta_L \sin 2\theta_T + \xi^2 \cos^2 2\theta_T},$$
<caption>Table 1</caption>

Values of the effective coefficients entering in (4) for arrays characterized by $a/h = 3/5$ and increasing density $\varphi = e/h$.

<table><thead><tr><th>e/h</th><th>B<sub>z</sub></th><th>B</th><th>B<sub>x</sub></th><th>C<sub>z</sub></th><th>C<sub>x</sub></th></tr></thead><tbody><tr><td>2/5</td><td>1.40</td><td>2.57</td><td>0.40</td><td>-0.16</td><td>-0.87</td></tr><tr><td>4/5</td><td>4.98</td><td>23.38</td><td>1.63</td><td>-0.03</td><td>-0.32</td></tr><tr><td>4.8/5</td><td>20.73</td><td>1500</td><td>6.91</td><td>-0.005</td><td>-0.06</td></tr></tbody></table>

(29)

(Chapter 5 in [23]). This reference solution is naturally recovered from the effective model (24)–(26) by setting $d = d_a$ (i.e., $a = 0$), $e = 0$, and letting $B_i = C_i = 0$, $i = x, z$, $B = 0$. This yields the free-field solution with $b_L^+ = A_L^{\text{inc}}$, $b_T^+ = A_T^{\text{inc}}$, $b_L^- = R_L^{\text{free}}$ and $b_T^- = R_T^{\text{free}}$.

### 3. Validation of the homogenized model: Shielding effects and resonances

From now on, we consider a typical shale-like elastic material characterized by the following properties: $\lambda = \mu = 3.10^9$ Pa, $\rho = 2000$ kg.m$^{-3}$, where $(\lambda, \mu)$ are the Lamé parameters, and $\rho$ is the density.

The dimensions of the cavities and their burial depths are on the order of a few meters, with potential applications to anthropogenic underground cavities in mind. In this section, we explore frequency ranges up to 100 Hz in order to gain a comprehensive understanding of the underlying phenomena. A more focused analysis in the seismic frequency band (0–10) Hz is conducted in Section 4.

Throughout this section, the displacement fields $(u_x, u_z)$ are computed for incident waves of the form (8), with coefficients $(A_L^{\text{inc}}, A_T^{\text{inc}})$ chosen to produce a unit-amplitude displacement at the free surface ($z = 0$), namely:

$$u_x^{\text{inc}2} + u_z^{\text{inc}2} = k_L^2 A_L^{\text{inc}2} + k_T^2 A_T^{\text{inc}2} + 2A_L^{\text{inc}} A_T^{\text{inc}} (\alpha_L - \alpha_T) = 1. \quad (30)$$

#### 3.1. A typical example

We begin by analyzing the wavefield response for cavity arrays with varying cavity density $e = 2, 4$ and $4.8$ m, within a fixed lattice period $h = 5$ m and height $a = 3$ m. The array is buried at depth $d = 40$ m.

Fig. 3 shows displacement fields $u_x(x, z)$ and $u_z(x, z)$ at a frequency of 30 Hz for an incident wave at angle $\theta_L = 40^\circ$ (corresponding to $\theta_T = 21.8^\circ$), with $A_L^{\text{inc}} = A$, $A_T^{\text{inc}} = 2A$, where A is a normalization factor ensuring unit surface amplitude (30). Both numerical simulations and the analytical effective model are shown. Let us recall that, in the actual (full) problem, the near field in the vicinity of the cavities consists of a multitude of evanescent modes, whereas in the homogenized model, their effect has been captured through the jump conditions. In this latter case, the displacements are given by

$$ \begin{aligned} z \in (-\infty, -d_a), \quad \left\{ \begin{array}{l} u_x(x, z) = i \left[ \beta (A_L^{\text{inc}} e^{i\alpha_L z} + R_L e^{-i\alpha_L z}) - \alpha_T (A_T^{\text{inc}} e^{i\alpha_T z} - R_T e^{-i\alpha_T z}) \right] e^{i\beta x}, \\ u_z(x, z) = i \left[ \alpha_L (A_L^{\text{inc}} e^{i\alpha_L z} - R_L e^{-i\alpha_L z}) + \beta (A_L^{\text{inc}} e^{i\alpha_L z} + R_L e^{-i\alpha_L z}) \right] e^{i\beta x}, \end{array} \right. \\ \\ z \in (-\infty, -d_a), \quad \left\{ \begin{array}{l} u_x(x, z) = i \left[ \beta (b_L^+ e^{i\alpha_L z} + b_L^- e^{-i\alpha_L z}) - \alpha_T (b_T^+ e^{i\alpha_T z} - b_T^- e^{-i\alpha_T z}) \right] e^{i\beta x}, \\ u_z(x, z) = i \left[ \alpha_L (b_L^+ e^{i\alpha_L z} - b_L^- e^{-i\alpha_L z}) + \beta (b_L^+ e^{i\alpha_L z} + b_L^- e^{-i\alpha_L z}) \right] e^{i\beta x}, \end{array} \right. \end{aligned} \quad (31) $$

where the coefficients $(R_L, R_T)$ and $(b_L^\pm, b_T^\pm)$ are given by (24)–(26) and in the present case, the values of effective parameters are given in Table 1.

The comparison shows excellent qualitative agreement between the numerical and effective solutions. Additionally, the effect of the cavity array is illustrated by comparing with the free-field solution (no cavities) $(u_x^{\text{free}}, u_z^{\text{free}})$. It is observed that the density of the cavity array, expressed by the ratio $\varphi = e/h$, significantly affects the soil response. Increasing it enhances wave reflection by the array and thus, reduces displacements in the top layer — defined as the region between the surface and the top of the cavities. This indicates a shielding effect, where the cavity array acts as a seismic shield by partially blocking wave energy propagation into the soil beneath.

Going further, Fig. 4 presents a more quantitative analysis of the reflection coefficients $(R_L, R_T)$ and the upward and downward propagating wave amplitudes $(b_L^\pm, b_T^\pm)$ within the top layer, across frequencies and for the three values of $e$. At zero frequency, all configurations recover the free-field values, which remain constant across frequencies, since the cavity array is transparent to long wavelengths. These reference values are shown with yellow dashed lines for $(R_L, R_T)$.

The curves show both the numerical results (solid lines) and the effective model predictions (black dashed lines). The agreement is excellent overall, although a slight decrease in accuracy is noted at higher frequencies. This is expected because the model has a range of validity in the low frequency regime, measured by low values of $k_L h, k_T h$. In the reported curves, with $f$ up to 60 Hz, these parameters reach values up to 0.9 and 1.5.

Importantly, the previous qualitative trend (denser arrays yield more shielding) is now nuanced: resonance peaks are clearly visible. These resonances are more pronounced with increasing cavity density. Yet between resonances, the amplitudes $(b_L^\pm, b_T^\pm)$ in the top layer remain well below the free-field case, further confirming the shielding effect.

We conclude with two remarks of different nature:
<figure><img src="image_6.png" /><figcaption>Fig. 3. Real parts of the displacement fields $u_x(x, z)$ and $u_z(x, z)$ — Left panels show the free-fields (no cavities); the other panels show the fields for a periodic array buried at depth $d = 40$ m with $a = 3$ m, $h = 5$, and increasing cavity densities (a) $e = 2$ m, (b) $e = 4$ m, and (c) $e = 4.8$ m. Results are shown for frequency $f = 30$ Hz over the domain {$x \in (0, 140)$ m, $z \in (-120, 0)$ m}. A shielding effect is visible, reducing displacement amplitudes in the top layer.</figcaption></figure>

– First, we note in Fig. 4(c) the presence of a resonance (marked with yellow arrows) for $e = 4.8$ m at $f = 40$ Hz, which is not predicted by the model. This resonance corresponds to the flexural resonance of the thin elastic ligament, with slenderness $d/(h-e)$. Accounting for such resonances can be done, see [25,26]. However, it lies beyond the scope of this study and would make the model less tractable. Hence in the present study, it has not been done.

– Next, in the plots of Fig. 3, the observed resonances have been identified as resonances of the top layer. These correspond to a class of resonances whose limiting case occurs when $e = h$, i.e., when the top layer is completely isolated from the half-space beneath, subject to stress-free boundary conditions at $z = 0$ and $z = -d$. For arrays with $e < h$, leakage occurs and the top layer is not perfectly isolated; resonances are excited, increasingly resembling perfect resonances as the cavity density grows. However, the associated resonance frequencies remain significantly shifted from the nominal values, expected at

$$ \alpha_L d \approx n\pi, \quad \alpha_T d \approx n\pi, \qquad (32) $$

with $n$ being an integer. Indeed, in addition to the resonance expected close to zero frequency, these estimates give a longitudinal resonance at $f = 37$ Hz and three transverse resonances at $f = 17.8$, 35.6 and 53.5 Hz, which is quite far from the resonance frequencies observed in Fig. 4. This is analyzed in the forthcoming section.
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Scattering coefficients against frequency *f* for the same configurations as in Fig. 3. Solid lines: direct numerical results, (21); dashed black lines: effective model predictions using (24)–(26), see also (31).</figcaption></figure>

## 3.2. Analysis of the resonances

We now turn to a more detailed analysis of the resonances associated with the top layer. In this subsection, we consider a fixed cavity configuration with *e* = 4 m, periodicity *h* = 5 m and cavity height *a* = 3 m, as in Figs. 3(b) and 4(b). However, we now reduce the burial depth to *d* = 20 m and extend the frequency range up to 100 Hz. This configuration preserves the five resonances previously observed and allows us to test the robustness of the effective model beyond its strict asymptotic regime: the dimensionless frequencies reach values up to *k*<sub>T</sub>*h* = 1.5 and *k*<sub>L</sub>*h* = 2.6, thus probing the limits of validity of the homogenized model.

Figs. 5 and 6 show the amplitudes of the surface displacements *u*<sub>x</sub>(*x*, 0) and *u*<sub>z</sub>(*x*, 0) as functions of frequency and angle of incidence, revealing the trajectories of the top-layer resonances. For clarity, we distinguish between the cases of purely longitudinal and purely transverse incident waves, characterized by angles of incidence *θ*<sub>L</sub> and *θ*<sub>T</sub>, respectively and frequencies are made dimensionless with respect to the corresponding wavenumbers *k*<sub>L</sub> and *k*<sub>T</sub>. In the absence of cavities, the amplitudes of surface displacement are independent of frequency and depend only on the angle of incidence. On the color maps shown in Figs. 5 and 6 (with cavities present), the free-field displacement values correspond to the zero-frequency limit and serve as a baseline for comparison.

We now comment on the resonance features visible in the color plots. For a longitudinal incident wave at normal incidence (*θ*<sub>L</sub> = 0), the displacements are purely vertical (i.e., *u*<sub>x</sub> = 0). In this case, two clear longitudinal resonances appear in the vertical displacement component *u*<sub>z</sub>, corresponding to the conditions *k*<sub>L</sub>*d* ≃ *nπ*, with *n* = 0 and 1, as in (32). Conversely, for a transverse
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Surface displacement amplitudes $|u_x(x, 0)|$ and $|u_z(x, 0)|$ for an incident longitudinal (L) wave ($A_L^{\text{inc}} = 0$), as functions of the angle of incidence $\theta_L$ and of the dimensionless frequency $k_L d/\pi$. (a) Numerical solution of the full problem. (b) Prediction from the effective model (31), using coefficients from (24)–(26).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Surface displacement amplitudes $|u_x(x, 0)|$ and $|u_z(x, 0)|$ for an incident transverse (T) wave ($A_L^{\text{inc}} = 0$), as functions of the angle of incidence $\theta_T$ and of the dimensionless frequency $k_T d/\pi$. Same representation as in Fig. 5. The range $\theta_T \in (0, \theta_T^c = 35^\circ)$ corresponds to $\theta_L \in (0, 90^\circ)$.</figcaption></figure>

incident wave at normal incidence ($\theta_T = 0$), the displacements are purely horizontal (i.e., $u_x = 0$), and four transverse resonances are observed, following the conditions $k_T d \approx n\pi$ for $n = 0,1,2,3$. As the incidence becomes oblique in either case, mode conversion between longitudinal and transverse waves occurs, causing both types of resonances to appear, regardless of the nature of the incident wave. In the case of a T-incident wave, one clearly observes the critical angle for T-L conversion at $\theta_T^c = 35^\circ$ (corresponding to $\theta_L = 90^\circ$), beyond which longitudinal waves cannot propagate. Consequently, the associated L-wave resonance disappear.

It can be seen that the resonance trajectories observed in the full numerical simulations do not simply follow the $1/\cos\theta$ scaling predicted by the half-wavelength resonance criterion in (32), for $\theta = \theta_L$ or $\theta_T$. Instead, they exhibit complex interactions that are generally well captured by the effective model. As in previous comparisons, the accuracy gradually decreases at higher frequencies,
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Finite size effect — Real part (a) and absolute value (b) of the vertical displacement field $u_z(x, z)$, normalized by $u_z^{\text{free}}(0, 0)$, obtained in the absence of any array under the same source conditions, at frequency $f = 14$ Hz. The left three panels show the results from Comsol for $N = 1, 5$ and $11$ cavities with a point-source at $z = -80$ m; the rightmost panels show the wavefields obtained for an infinite array. The domain is $\{x \in (-40, 40) \text{ m}, z \in (-100, 0) \text{ m}\}$.</figcaption></figure>

consistent with the asymptotic nature of the homogenized formulation. Moreover, the effective model does not accurately reproduce the regions of avoided crossing, as these result from subtle coupling mechanisms that would only be captured by higher-order corrections in the homogenized model, see [27].

### 3.3. Finite size effects: From single to multiple cavities

We now turn to the question of finite-size effects. Thus far, our effective model has been validated for an infinite periodic array of cavities. This configuration has the advantage of being amenable to an explicit, analytical treatment, which is both computationally efficient and insightful, especially when compared to full numerical simulations that can become computationally expensive. However, in real-world-applications – particularly in seismic contexts – arrays of underground cavities are always of finite extent, even though many practical configurations, such as mining galleries, may have a large extent with a large number of cavities. It is therefore important to assess whether the effects predicted for an infinite array persist in finite configurations, a question that we address in this section by considering cavity arrays with $N = 1, 5$ and $11$ cavities.

The corresponding configurations are illustrated in Figs. 7 and 8, for two representative frequencies: the lowest resonance at 14 Hz (relevant for seismic applications) and a higher resonance at 55 Hz. These simulations use the same network geometry as in the previous sections, with $a = 3$ m, $e = 4$ m, $h = 5$ m, and $d = 20$ m. We report the vertical displacement fields, both in real part and in magnitude. For the infinite array (right panels), the wavefields are computed using the modal method with a normally incident L plane wave. For the finite-size arrays, the fields are computed with Comsol, using a quasi-point source located 60 m below the array (i.e., at $z = -80$ m), centered at $x = 0$ and emitting isotropically. Perfectly Matched Layers (PMLs) are applied on the vertical and bottom boundaries of the computational domain, placed sufficiently far from both the array and the source. To allow meaningful comparisons, all displacements are normalized by the free-field vertical surface displacement $u_z^{\text{free}}(0, 0)$, computed in the absence of any cavities, under the same source conditions. (Note that for a plane wave, the free-field amplitude at $z = 0$ is constant along $x$, whereas for a point-source it depends on $x$ and is maximum at $x = 0$).

Qualitatively, the results show that the amplification effect within the top layer is already visible for as few as $N = 5$ cavities, with amplitudes comparable to those observed in the infinite case. This suggests that the essential features of amplification and shielding observed in finite cavity arrays are well described by the infinite-array model.
<figure><img src="image_11.png" /><figcaption>Fig. 8. Finite size effect — same representation as in Fig. 7 for frequency $f = 55$ Hz.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 9. Finite size effect — Amplification profiles at (a) $f = 14$ Hz and (b) $f = 55$ Hz for $N = 0$ (free-field) and $N = 1, 5$ and $11$ cavities, from Figs. 7 and 8. Panel (c) shows the amplification at $x = 0$ against frequency for the finite size arrays ($N = 1, 5, 11$) and for reference, the case of an infinite array (blue line from numerics and dashed black line from the model). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

These observations are further supported by Fig. 9. Panel (c) shows the vertical amplification at the center point ($x = z = 0$) as a function of frequency for the different configurations: $N = 1, 5, 11$, and the infinite array. Panels (a) and (b) show the profile along $x$ of the surface amplification at the first two resonances, which help illustrate the influence of wavefront curvature caused by the point source and diffraction effects at the array edges — effects that are absent in the plane wave configuration. The free-field responses, shown in gray, provide a baseline for comparison. In Panel (c), the variation of amplification with frequency is smoother in the infinite case but remains very similar in terms of both resonance peaks and amplitude levels for $N = 5$ and $N = 11$ (note that, for $N = 11$, a small additional peak appears around 90 Hz, corresponding to a transverse resonance arising from mode conversion between L- and T-waves; this conversion becomes non-negligible at high frequencies).
<figure><img src="image_13.png" /><figcaption>Fig. 10. Horizontal amplification factor $U_x$, Eq. (33), against frequency $f$ for different incidence angles: (a) incident L wave and (b) incident T wave. In panel (b), the T-wave incidence angles $\theta_T$ are chosen such that the corresponding $\theta_L$ = angles match those in (a): 0, 15°, 30° and 45°.</figcaption></figure>

Finally, we note that the lowest resonance remains especially stable across configurations, both in terms of resonance frequency and associated, maximum, amplitude. This makes it a particularly promising feature for applications in seismic shielding or vibration mitigation. This lowest resonance will be examined more closely in the following section.

## 4. Earthquake engineering: Focus on first resonance

We now turn to the study of the first resonance, located in the low-frequency range (0-10) Hz, which is particularly relevant for earthquake engineering applications. As observed in Fig. 9, this resonance exhibits remarkably stable features when transitioning from an infinite to a finite cavity array. In this section, we analyze the sensitivity of its key characteristics – namely the resonant frequency and associated maximum amplification – with respect to the geometry of the cavities (parameters $a$ and $e$), while keeping the period $h = 5$ m constant. We also examine the influence of the burial depth $d$. Since horizontal ground motion plays a major role in evaluating potential structural damage, we focus our attention on the horizontal displacement component $u_x$. Throughout this section, we refer to the horizontal amplification factor, $U_x$, defined as the ratio between the displacement amplitude in the presence of cavities and the free-field amplitude, at the free surface:

$$U_x = \frac{|u_x(x, 0)|}{|u_x^{\text{free}}(x, 0)|} \qquad (33)$$

We begin by fixing the parameters $a = 3$ m, $e = 4$ m, and $d = 20$ m, and explore the influence of the angle of incidence for both longitudinal (L) and transverse (T) incident waves. The horizontal displacement amplification is shown in Fig. 10. For incidence angles up to 45°, the maximum amplification remains nearly constant, around 2.5 in this configuration. The resonance frequency shifts only slightly – from 4 Hz to 5 Hz – with the most notable variations observed above 30°.

This observation is key, as it suggests that a simple estimate of the resonance frequency and the associated amplification can be obtained under the assumption of normal incidence. We therefore focus on a normally incident T wave, for which the horizontal displacement in the top layer can be approximated by:

$$u_x(\mathbf{r}) \approx \frac{1 - ik_T d}{1 - ik_T d - k_T^2 d h \mathbf{B}}, \quad u_x^{\text{free}}(\mathbf{r}) = 2 \cos k_T z. \qquad (34)$$

From this approximation, we derive the following explicit expressions for the resonance frequency $f_R$ and the associated maximum amplification $U_R$:

$$f_R = \frac{c}{2\pi} \frac{1}{\sqrt{h\mathbf{B}d}}, \quad U_R = \sqrt{\frac{h\mathbf{B}}{d}}. \qquad (35)$$

These formulas reveal that both $f_R$ and $U_R$ decrease proportionally to $1/\sqrt{d}$, highlighting the impact of burial depth. The influence of cavity geometry enters through the term $h\mathbf{B}$, which increases with the periodicity $h$ and both the cavity aspect ratios $e/h$ and $a/h$. Hence, in simple terms: the larger the void volume, the lower the resonance frequency and the stronger the amplification.

These trends are confirmed in Fig. 11, which reports the amplification factor $U_x$ as a function of frequency while varying $a$ and $e$ separately, and then varying $d$. In all the reported cases, the model accurately reproduces these variations (the corresponding effective coefficients are listed in Tables 2 and 3). While the agreement between the effective model and full numerics is excellent in Fig. 11, the model used here is the full homogenized model (from (24)–(26)), not the simplified approximation in (35). To test the validity of the latter – which is intended for rapid estimates – we report in Fig. 12 the variations of $f_R$ and $U_R$ computed via direct numerical simulations (solid blue lines) and the simplified analytical prediction in (35) (black lines). The predictions of (35)
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Sensitivity of the horizontal amplification factor $U_x$, Eq. (33), to geometric parameters: (a) cavity height $a$, (b) array density $e$ and (c) burial depth $d$. Colored lines: direct numerical simulations; dashed black lines: predictions from the effective model.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Variations of the resonance frequency $f_R$ (top rows) and the corresponding horizontal amplification $U_R$ (bottom rows) as functions of: (a) cavity height $a$, (b) array density $e$ and (c) burial depth $d$. Blue lines: results from direct numerics using profiles as reported in Fig. 11. Black lines: simplified estimates from (35). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<caption>**Table 2**</caption>

Values of the effective coefficients entering in (4) for arrays characterized by $e/h = 3/5$ and increasing $a/h$.

<table><thead><tr><th>a/h</th><th>1/5</th><th>2/5</th><th>3/5</th><th>4/5</th><th>5/5</th><th>6/5</th><th>7/5</th><th>8/5</th><th>9/5</th></tr></thead><tbody><tr><th>B<sub>z</sub></th><td>1.29</td><td>1.84</td><td>2.41</td><td>2.97</td><td>3.53</td><td>4.09</td><td>4.66</td><td>5.22</td><td>5.78</td></tr><tr><th>B</th><td>1.25</td><td>2.91</td><td>5.57</td><td>9.51</td><td>15.01</td><td>22.34</td><td>31.79</td><td>43.66</td><td>58.20</td></tr><tr><th>B<sub>x</sub></th><td>0.38</td><td>0.56</td><td>0.75</td><td>0.94</td><td>1.13</td><td>1.31</td><td>1.50</td><td>1.69</td><td>1.88</td></tr><tr><th>C<sub>z</sub></th><td>0.06</td><td>-0.007</td><td>-0.09</td><td>-0.17</td><td>-0.25</td><td>-0.33</td><td>-0.41</td><td>-0.49</td><td>-0.57</td></tr><tr><th>C<sub>x</sub></th><td>-0.13</td><td>-0.37</td><td>-0.61</td><td>-0.85</td><td>-1.09</td><td>-1.33</td><td>-1.57</td><td>-1.81</td><td>-2.05</td></tr></tbody></table>

are reasonably accurate as long as the resonance remains in the low-frequency regime, which is consistent with the assumptions $k_T d \ll 1$ used to derive them (this condition is more restrictive than $k_T h \ll 1$, which governs the validity of the effective model itself). Nonetheless, these simplified formulas remain useful for quickly estimating the impact of geometric parameters on resonance characteristics, especially in early-stage seismic design or assessment.

## 5. Conclusion

In this work, we have investigated the seismic response of a shallow subsurface layer containing a periodic array of voids, focusing on the development and validation of an effective model based on homogenized jump conditions. Starting from a reference full-space elasticity problem, we showed that regularly spaced cavities can significantly affect wave propagation, notably through
<table><caption>Table 3</caption><caption>Values of the effective coefficients entering in (4) for arrays characterized by a/h = 3/5 and increasing density φ = e/h.</caption><thead><tr><th>e/h</th><th>1/5</th><th>1.5/5</th><th>2/5</th><th>2.5/5</th><th>3/5</th><th>3.5/5</th><th>4/5</th><th>4.5/5</th><th>4.8/5</th></tr></thead><tbody><tr><th>B<sub>z</sub></th><td>0.89</td><td>1.10</td><td>1.40</td><td>1.81</td><td>2.41</td><td>3.23</td><td>4.98</td><td>9.32</td><td>20.73</td></tr><tr><th>B</th><td>1.45</td><td>1.90</td><td>2.57</td><td>3.64</td><td>5.57</td><td>9.84</td><td>23.38</td><td>124.48</td><td>1500</td></tr><tr><th>B<sub>x</sub></th><td>0.23</td><td>0.30</td><td>0.40</td><td>0.54</td><td>0.75</td><td>1.07</td><td>1.63</td><td>3.09</td><td>6.91</td></tr><tr><th>C<sub>z</sub></th><td>-0.27</td><td>-0.21</td><td>-0.16</td><td>-0.12</td><td>-0.09</td><td>-0.06</td><td>-0.03</td><td>-0.01</td><td>-0.005</td></tr><tr><th>C<sub>x</sub></th><td>-1.09</td><td>-0.99</td><td>-0.87</td><td>-0.74</td><td>-0.61</td><td>-0.47</td><td>-0.32</td><td>-0.16</td><td>-0.06</td></tr></tbody></table>

partial reflection and the emergence of resonance phenomena in the soil layer between the surface and the top of the cavity network. The proposed effective model, derived through asymptotic homogenization, was validated against full numerical simulations and shown to accurately reproduce both the near-field displacement fields and far-field scattering behavior over a broad frequency range. Despite being developed under low-frequency assumptions (long-wavelength excitation), the model remains robust up to moderate dimensionless frequencies ($k_T h \approx 1.5$, $k_L h \approx 2.6$), providing a substantial computational advantage compared to direct simulations.

Particular attention was devoted to the first resonance, which lies within the seismic frequency range and was found to be especially robust and predictable, even when transitioning from infinite to finite-size cavity arrays. A simple formula was proposed to estimate both the resonance frequency and the associated amplification, with explicit dependence on burial depth and cavity geometry. This provides a practical tool for assessing the impact of shallow subsurface cavity networks, particularly in contexts involving seismic site effects or anthropogenic voids such as abandoned mines.

The present formulation, however, relies on several simplifying assumptions that delimit its range of validity. The jump conditions are derived for an infinite and perfectly periodic array and under long-wavelength conditions ($k_T h \ll 1$, $k_L h \ll 1$) within the framework of linear elasticity, neglecting soil damping, stratification, and potential nonlinear effects. Moreover, the analysis is restricted to a two-dimensional geometry. Numerical comparisons with finite arrays show that the main amplification and resonance features persist for realistic configurations, but a full theoretical treatment of finite-size and disorder effects remains to be developed.

Beyond these assumptions, several natural extensions can be envisioned. The present study focused on in-plane body waves, whereas a source located at the free surface would also generate Rayleigh waves that interact with the cavity network, in analogy with the interaction of surface waves with arrays of cylindrical inclusions in soft soils [18,28]. Future work will also aim to extend the present framework to configurations more representative of engineering practice, where cavities are lined or reinforced to prevent collapse. In a three-dimensional setting (e.g., tunnels, galleries, or interconnected voids), the homogenized interface approach can be enriched to include thin elastic linings, compliant interphases, or partially filled cavities, thereby accounting for the additional stiffness and damping introduced by support systems. Such developments will not only improve the realism of the model but also contribute to the broader field of seismic-wave control through engineered subsurface structures, at the interface between seismic metamaterials and earthquake engineering.

### CRediT authorship contribution statement

**Agnès Maurel:** Writing – review & editing, Writing – original draft, Methodology, Formal analysis, Software, Numerical investigation, Validation, Visualization. **Stéphane Brulé:** Writing – review & editing, Conceptualization, Methodology, Validation. **Sébastien Guenneau:** Writing – review & editing, Methodology, Validation. **Kim Pham:** Writing – review & editing, Methodology, Formal analysis, Software, Numerical investigation, Validation.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

Kim Pham acknowledges support from the French Agence Nationale de la Recherche (ANR) under grant ANR-19-CE08-0006 MetaReso.

### Appendix A. Elementary problems

The jump conditions presented in Eq. (4) are derived from the asymptotic homogenization framework developed in [22], with simplifications arising from the fact that the cavities considered here are rectangular and thus symmetric with respect to both the x- and z-axes. For completeness, we recall the notation introduced in (4), namely ($B, B_x, B_z, C_x, C_z$), and we detail below the elementary problems from which these effective coefficients are computed. The elementary problems are formulated in rescaled coordinates $r_m = r/h$, over the domain:

$$ \mathcal{Y}_p = \{r_m = (x_m, z_m), x_m \in (-1/2, 1/2), z_m \in (-\infty, +\infty)\}, $$
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 13. Effective coefficients as functions of the filling ratio $\varphi$ for $a/h = 3/5$, see also Table 3. (a) Diverging coefficients ($B, B_x, B_z$) plotted in logarithmic scale versus $(1-\varphi)$; dashed lines indicate the fitted power-law trends $B \sim (1-\varphi)^{-2.5}$ and $B_x = \lambda B_z/(\lambda + 2\mu) \sim (1-\varphi)^{-1}$, in agreement with (41) for $\varphi \to 1$. (b) Vanishing coefficients ($C_x, C_z$) as functions of $(1-\varphi)$; dashed lines show the fitted dependences $C_x \sim (1-\varphi)$ and $C_z \sim (1-\varphi)$.</figcaption></figure>

where $C$ is the region of the cavity. For each pair of directions $i, j \in \{x, z\}$, we introduce the fields ($V^{ij}, \Sigma^{i,j}$), which are periodic in $x_m$ and satisfy the following static problem:

$$ \left\{ \begin{array}{l} \mathrm{div} \Sigma^{ij} = 0, \quad \Sigma^{ij} = A(\epsilon(V^{ij}) + e^{ij}), \quad \text{in } \mathcal{Y}_p, \\ \Sigma^{ij} n = 0 \quad \text{on the cavity boundaries,} \\ \lim_{z_m \to \pm\infty} \epsilon(V^{ij}) = 0. \end{array} \right. \qquad (36) $$

where $A(\epsilon) = \lambda \mathrm{tr}(\epsilon) + 2\mu\epsilon$ and $e^{ij} = e_i \otimes e_j$. We thus need to solve three problems: $V^{xx}, V^{yy}$ and $V^{xy}$ (avec $V^{yx} = V^{xy}$). The fields $V^{ij}$ are defined up to constant which can be fixed to

$$ \lim_{z_m \to \pm\infty} V^{ij} = \pm \frac{b^{ij}}{2}. \qquad (37) $$

We also introduce the vectors:

$$ c^{ij} = \int_{\mathcal{Y}_p} A(\epsilon(V^{ij})) e_x dx_m. \qquad (38) $$

Due to the symmetry of the cavity geometry, the following relations hold: $b_x^{xx} = b_z^{xz} = b_x^{zz} = 0$ and we define

$$ B_z = b_z^{zz} + a/h, \quad B_x = b_x^{xx}, \quad B = b_x^{xz} + a/h. \qquad (39) $$

We also have $c_z^{zz} = c_z^{xx} = 0$, $c_z^{xz} = \mu\varphi a/h$ and for symmetric cavities, $c_x^{xz} = 0$, and we define

$$ C_z = \frac{c_x^{zz}}{\mu}, \quad C_x = \frac{c_x^{xx}}{\mu}, \qquad (40) $$

These coefficients satisfy the relation (corresponding to Eq. (3.23) in [22]):

$$ \lambda B_z - \mu C_z = (\lambda + 2\mu)B_x + \lambda(1 - \varphi)a/h. \qquad (41) $$

By applying this effective behavior to the asymptotic analysis of Eq. (2.2) in [22] (after accounting for changes in notation), we recover the jump conditions given in Eq. (4).

## Appendix B. The near-clogging and clogging cases

The near-clogging regime corresponds to $\varphi \to 1$, with $\varphi = 1$ representing the limiting configuration of a homogeneous air layer of thickness $a$ that fully separates the top layer from the underlying elastic half-space. A natural question is whether our effective jump conditions (4) correctly recover this limiting case, for which the boundary conditions must be stress-free at both interfaces $z = -d$ and $z = -d_a$.

To address this point, Fig. 13 shows the variations of the effective coefficients ($B, B_z, B_x, C_z, C_x$) with respect to the filling ratio $\varphi$. As already visible in Table 3, the parameters ($B, B_z, B_x$) increase with $\varphi$, whereas ($C_z, C_x$) tend to zero. More specifically, we observe that in the limit $\varphi \to 1$, the effective parameters are well described by the following power-law fits $C_x = c_x(1-\varphi)$, $C_z = c_z(1-\varphi)$, $B_x = b_x(1-\varphi)^{-1}$, and $B = b(1-\varphi)^{-2.5}$, where ($c_x, c_z, b_x, b$) are independent of $\varphi$ (they depend on $a$; in particular $c_x = c_z = 0$ for $a = 0$, i.e., for cracks, while $b_x$ and $b$ remain nonzero). These trends are consistent with the expected transition, as $\varphi$ approaches 1, from jump-type interface conditions to stress-free boundary conditions.
To make this transition explicit, we reorganize the set of jump conditions (4) as

$$ \langle \sigma_{xz} \rangle = \frac{\lambda}{hB} \left( [u_x] + a \frac{\partial \langle u_z \rangle}{\partial x} \right), \qquad \langle \sigma_{zz} \rangle = \frac{\lambda}{hB_x} \left( [u_z] - \left( \frac{h\mu}{\lambda} C_z + (1-\varphi)a \right) \left( \frac{\partial u_z}{\partial z} \right) \right), \\ [\sigma_{xz}] = (1-\varphi)a \left( \frac{\partial \sigma_{xz}}{\partial z} \right) - \mu h \frac{\partial}{\partial x} \left( C_x \frac{\partial \langle u_x \rangle}{\partial x} + C_z \left( \frac{\partial u_z}{\partial z} \right) \right), \qquad [\sigma_{zz}] = (1-\varphi)a \left( \frac{\partial \sigma_{zz}}{\partial z} \right) - \varphi a \frac{\partial \langle \sigma_{xz} \rangle}{\partial x}. \qquad (42) $$

(using relation (41) to eliminate $B_z$).

In the limit $\varphi \to 1$, these expressions reduce to

$$ \langle \sigma_{xz} \rangle = (1-\varphi)^{2.5} \frac{\lambda}{hb} \left( [u_x] + a \frac{\partial \langle u_z \rangle}{\partial x} \right), \qquad [\sigma_{xz}] = (1-\varphi) \left[ a \left( \frac{\partial \sigma_{xz}}{\partial z} \right) - \mu h \frac{\partial}{\partial x} \left( c_x \frac{\partial \langle u_x \rangle}{\partial x} + c_z \left( \frac{\partial u_z}{\partial z} \right) \right) \right], \qquad (43) $$

which shows that, as $\varphi \to 1$, $\sigma_{xz}$ vanishes at both interfaces $z = -d$ and $z = -d_u$, with no singular behavior at $\varphi = 1$. Similarly,

$$ \langle \sigma_{zz} \rangle = \frac{(1-\varphi)}{hb_x} \left( [u_z] - (1-\varphi) \left( \frac{h\mu}{\lambda} c_z + a \right) \left( \frac{\partial u_z}{\partial z} \right) \right), \qquad [\sigma_{zz}] = (1-\varphi)a \left( \frac{\partial \sigma_{zz}}{\partial z} \right) - \varphi a \frac{\partial \langle \sigma_{xz} \rangle}{\partial x}. \qquad (44) $$

Since $\langle \sigma_{xz} \rangle = 0$ for all $x$, we also have $\frac{\partial \langle \sigma_{xz} \rangle}{\partial x} = 0$, and hence $\sigma_{zz} = 0$ at both $z = -d$ and $z = -d_u$.

Consequently, in the near-clogging limit, the jump conditions continuously evolve toward the physically expected stress-free boundary conditions at both interfaces, confirming the internal consistency of the effective model. A similar behavior is found in acoustics, where – mutatis mutandis (the cavities are sound-hard inclusions and the elastic medium is replaced by a fluid) – the interface is characterized by a single “blockage coefficient” $B$. There too, one observes a smooth transition from continuity of pressure and normal velocity for $\varphi \to 0$ to vanishing normal velocity for $\varphi \to 1$, without any singular limit [29].

## Data availability

The results presented in this article are generated by numerical codes developed by the authors, as described in Section 2.3 and in [24]. External datasets were not used. The numerical data supporting the findings of this study are available from the corresponding author upon reasonable request.

## References

[1] J. Jones, J. Whittier, Waves at a flexibly bonded interface, J. Appl. Mech. 34 (4) (1967) 905-909.

[2] G.S. Murty, Reflection, transmission and attenuation of elastic waves at a loosely-bonded interface of two half spaces, Geophys. J. Int. 44 (2) (1976) 389-404.

[3] J.-M. Baik, R.B. Thompson, Ultrasonic scattering from imperfect interfaces: a quasi-static model, J. Nondestruct. Eval. 4 (3-4) (1984) 177-196.

[4] P. Martin, Boundary integral equations for the scattering of elastic waves by elastic inclusions with thin interface layers, J. Nondestruct. Eval. 11 (3-4) (1992) 167-174.

[5] Y. Angel, J.D. Achenbach, Reflection and transmission of elastic waves by a periodic array of cracks, J. Appl. Mech. 52 (1) (1985) 33-41.

[6] Y. Angel, J. Achenbach, Reflection of ultrasonic waves by an array of microcracks, in: Review of Progress in Quantitative Nondestructive Evaluation, Springer, 1985, pp. 83-89.

[7] Y. Mikata, J. Achenbach, Interaction of harmonic waves with a periodic array of inclined cracks, Wave Motion 10 (1) (1988) 59-72.

[8] J.D. Achenbach, M. Kitahara, Reflection and transmission of an obliquely incident wave by an array of spherical cavities, J. Acoust. Soc. Am. 80 (4) (1986) 1209-1214.

[9] J.D. Achenbach, M. Kitahara, Y. Mikata, D. Sotiropoulos, Reflection and transmission of plane waves by a layer of compact inhomogeneities, Pure Appl. Geophys. 128 (1-2) (1988) 101-118.

[10] F.T. Gizzi, N. Masini, Historical damage pattern and differential seismic effects in a town with ground cavities: a case study from Southern Italy, Eng. Geol. 88 (1-2) (2006) 41-58.

[11] F. Gizzi, A. Loperte, A. Satriani, V. Lapenna, N. Masini, M. Proto, Georadar investigations to detect cavities in a historical town damaged by an earthquake of the past, Adv. Geosci. 24 (2010) 15-21.

[12] C. Smerzini, J. Aviles, R. Paolucci, F. Sánchez-Sesma, Effect of underground cavities on surface earthquake ground motion under SH wave propagation, Earthq. Eng. Struct. Dyn. 38 (12) (2009) 1441-1460.

[13] H. Alielahi, M. Kamalian, M. Adampira, Seismic ground amplification by unlined tunnels subjected to vertically propagating SV and P waves using BEM, Soil Dyn. Earthq. Eng. 71 (2015) 63-79.

[14] H. Alielahi, M. Adampira, Site-specific response spectra for seismic motions in half-plane with shallow cavities, Soil Dyn. Earthq. Eng. 80 (2016) 163-167.

[15] L. Huang, Z. Liu, C. Wu, J. Liang, The scattering of plane p, SV waves by twin lining tunnels with imperfect interfaces embedded in an elastic half-space, Tunn. Undergr. Space Technol. 85 (2019) 319-330.

[16] Z. Liu, J. Liu, S. Meng, X. Sun, Diffraction of elastic waves by a fluid-filled crack in a fluid-saturated poroelastic half-space, Geophys. J. Int. 225 (3) (2021) 1530-1553.

[17] D. Feizi, J. Asgari Marnani, H. Alielahi, M. Panji, Seismic ground amplification induced by box-shaped tunnels, Earthq. Eng. Eng. Vib. 21 (3) (2022) 697-714.

[18] S. Brûlé, E. Javelaud, S. Enoch, S. Guenneau, Experiments on seismic metamaterials: molding surface waves, Phys. Rev. Lett. 112 (13) (2014) 133901.

[19] A. Colombi, P. Roux, S. Guenneau, P. Gueguen, R.V. Craster, Forests as natural seismic metamaterials, Sci. Rep. 6 (2016) 19238.

[20] S. Brûlé, S. Enoch, S. Guenneau, Emergence of seismic metamaterials: Current state and future perspectives, Phys. Lett. A 384 (1) (2020) 126034.

[21] M. Miniaci, A. Krushynska, F. Bosia, N.M. Pugno, Metabarriers for seismic protection based on locally resonant metamaterials, Phys. Rev. B 104 (9) (2021) 094305.

[22] K. Pham, A. Maurel, J.J. Marigo, Revisiting imperfect interface laws for two-dimensional elastodynamics, Proc. R. Soc. A 477 (2245) (2021) 20200519.

[23] J. Achenbach, Wave Propagation in Elastic Solids, Elsevier, 2012.

[24] A. Maurel, K. Pham, Multimodal method for the scattering by an array of plates connected to an elastic half-space, J. Acoust. Soc. Am. 146 (6) (2019) 4402-4412.
[25] J.J. Marigo, K. Pham, A. Maurel, S. Guenneau, Effective model for elastic waves in a substrate supporting an array of plates/beams with flexural and longitudinal resonances, J. Elasticity 146 (1) (2021) 143–177.

[26] J.J. Marigo, K. Pham, A. Maurel, S. Guenneau, Effective model for elastic waves propagating in a substrate supporting a dense array of plates/beams with flexural resonances, J. Mech. Phys. Solids 143 (2020) 104029.

[27] K. Pham, A. Maurel, S. Félix, S. Guenneau, Hybridized love waves in a guiding layer supporting an array of plates with decorative endings, Materials 13 (7) (2020) 1632.

[28] S. Enoch, S. Brûlé, E. Javelaud, S. Guenneau, Flat lens effect on seismic waves propagation in the subsoil, Sci. Rep. 7 (1) (2017) 18066.

[29] J.J. Marigo, A. Maurel, Homogenization models for thin rigid structured surfaces and films, J. Acoust. Soc. Am. 140 (1) (2016) 260–273.