1. Introduction
Graded metamaterials, i.e. media composed of a spatially varying resonant substructure, have been shown to exhibit remarkable properties across a wide range of wave sciences, including electromagnetic waves (Tsakmakidis, Boardman & Hess Reference Tsakmakidis, Boardman and Hess2007), acoustic waves (Zhu et al. Reference Zhu, Chen, Zhu, Garcia-Vidal, Yin, Zhang and Zhang2013; Ni et al. Reference Ni, Wu, Chen, Zheng, Xu, Nayar, Liu, Lu and Chen2014; Jiménez et al. Reference Jiménez, Romero-García, Pagneux and Groby2017; Bennetts, Peter & Craster Reference Bennetts, Peter and Craster2019; Zhao & Zhou Reference Zhao and Zhou2019), elastic waves (Skelton et al. Reference Skelton, Craster, Colombi and Colquitt2018; Arreola-Lucas et al. Reference Arreola-Lucas, Báez, Cervera, Climente, Méndez-Sánchez and Sánchez-Dehesa2019), seismic waves (Colombi et al. Reference Colombi, Colquitt, Roux, Guenneau and Craster2016) and water waves; the latter area is the focus of this paper. Modelling has shown that water wave metamaterials can locally amplify wave energy; the subsequent absorption or reflection of this energy may have important applications in the design of wave energy converters and coastal breakwaters, respectively (Bennetts, Peter & Craster Reference Bennetts, Peter and Craster2018; Wilks, Montiel & Wakes Reference Wilks, Montiel and Wakes2022). Grading the geometric or physical properties of the resonators in the metamaterial can result in broadband effects. Such effects were also demonstrated experimentally in the context of water waves by Archer et al. (Reference Archer, Wolgamot, Orszaghova, Bennetts, Peter and Craster2020).
The most commonly used conceptual tool for understanding wave propagation through graded metamaterials has been what we shall refer to as the local Bloch wave approximation (LBWA). This approach assumes that if the grading of the metamaterial is sufficiently weak, then the wave field inside the metamaterial can locally be approximated as a sum of propagating solutions of the Bloch–Floquet problem, in which the unit cell of the infinite array matches the local geometry of the metamaterial. These propagating solutions, which are associated with a frequency-dependent group velocity, only exist in frequency intervals known as passbands. In the complementary frequency intervals (known as stopbands), waves cannot propagate, i.e. only decaying modes are present. The LBWA explains the properties of graded metamaterials as follows.
(i) Waves of a given frequency enter the metamaterial and gradually slow down because the group velocity of the local Bloch wave is decreasing;
(ii) When the group velocity becomes zero at the so-called turning point, propagation farther into the metamaterial is forbidden because the local geometry only supports decaying modes;
(iii) Since propagation beyond the turning point is impossible, the local forward-propagating wave must couple with the local backward-propagating wave at the turning point, which continues to propagate backwards through the array, ultimately leading to reflection;
(iv) The near-zero group velocity of the forward-propagating and backward-propagating Bloch wave near the turning point results in local energy amplification;
(v) The grading causes the location of the turning point to be frequency dependent, therefore allowing the metamaterial to induce the spatial separation of broadband signals into their frequency components.
Variations of the above argument can be found in a large number of research papers on the topic of graded metamaterials (Tsakmakidis et al. Reference Tsakmakidis, Boardman and Hess2007; Zhu et al. Reference Zhu, Chen, Zhu, Garcia-Vidal, Yin, Zhang and Zhang2013; Cebrecos et al. Reference Cebrecos, Picó, Sánchez-Morcillo, Staliunas, Romero-García and Garcia-Raffi2014; Bennetts et al. Reference Bennetts, Peter and Craster2018; Wilks et al. Reference Wilks, Montiel and Wakes2022). The interaction of waves with graded metamaterials is often termed rainbow trapping, due to the ability of the metamaterial to separate frequencies. As noted by Chaplain et al. (Reference Chaplain, Pajer, De Ponti and Craster2020), a better term for this interaction is rainbow reflection, since local amplifications of energy gradually decay, while trapped modes do not (also see He, He & Jin (Reference He, He and Jin2012)).
While there are many examples of the qualitative version of the LBWA being successfully used to explain wave interaction with graded metamaterials, few studies have examined the approximation quantitatively. In particular, it is not known how the accuracy of the approximation depends on the grading strength of the metamaterial, i.e. the rate at which the resonant properties of the metamaterial vary spatially. Some discussions based on the LBWA assume that small reflections of the local Bloch wave – which can occur between adjacent unit cells with different properties – may be ignored in weakly graded metamaterials. In other words, it is assumed that the transmission of Bloch waves is perfect except at the turning point, where instead perfect reflection occurs. Tsakmakidis et al. (Reference Tsakmakidis, Boardman and Hess2007), who first introduced the concept of rainbow reflection in their study of an adiabatically tapered electromagnetic waveguide, found that an excellent approximation of the wave-power transport into the device was obtained by assuming that these reflections were negligible. However, the extent to which these reflections remain negligible in other metamaterials with more significant grading is unclear.
Another key assumption of the LBWA is that decaying modes – i.e. non-propagating solutions to the Bloch–Floquet problem which are excited in graded arrays – may be ignored. Cebrecos et al. (Reference Cebrecos, Picó, Sánchez-Morcillo, Staliunas, Romero-García and Garcia-Raffi2014), who studied chirped sonic crystals consisting of cylinders, used coupled mode theory to approximate the array as a multilayered structure which only supports wave propagation along the direction of wave incidence. By further assuming that the dispersion relation of the crystal varies smoothly with the spatial parameter, the authors derived a system of ordinary differential equations that accurately predicts the wave amplitude envelope in the crystal. In their study of water wave interaction with line arrays of split-ring resonators with graded radii, Bennetts et al. (Reference Bennetts, Peter and Craster2018) showed that the solution was locally dominated by Rayleigh–Bloch waves, i.e. propagating homogeneous solutions to the problem of diffraction by an infinite periodic line array of scatterers (for details see Porter & Evans (Reference Porter and Evans1999)). This was done by comparing the full numerical solution with an approximation, in which the velocity potential in each unit cell is represented using solely the Rayleigh–Bloch eigenfunctions of the unit cell's transfer matrix. However, none of these papers have studied the LBWA mechanistically. Here, we will investigate wave propagation in a graded metamaterial by exploring the coupling mechanism between Bloch waves in adjacent unit cells. To do this, we consider the scattering of Bloch waves across the interface between two semi-infinite uniform arrays positioned end-to-end.
The idea of using the semi-infinite problem to understand complex wave phenomena in finite arrays is not new. The present paper is inspired by the work of Thompson, Linton & Porter (Reference Thompson, Linton and Porter2008) who studied the interaction of water waves by long finite line arrays of cylinders. These authors sought to understand the resonances of these arrays further, which had previously been discussed by Maniar & Newman (Reference Maniar and Newman1997) and are related to Rayleigh–Bloch waves. Thompson et al. (Reference Thompson, Linton and Porter2008) decomposed the wave scattered by each cylinder into (i) the diffracted wave assuming the cylinder belongs to an infinite line array, (ii) the left and right-propagating Rayleigh–Bloch waves and (iii) decaying circular waves which originate at the array ends. The amplitudes of the Rayleigh–Bloch waves, which are excited by the incident waves at the array ends, were computed from the solution to the semi-infinite problem derived by Linton, Porter & Thompson (Reference Linton, Porter and Thompson2007). The so-called large array approximation of Thompson et al. (Reference Thompson, Linton and Porter2008), which assumes that the interactions of the decaying circular waves at the opposite array ends can be ignored, yielded numerically accurate predictions of the scattered wave coefficients. We will apply a similar modelling approach to study graded arrays of vertical barriers in a two-dimensional fluid.
In our previous work (Wilks et al. Reference Wilks, Montiel and Wakes2022), we demonstrated rainbow reflection in graded arrays of surface-piercing vertical barriers of point-thickness in a two-dimensional fluid of constant and finite depth. We claimed the LBWA applied to this problem and used it to explain the rainbow reflection, although we did not verify this claim. Subsequently, we augmented the boundary value problem for the graded array of vertical barriers to include an energy absorption mechanism (i.e. damped pistons) in the interbarrier regions. By optimising the parameters of the device, we showed that it could achieve near-perfect energy absorption over a prescribed frequency interval. More recently, Huang & Porter (Reference Huang and Porter2023) considered the related problem of a water wave metamaterial consisting of an equally spaced array of identical vertical barriers. The authors solved the problem (i) numerically using an efficient transfer-matrix method and (ii) via homogenisation in the small spacing limit. The agreement between these methods is excellent if this small-spacing assumption is valid, although the homogenised model fails as the resonant frequency of the array is approached. Here, we will also consider the scattering of water waves by arrays of vertical barriers. This setting is advantageous because the scattering problem is simple to formulate and solve using computationally efficient methods, yet it is sufficiently non-trivial that metamaterials can support both propagating and decaying solutions to the Bloch–Floquet problem at the same frequency.
In § 2, we solve the scattering of Bloch waves across the interface between two semi-infinite arrays of vertical barriers. We arrive at this solution hierarchically, after discussing (i) scattering by a single vertical barrier, (ii) scattering by finite arrays of vertical barriers, (iii) propagation of Bloch waves through infinite periodic arrays of vertical barriers and (iv) scattering by semi-infinite arrays of vertical barriers. In § 3, we present results for the transmission of Bloch waves across the interface between two semi-infinite arrays of vertical barriers. In § 4, we numerically implement the LBWA based on the assumption that the coupling of Bloch waves between adjacent subregions in a graded array behaves like the coupling across the interface between two semi-infinite arrays. Our results suggest that the LBWA accurately predicts the free surface amplitude at most frequencies, with the exception of intervals which lie just above the cutoff frequencies associated with each of the vertical barriers in the array. Specifically, these cutoff frequencies are the maximum frequency at which the corresponding infinite array problem supports propagating Bloch solutions. We argue that these errors are due to the excitation of slowly attenuating Bloch waves in the region immediately beyond the turning point.
2. Scattering by arrays of vertical barriers
2.1. A single vertical barrier
We begin by considering the problem of wave scattering by a single surface-piercing vertical barrier, which is a classic problem in linear water wave theory (Ursell Reference Ursell1947; John Reference John1948). Here, we consider the finite-bathymetry version of the problem. In our notation, the fluid occupies the region 
$\varOmega =\{(x,z)|x\in \mathbb {R},-H< z<0\}\setminus \varGamma$, where 
$H$ is the depth, the sea-bed being situated at 
$z=-H$. The mean position of the free surface is situated at 
$z=0$ and 
$\varGamma =\{(0,z)|-d< z<0\}$ describes the barrier, which has submergence depth of 
$d$. Under the usual assumptions of time-harmonic linear water wave theory (Linton & McIver Reference Linton and McIver2001; Mei, Stiassnie & Yue Reference Mei, Stiassnie and Yue2005), the velocity potential of the fluid is described by 
$\mathrm {Re}(\phi (x,z)\, {\rm e}^{-\mathrm {i}\omega t})$, where 
$t$ is time, 
$\omega$ is the angular frequency, and the complex valued function 
$\phi$ satisfies the following boundary value problem:
\begin{gather} \bigtriangleup \phi =0 \quad (x,z)\in\varOmega, \end{gather}
\begin{gather}\partial_z\phi=\frac{\omega^2}{g}\phi \quad z=0, \end{gather}
\begin{gather}\partial_z \phi=0 \quad z={-}H, \end{gather}
\begin{gather}\left(\frac{\partial}{\partial x} \mp \mathrm{i} k_0\right)(\phi-\phi^{In})\to 0 \quad \mbox{as\ }k_0x\to\pm\infty, \end{gather}
\begin{gather}\partial_x\phi=0 \quad (x,z)\in\varGamma, \end{gather}
\begin{gather}\sqrt{x^2+(z+d)^2}\|\nabla \phi\|\to 0 \quad \mbox{as\ }\sqrt{x^2+(z+d)^2}\to 0, \end{gather}
where 
$g$ is acceleration due to gravity and 
$\phi ^{In}$ is the potential of the prescribed incident wave. Moreover, the quantities 
$k_m$ are the solutions to the dispersion relation 
$k\tanh kH=\omega ^2/g$, with 
$k_0\in \mathbb {R}^+$ and 
${-\mathrm {i} k_m\in ((m-1){\rm \pi} /H,m{\rm \pi} /H)}$ for all 
${m\in \mathbb {N}}$. After applying separation of variables, we adopt a piecewise ansatz for the solution to (2.1) of the form
\begin{equation} \phi(x,z)=\begin{cases} \displaystyle \sum_{m=0}^\infty (A_m^{(0)}\exp(\mathrm{i} k_m x)+B_m^{(0)}\exp(-\mathrm{i} k_m x))\psi_m(z), & x<0,\\ \displaystyle \sum_{m=0}^\infty (A_m^{(1)}\exp(\mathrm{i} k_m x)+B_m^{(1)}\exp(-\mathrm{i} k_m x))\psi_m(z), & x>0, \end{cases} \end{equation}where the vertical eigenfunctions have been defined as
\begin{equation} \psi_m(z) = \left(\frac{\sinh(2k_m H)}{4k_m H}+\frac{1}{2}\right)^{{-}1/2}\cosh(k_m(z+H)). \end{equation}
The expansions in (2.2) automatically satisfy (2.1a–c). The Sommerfeld radiation condition (2.1d) allows us to formulate the boundary value problem as a scattering problem, in which the unknown coefficient vectors of the scattered field (
$\boldsymbol {A}^{(1)}$ and 
$\boldsymbol {B}^{(0)}$) are sought in terms of the prescribed coefficient vectors of the incident field (
$\boldsymbol {A}^{(0)}$ and 
$\boldsymbol {B}^{(1)}$). The entries of these coefficient vectors are the coefficients in (2.2) after the infinite sums have been truncated at 
$m=N_{sol}$. The coefficient vectors of the scattered field must be chosen so that 
$\phi (x,z)$ satisfies the boundary condition on the barrier (2.1e) and is continuously differentiable across the gap beneath the barriers. This is accomplished using the integral equation/Galerkin method of Porter & Evans (Reference Porter and Evans1995), which leverages prior knowledge of the singularity at the tip of the barrier (2.1f) to enable accurate computation of the coefficient vectors. The method yields a scattering matrix (which we express in blockwise form) that satisfies
\begin{equation} \begin{bmatrix} {\mathsf{S}}_{11} & {\mathsf{S}}_{12}\\ {\mathsf{S}}_{21} & {\mathsf{S}}_{22} \end{bmatrix} \begin{bmatrix} \boldsymbol{A}^{(0)}\\ \boldsymbol{B}^{(1)} \end{bmatrix} = \begin{bmatrix} \boldsymbol{A}^{(1)}\\ \boldsymbol{B}^{(0)} \end{bmatrix}. \end{equation}
Other than 
$N_{sol}$, the accuracy of the Galerkin solution depends on two other truncation parameters, which we have tuned so that at least five-figure accuracy of 
$B_0^{(0)}$ and 
$A_0^{(1)}$ is achieved across a range of test frequencies. The reader is referred to Porter & Evans (Reference Porter and Evans1995) for further details.
2.2. Finite arrays of vertical barriers
Next, we consider the scattering by 
$N+1$ vertical barriers which are positioned at 
$x=nW$ for 
$n\in \{0,\ldots,N\}$. Although the treatment here is essentially identical to that in our earlier paper (Wilks et al. Reference Wilks, Montiel and Wakes2022), we use this subsection to introduce terms which will be used subsequently in this paper. We restrict the problem to the case where the horizontal distance between barriers 
$W$ and their submergence depth 
$d$ are constant. Both restrictions could be easily relaxed (the latter by incorporating different scattering matrices at each of the different barriers) but this is not done here for ease of exposition. The solution to the multiple scattering problem may be written in the form
\begin{equation} \phi(x,z)=\sum_{m=0}^\infty (A_m^{(n)} \exp(\mathrm{i} k_m (x-nW))+B_m^{(n)}\exp(-\mathrm{i} k_m (x-nW)))\psi_m(z), \end{equation}
for 
$n\in \{0,\ldots,N+1\}$ and
\begin{equation} x\in I_n := \begin{cases} (-\infty,0), & n=0,\\ ((n-1)W,nW), & 0< n\leq N,\\ (NW,\infty), & n=N+1. \end{cases} \end{equation}
A scattering matrix equation analogous to (2.4), modified to incorporate the phase-shift matrix 
$\boldsymbol{\mathsf{L}}=\mathrm {diag}(\exp (\mathrm {i} k_m W))$, in which 
$0\leq m\leq N_{sol}$, can be written as
\begin{equation} \begin{bmatrix} \boldsymbol{\mathsf{L}} {\mathsf{S}}_{11} & \boldsymbol{\mathsf{L}} {\mathsf{S}}_{12} \boldsymbol{\mathsf{L}}\\ {\mathsf{S}}_{21} & {\mathsf{S}}_{22}\boldsymbol{\mathsf{L}} \end{bmatrix} \begin{bmatrix} \boldsymbol{A}^{(n-1)}\\ \boldsymbol{B}^{(n)} \end{bmatrix} = \begin{bmatrix} \boldsymbol{A}^{(n)}\\ \boldsymbol{B}^{(n-1)} \end{bmatrix}, \end{equation}
for all 
$n\in \{1,\ldots,N+1\}$. In the multiple scattering problem, 
$\boldsymbol {A}^{(0)}$ and 
$\boldsymbol {B}^{(N+1)}$ are known. To compute the remaining unknown coefficients, we use the scattering matrix method (Ko & Sambles Reference Ko and Sambles1988; Bennetts & Squire Reference Bennetts and Squire2009; Montiel, Squire & Bennetts Reference Montiel, Squire and Bennetts2015) to obtain a scattering matrix for the whole array and compute the unknown coefficients recursively. This method remains numerically stable for large 
$N$.
2.3. Infinite arrays of vertical barriers
We now seek propagating solutions to the infinite array problem, i.e. the case where the barriers are positioned at 
$nW$ for all 
$n\in \mathbb {Z}$. Despite our treatment here being again essentially identical to that in our previous paper (Wilks et al. Reference Wilks, Montiel and Wakes2022), we introduce terms that will be important in what follows. Our formulation and method are also analogous to those of Peter & Meylan (Reference Peter and Meylan2009). In the infinite array problem, we require that (2.5) and (2.7) hold for all 
$n\in \mathbb {Z}$ after redefining the horizontal intervals as 
$I_n=((n-1)W,nW)$. Bloch's theorem motivates seeking quasiperiodic solutions of the form
\begin{equation} \phi(x+nW,z)=\exp({\pm}\mathrm{i} q n W)\phi(x,z), \end{equation}
for all 
$n\in \mathbb {Z}$, where 
$q (-q)$ is the unknown Bloch wavenumber of the forward-propagating (backward-propagating) Bloch wave. Periodicity arguments allow us to restrict 
$q\in (0,{\rm \pi} /W)$. Solutions to (2.8) may easily be shown to satisfy
\begin{gather} A_m^{(n)}=\exp({\pm}\mathrm{i} q n W)A_m^{{\pm}}, \end{gather}
\begin{gather}B_m^{(n)}=\exp({\pm}\mathrm{i} q n W)B_m^{{\pm}}, \end{gather}
i.e. the coefficients in different interbarrier regions are related via a phase-shift which depends on the direction of the Bloch wave. Combining equations (2.9) and (2.7) for 
$n=1$ yields the following generalised eigenvalue equation:
\begin{equation} \begin{bmatrix} \boldsymbol{\mathsf{L}} {\mathsf{S}}_{11} & 0_{M}\\{\mathsf{S}}_{21} & -\boldsymbol{\mathsf{I}}\end{bmatrix} \begin{bmatrix} \boldsymbol{A}^{{\pm}}\\ \boldsymbol{B}^{{\pm}} \end{bmatrix} = \exp({\pm}\mathrm{i} q W) \begin{bmatrix}\boldsymbol{\mathsf{I}} & -\boldsymbol{\mathsf{L}} {\mathsf{S}}_{12}\boldsymbol{\mathsf{L}}\\0_{M} & -{\mathsf{S}}_{22}\boldsymbol{\mathsf{L}}\end{bmatrix} \begin{bmatrix} \boldsymbol{A}^{{\pm}}\\ \boldsymbol{B}^{{\pm}} \end{bmatrix}, \end{equation}
where 
$\boldsymbol{\mathsf{I}}$ and 
$0_{M}$ denote the identity and zero matrices of dimension 
$N_{sol}+1$, respectively. If generalised eigenvalues of (2.10) exist (do not exist) on the unit circle at a given angular frequency 
$\omega$, then 
$\omega$ is said to lie in a passband (stopband). The eigenvectors of (2.10) are normalised so that
\begin{equation} |A_0^\pm|^2-|B_0^\pm|^2+2\mathrm{i}\sum_{m=1}^\infty \frac{k_m}{k_0}\mathrm{Im}(A_m^\pm \overline{B_m^\pm})=1, \end{equation}
in which 
$\bar {w}$ denotes the complex conjugate of 
$w$. This normalisation ensures that the conservation of energy identities given in §§ 2.4 and 2.5 take a familiar form.
2.4. Semi-infinite arrays of vertical barriers
Next, we consider the problem of scattering by a semi-infinite array of vertical barriers, which are positioned at 
$x=nW$ for 
$n\in \mathbb {N}_0$ (i.e. 
$n$ is a non-negative integer) and have identical submergence depth 
$d$. To the best of our knowledge, this particular problem has not been considered before, although the scattering of water waves by semi-infinite arrays has previously been considered in both two and three dimensions for different scatterer geometries (Porter & Evans Reference Porter and Evans2006; Linton et al. Reference Linton, Porter and Thompson2007; Peter & Meylan Reference Peter and Meylan2007). Our treatment here, which uses the so-called filtering method when Bloch waves are supported by the semi-infinite array, is based on these papers. For other solution methods to the semi-infinite diffraction problem, see Martin, Abrahams & Parnell (Reference Martin, Abrahams and Parnell2015) and Joseph & Craster (Reference Joseph and Craster2015).
The separation of variables expansion for the problem is given by (2.5), in which 
$n\in \mathbb {N}_0$ and
\begin{equation} I_n=\begin{cases} (-\infty,0), & n=0,\\ ((n-1)W,nW), & n>0. \end{cases} \end{equation}
Equation (2.7) then holds for all 
$n\in \mathbb {N}$. We first consider the case where 
$\omega$ is in a passband of the semi-infinite array. Anticipating the use of the solution to solve the coupling between two semi-infinite arrays, we permit the incident wave to have both left-travelling and right-travelling components. The right-travelling component is described generally by the known coefficient vector 
$\boldsymbol {A}^{(0)}$, whereas the left-travelling component is restricted to consist only of a leftward-propagating Bloch wave of known complex amplitude 
$\beta$. The scattered wave in 
$x<0$ is described by the coefficient vector 
$\boldsymbol {B}^{(0)}$. In the semi-infinite array, the scattered field consists of a rightward-propagating Bloch wave of unknown complex amplitude 
$\alpha$, as well as a component made up of non-propagating Bloch waves that decay as 
$x\to +\infty$. This motivates decomposing the coefficient vectors as
\begin{gather} \boldsymbol{A}^{(n)} = \alpha \exp(\mathrm{i} q n W)\boldsymbol{A}^+ + \beta\exp(-\mathrm{i} q n W)\boldsymbol{A}^- + \boldsymbol{C}^{(n)}, \end{gather}
\begin{gather}\boldsymbol{B}^{(n)} = \alpha \exp(\mathrm{i} q n W)\boldsymbol{B}^+ + \beta\exp(-\mathrm{i} q n W)\boldsymbol{B}^- + \boldsymbol{D}^{(n)}, \end{gather}
for all 
$n\geq 1$, where 
$\boldsymbol {C}^{(n)}$ and 
$\boldsymbol {D}^{(n)}$ describe the component of the wave that decays as 
$x\to \infty$. This implies that 
$\boldsymbol {C}^{(n)},\boldsymbol {D}^{(n)}\to \boldsymbol {0}_{V}$ as 
$n\to \infty$, where 
$\boldsymbol {0}_{V}$ is the 
$(N_{sol}+1)$-dimensional zero column vector.
Next, we obtain the unknown coefficients using the filtering method, which leverages information about the phase of the scattered wave (implied by the Bloch wavenumber) to achieve a significant reduction of the numerical error of the amplitude of the scattered Bloch wave compared with more direct methods (Linton et al. Reference Linton, Porter and Thompson2007). To do this, we define
\begin{gather} \boldsymbol{\tilde{A}}^{(n)} := \begin{cases} \boldsymbol{A}^{(n)}-\beta\exp(-\mathrm{i} q n W)\boldsymbol{A}^- , & n>0,\\ \boldsymbol{A}^{(0)}, & n=0, \end{cases} \end{gather}
\begin{gather}\boldsymbol{\tilde{B}}^{(n)} := \begin{cases} \boldsymbol{B}^{(n)}-\beta\exp(-\mathrm{i} q n W)\boldsymbol{B}^- , & n>0,\\ \boldsymbol{B}^{(0)}, & n=0, \end{cases} \end{gather}
which approach the coefficients of the scattered Bloch wave as 
$n\to \infty$, i.e. 
$\boldsymbol {\tilde {A}}^{(n)}\to \alpha \exp (\mathrm {i} q n W)\boldsymbol {A}^+$ and 
$\boldsymbol {\tilde {B}}^{(n)}\to \alpha \exp (\mathrm {i} q n W)\boldsymbol {B}^+$. We also define
\begin{gather} \boldsymbol{\tilde{C}}^{(n)} := \begin{cases} \boldsymbol{\tilde{A}}^{(n)}-\exp(\mathrm{i} q W)\boldsymbol{\tilde{A}}^{(n-1)} , & n>0,\\ \boldsymbol{A}^{(0)} , & n=0, \end{cases} \end{gather}
\begin{gather}\boldsymbol{\tilde{D}}^{(n)} := \begin{cases} \boldsymbol{\tilde{B}}^{(n)}-\exp(\mathrm{i} q W)\boldsymbol{\tilde{B}}^{(n-1)}, & n>0,\\ \boldsymbol{B}^{(0)} , & n=0, \end{cases} \end{gather}
which both decay to 
$\boldsymbol {0}_{V}$ as 
$n\to \infty$. The quantities in (2.14) and (2.15) are related by the following telescoping sum identities:
\begin{gather} \boldsymbol{\tilde{A}}^{(n)}= \sum_{j=0}^n \exp(\mathrm{i} q (n-j)W)\boldsymbol{\tilde{C}}^{(j)}, \end{gather}
\begin{gather}\boldsymbol{\tilde{B}}^{(n)} = \sum_{j=0}^n \exp(\mathrm{i} q (n-j)W)\boldsymbol{\tilde{D}}^{(j)}. \end{gather}
Next, we substitute (2.16) into (2.7), in order to reformulate the problem in terms of the unknowns 
$\boldsymbol {\tilde {C}}^{(n)}$ and 
$\boldsymbol {\tilde {D}}^{(n)}$. After some algebra, we eventually obtain
\begin{gather} \boldsymbol{\mathsf{L}}{\mathsf{S}}_{11}\boldsymbol{\tilde{C}}^{(0)}+\sum_{j=0}^1 \exp(\mathrm{i} q (1-j)W)[\boldsymbol{\mathsf{L}}{\mathsf{S}}_{12}\boldsymbol{\mathsf{L}}\boldsymbol{\tilde{D}}^{(j)}- \boldsymbol{\tilde{C}}^{(j)}]=\beta\exp(-\mathrm{i} q W)[\boldsymbol{A}^- -\boldsymbol{\mathsf{L}}{\mathsf{S}}_{12}\boldsymbol{\mathsf{L}}\boldsymbol{B}^-], \end{gather}
\begin{gather}{\mathsf{S}}_{21}\boldsymbol{\tilde{C}}^{(0)}+\exp(\mathrm{i} q W)[{\mathsf{S}}_{22}\boldsymbol{\mathsf{L}}-\boldsymbol{\mathsf{I}}]\boldsymbol{\tilde{D}}^{(0)}+{\mathsf{S}}_{22}\boldsymbol{\mathsf{L}} \boldsymbol{\tilde{D}}^{(1)}={-}\beta\exp(-\mathrm{i} q W){\mathsf{S}}_{22}\boldsymbol{\mathsf{L}}\boldsymbol{B}^-, \end{gather}
in the 
$n=1$ case, and
\begin{gather} \sum_{j=0}^{n-1}\exp(\mathrm{i} q(n-j)W)[(\exp(-\mathrm{i} q W)\boldsymbol{\mathsf{L}}{\mathsf{S}}_{11}-\boldsymbol{\mathsf{I}})\boldsymbol{\tilde{C}}^{(j)}+\boldsymbol{\mathsf{L}}{\mathsf{S}}_{12} \boldsymbol{\mathsf{L}}\boldsymbol{\tilde{D}}^{(j)}]-\boldsymbol{\tilde{C}}^{(n)} +\boldsymbol{\mathsf{L}}{\mathsf{S}}_{12}\boldsymbol{\mathsf{L}}\boldsymbol{\tilde{D}}^{(n)}=0, \end{gather}
\begin{gather}\sum_{j=0}^{n-1}\exp(\mathrm{i} q(n-j)W)[\exp(-\mathrm{i} q W){\mathsf{S}}_{21}\boldsymbol{\tilde{C}}^{(j)}+({\mathsf{S}}_{22}\boldsymbol{\mathsf{L}}-\exp(-\mathrm{i} q W))\boldsymbol{\tilde{D}}^{(j)}]+{\mathsf{S}}_{22} \boldsymbol{\mathsf{L}} \boldsymbol{\tilde{D}}^{(n)}=0, \end{gather}
in the 
$n>1$ case. In particular, (2.10) was used to show that the right-hand sides of (2.17c,d) are zero, which implies that the incident Bloch wave can only excite other wave modes at 
$x=0$.
To obtain the numerical solution, we assume that for some sufficiently large 
$P$, 
$\boldsymbol {\tilde {C}}^{(p)}=\boldsymbol {\tilde {D}}^{(p)}= \boldsymbol {0}_{V}$ for all 
$p>P$. A 
$2P(N_{sol}+1)$-dimensional linear system for the remaining unknown coefficients is assembled using (2.17a,b), (2.17c) for 
$n\in \{2,\ldots,P-1\}$, (2.17d) for 
$n\in \{2,\ldots,P\}$ and the statement 
$\boldsymbol {\tilde {C}}^{(0)}=\boldsymbol {A}^{(0)}$. After solving this linear system, the original desired coefficients 
$\boldsymbol {A}^{(n)}$ and 
$\boldsymbol {B}^{(n)}$ can then be obtained from (2.13), (2.14) and (2.16). Moreover, some matrix algebra can be used to obtain the transmission and reflection operators 
${\mathsf{T}}_f^+$, 
${\mathsf{T}}_b^+$, 
${\mathsf{R}}_f^+$ and 
${\mathsf{R}}_b^+$, which satisfy
\begin{gather} {\mathsf{T}}_f^+\boldsymbol{A}^{(0)}+{\mathsf{R}}_b^+\beta = \alpha , \end{gather}
\begin{gather}{\mathsf{R}}_f^+\boldsymbol{A}^{(0)}+{\mathsf{T}}_b^+\beta =\boldsymbol{B}^{(0)}, \end{gather}
in which the superscript 
$+$ indicates that the metamaterial occupies the positive half-line, while the subscripts 
$f$ and 
$b$ refer to the forward and backward scattering problems, respectively. In (2.18), 
${\mathsf{T}}_f^+$ is a 
$(N_{sol}+1)$-dimensional row vector, 
${\mathsf{R}}_b^+$ is a scalar, 
${\mathsf{R}}_f^+$ is a 
$(N_{sol}+1)$-dimensional square matrix and 
${\mathsf{T}}_b^+$ is a 
$(N_{sol}+1)$-dimensional column vector. If only propagating modes are incident to the system, (2.18) can be reduced to a 
$2\times 2$ matrix equation of the form
\begin{equation} \boldsymbol{\mathsf{S}}(0,d)\begin{bmatrix} A_0^{(0)}\\ \beta \end{bmatrix}= \begin{bmatrix} \alpha \\ B_0^{(0)} \end{bmatrix}, \end{equation}
i.e. it relates the amplitudes of the propagating waves. The reason for using the notation 
$\boldsymbol{\mathsf{S}}(0,d)$ for the 
$2\times 2$ matrix will become apparent in § 2.5.
Next, we consider conservation of energy. In the forward scattering problem, we set 
$A_m^{(0)}=\delta _{0m}$ (in which 
$\delta _{ij}$ denotes the Kronecker delta) and 
$\beta =0$. The reflection and transmission coefficients are then given by 
$R=B_0^{(0)}$ and 
$T=\alpha$. In the backward scattering problem, we set 
$\beta =1$ and 
$A_m^{(0)}=0$ for all 
$0< m< N_{sol}$. In this case, the reflection and transmission coefficients are given by 
$T=B_0^{(0)}$ and 
$R=\alpha$. In both scattering problems, the conservation of energy identity 
$|R|^2+|T|^2=1$ can be shown to hold by applying Green's second identity. The familiar form of this identity is guaranteed by the normalisation of the Bloch waves given in (2.11).
We remark that when 
$\omega$ is in a stopband, no propagating Bloch waves exist in the metamaterial region and the filtering method is not required. In this case, the only required diffraction matrix is 
${\mathsf{R}}_f^+$, which can be approximated accurately using the solution for a finite, regularly spaced array of vertical barriers, that is, the method outlined in § 2.2. In the absence of propagating Bloch waves, all of the matrices 
${\mathsf{T}}_f^+$, 
${\mathsf{T}}_b^+$ and 
${\mathsf{R}}_b^+$ are not well defined. For consistency of notation in this case, we write a 
$1\times 1$ scattering matrix equation of the form 
$\boldsymbol{\mathsf{S}}(0,d)A_0^{(0)}=B_0^{(0)}$. The conservation of energy identity for the forward scattering problem in which 
$A_m^{(0)}=\delta _{0m}$ is simply 
$|R|=1$, where 
$R= B_0^{(0)}$.
We also require the solution to the problem where the barriers occupy the negative half-line. Specifically, the barriers are situated at 
$x=-nW$ for 
$n\in \mathbb {N}$. This solution is obtained by applying the transformation 
$x\mapsto -(x+W)$. Following some lengthy algebra, we obtain the following relationships between the diffraction matrices for the positive and negative half-line problems:
\begin{gather} {\mathsf{R}}_b^- = \boldsymbol{\mathsf{L}}{\mathsf{R}}_f^+\boldsymbol{\mathsf{L}} \quad {\mathsf{T}}_b^- = J{\mathsf{T}}_f^+\boldsymbol{\mathsf{L}}, \end{gather}
\begin{gather}{\mathsf{T}}_f^+ = \boldsymbol{\mathsf{L}}{\mathsf{T}}_b^-J \quad {\mathsf{R}}_f^- = J^2{\mathsf{R}}_b^+, \end{gather}
where the number 
$J$ has been defined as
\begin{equation} J:=\exp(-\mathrm{i} k_m W)\frac{A_m^+}{B_m^-}=\exp(\mathrm{i} k_m W)\frac{B_m^+}{A_m^-}, \end{equation}
which can be shown to be independent of the choice of 
$m$. The negative half-line diffraction matrices satisfy
\begin{gather} {\mathsf{T}}_f^-\gamma + {\mathsf{R}}_b^-\boldsymbol{B}^{(0)}=\boldsymbol{A}^{(0)}, \end{gather}
\begin{gather}{\mathsf{R}}_f^-\gamma + {\mathsf{T}}_b^-\boldsymbol{B}^{(0)}=\eta, \end{gather}
where 
$\gamma$ and 
$\eta$ denote the amplitudes of the incident (right-travelling) and scattered (left-travelling) Bloch waves in the metamaterial region, respectively. Similarly to the right semi-infinite array problem, (2.22) may be reduced to a 
$2\times 2$ scattering matrix equation of the form
\begin{equation} \boldsymbol{\mathsf{S}}(d,0)\begin{bmatrix} \gamma\\ B_0^{(0)} \end{bmatrix} = \begin{bmatrix} A_0^{(0)}\\ \eta \end{bmatrix}, \end{equation}when evanescent modes are discarded.
2.5. Two coupled semi-infinite arrays
Lastly, we solve the problem of two coupled semi-infinite arrays. The submergence depth of the barriers of the semi-infinite array on the negative (positive) half-line is denoted 
$d^-$ (
$d^+$). For simplicity we will only consider the case when the spacing 
$W$ is identical for both semi-infinite arrays, although this condition would be easy to relax at the cost of introducing more notation. A schematic of the geometry is given in figure 1.
Figure 1. A schematic of the geometry of the problem of two coupled semi-infinite arrays of surface-piercing vertical barriers in a fluid of finite depth.
In general, we consider the case where incident Bloch waves can originate from both semi-infinite arrays. The Bloch wave incident from the left (right) semi-infinite array has complex amplitude 
$\alpha ^-$ (
$\beta ^+$) and propagates towards the right (left). The unknown amplitude of the scattered wave travelling left (right) in the left (right) semi-infinite array is denoted 
$\beta ^-$ (
$\alpha ^+$). With reference to the diffraction matrices in § 2.4, these unknown amplitudes can be obtained from the following linear system:
\begin{equation} \begin{bmatrix} -1 & 0 & {\mathsf{T}}_f^+ & \boldsymbol{0}_{V}^\intercal\\ \boldsymbol{0}_{V} & \boldsymbol{0}_{V} & {\mathsf{R}}_f^+ & -\boldsymbol{\mathsf{I}}\\ \boldsymbol{0}_{V} & \boldsymbol{0}_{V} & -\boldsymbol{\mathsf{I}} & {\mathsf{R}}_b^-\\ 0 & -1 & \boldsymbol{0}_{V}^\intercal & {\mathsf{T}}_b^- \end{bmatrix} \begin{bmatrix} \alpha^+\\ \beta^-\\ \boldsymbol{A}^{(0)}\\ \boldsymbol{B}^{(0)} \end{bmatrix} =\begin{bmatrix} -{\mathsf{R}}_b^+\beta^+\\-{\mathsf{T}}_b^+\beta^+\\-{\mathsf{T}}_f^- \alpha^-\\ -{\mathsf{R}}_f^- \alpha^- \end{bmatrix}, \end{equation}
in which 
$\intercal$ denotes the transpose. When 
$\alpha ^-=1$ and 
$\beta ^+=0$, the reflection and transmission coefficients for this problem are 
$R=\beta ^-$ and 
$T=\alpha ^+$. Green's second identity and the Bloch wave normalisation (2.11) yields the conservation of energy identity 
$|R|^2+|T|^2=1$. When 
$\alpha ^-=0$ and 
$\beta ^+=1$, the aforementioned identity again holds with the roles of 
$R$ and 
$T$ reversed. After inverting the matrix in (2.25) and using matrix algebra to eliminate 
$\boldsymbol {A}^{(0)}$ and 
$\boldsymbol {B}^{(0)}$, the following 
$2\times 2$ linear system for the unknown Bloch amplitudes is obtained:
\begin{equation} \boldsymbol{\mathsf{S}}(d^-,d^+)\begin{bmatrix} \alpha^-\\ \beta^+ \end{bmatrix}=\begin{bmatrix} \alpha^+\\ \beta^- \end{bmatrix}. \end{equation} The previous discussion assumed that propagating Bloch waves are defined in both semi-infinite regions. If 
$\omega$ is in the stopband of the right semi-infinite array, then (2.25) reduces to
\begin{equation} \begin{bmatrix} \boldsymbol{0}_{V} & {\mathsf{R}}_f^+ & -\boldsymbol{\mathsf{I}}\\ \boldsymbol{0}_{V} & -\boldsymbol{\mathsf{I}} & {\mathsf{R}}_b^-\\ -1 & \boldsymbol{0}_{V}^\intercal & {\mathsf{T}}_b^- \end{bmatrix} \begin{bmatrix} \beta^-\\ \boldsymbol{A}^{(0)}\\ \boldsymbol{B}^{(0)} \end{bmatrix} =\begin{bmatrix} \boldsymbol{0}_{V}\\-{\mathsf{T}}_f^- \alpha^-\\ -{\mathsf{R}}_f^- \alpha^- \end{bmatrix}. \end{equation}
For 
$\alpha ^-=1$, the conservation of energy condition becomes 
$|R|=1$, where 
$R=\beta ^-$. For notational consistency, we remark that (2.27) reduces to the following 
$1\times 1$ linear system:
\begin{equation} \boldsymbol{\mathsf{S}}(d^-,d^+)\begin{bmatrix} \alpha^- \end{bmatrix}=\begin{bmatrix} \beta^- \end{bmatrix}. \end{equation}3. Bloch wave coupling in semi-infinite arrays
In the remainder of this paper, we fix the parameters 
$H=20$ m and 
$W=2$ m. With regards to the numerical parameters, we set the truncation of the solution expansion 
$N_{sol}=20$. Furthermore, the point in the array where decaying effects are assumed to have vanished is fixed at 
$P=50$ for the remainder of this paper. This value of 
$P$ was chosen so that at least four figures of convergence was obtained for 
$\alpha$, where 
$\alpha$ is the amplitude of the right-travelling Bloch wave in the problem of forward scattering by a single semi-infinite array described in § 2.4. The validation of our method for the problem of two coupled semi-infinite arrays is discussed in the next subsection.
3.1. Method validation
Several steps were taken to validate the method given in § 2.5 and its subsequent implementation. First, the absolute error of the conservation of energy identities was checked and found to be less than 
$5 \times 10^{-7}$ across the range of parameters and frequencies considered here. Second, we chose an arbitrary selection of barriers on which to verify the boundary and matching conditions. We confirmed that to a good approximation (i) the normal derivative of the potential on each barrier was zero and (ii) the potential was continuously differentiable beneath the barrier. Third, we recognise that the solution to the two coupled semi-infinite arrays problem must reduce to the infinite array problem when the barriers of both arrays have the same submergence depth, i.e. 
$d^+=d^-$. This means that Bloch waves should not undergo any interaction at the interface and that 
$\boldsymbol{\mathsf{S}}(d^-,d^+)$ is the 
$2\times 2$ identity matrix. This has been verified.
3.2. Transmission between two semi-infinite arrays
We first recall some results from Wilks et al. (Reference Wilks, Montiel and Wakes2022) about Bloch wave dispersion in an infinite array of vertical barriers. Periodic arrays of vertical barriers have a single finite passband at low frequencies, which lies beneath an infinite stopband. The so-called cutoff frequency which separates the passband and stopband is closely related to the resonant frequency of a pair of vertical barriers, which decreases as the submergence depth of the two barriers increases (Wilks et al. Reference Wilks, Montiel and Wakes2022).
Curves showing the effect of frequency on the transmission of Bloch waves from the left semi-infinite array to the right semi-infinite array are given in figure 2. In that figure, we fix the submergence depth of the barriers in the right semi-infinite array to be 
$d^+=5$ m, while varying the submergence depth of the barriers in the left semi-infinite array from 
$d^-=0$ m to 
$d^-=5$ m. We observe that at low frequencies, transmission is very close to unity regardless of the submergence depth of the barriers in the left semi-infinite array. This is because the interaction between the barriers and the wave is negligible when 
$k_0d\ll 1$ and 
$d\ll H$ (Ursell Reference Ursell1947). For 
$d^-< d^+$, transmission decreases as the frequency approaches the cutoff frequency from below. The rate of this decrease is higher when 
$d^+-d^-$ is larger. When 
$0\ {\rm m}< d^-<5\ {\rm m}$ and 
$\omega$ is in the stopband of the right semi-infinite array (beginning at approximately 
$1.342$ s
$^{-1}$), Bloch waves cannot propagate to the right, thus 
$T$ vanishes as only reflection can occur. In these cases, 
$\omega$ eventually enters the stopband of the left semi-infinite array as well. When this happens, the problem is not well defined and the transmission curve terminates, although several of these terminations occur outside of the frame of figure 2. When 
$d^-=0$ m, the open sea region 
$x<0$ supports propagating plane waves at all frequencies, thus the corresponding transmission curve never terminates.
Figure 2. The effect of frequency on the transmission of Bloch waves from the left semi-infinite array to the right semi-infinite array. The submergence depth of the barriers in the right semi-infinite array is fixed at 
$d^+=5$ m, while the submergence depth of the barriers in the left semi-infinite array ranges from 
$d^-=1$ m to 
$d^-=5$ m. The perfect transmission indicated by the curve for 
$d^+=d^-$ provides validation for the method of coupling two semi-infinite arrays – see the discussion in § 3.1 for further details. When 
$d^-=0$ m, the problem reduces to scattering by a single semi-infinite array – the corresponding transmission curve is also included.
The results from this section have implications for the LBWA. In any unit cell of a graded array, the local Bloch waves must couple with the local Bloch waves in adjacent unit cells. It is reasonable to expect that this coupling will behave like the coupling of Bloch waves across two semi-infinite arrays. When the grading is not sufficiently weak, the difference between the submergence depths of adjacent barriers may not be small. This means that the naïve version of the LBWA, which treats Bloch wave coupling before the turning point as perfect transmission, may not apply. Instead, these reflections should be accounted for and the problem should be treated like a multiple scattering problem. In the next section, we describe and implement a version of the LBWA which accounts for reflections of Bloch waves yet continues to ignore the contributions of decaying Bloch modes.
4. Numerical implementation of the LBWA
We consider the problem of scattering by a finite graded array of 
$N$ vertical barriers using an approximate method based on the LBWA. In particular, we assume that between every adjacent pair of barriers, the wave can be represented as the superposition of right- and left-travelling Bloch waves. This assumption reduces the problem to the well-known problem of one-dimensional multiple scattering of waves on a string (see Martin (Reference Martin2014) for an outline of these problems). In this analogy, the grading of the submergence depths of the barriers in the water wave metamaterial would be equivalent to the string having piecewise-constant material properties. To describe our scattering problem in this way, we must obtain the transmission and reflection coefficients that couple Bloch waves in adjacent regions – we derive these from the coupling of two semi-infinite arrays as discussed in § 2.5. The solution obtained using this method will be contrasted with the solution obtained using the semianalytic method introduced in § 2.2.
4.1. Method
As in § 2.2, the barriers are positioned at 
$x=nW$ and have submergence depth 
$d^{(n)}$ for 
$n\in \{0,\ldots,N\}$. In the region 
$x<0$, we assume a plane incident wave of the form 
$\exp (\mathrm {i} k_0 x)\psi _0(z)$ and also approximate the scattered field as a plane wave, which implies that
\begin{gather} A_m^{(0)}=\tilde{A}_{Inc}\delta_{0m}, \end{gather}
\begin{gather}B_m^{(0)}\approx\tilde{B}_{R}\delta_{0m}, \end{gather}
where 
$\tilde {B}_{R}$ is unknown and
\begin{equation} \tilde{A}_{Inc}\approx A_{Inc}\frac{g}{\omega}\left(\frac{\sinh(2k_0H)}{4k_0H} +\frac{1}{2}\right)^{1/2}, \end{equation}
where 
$A_{Inc}$ is the amplitude of the incident wave. To the right of the array, the transmitted wave is approximated as a plane wave of the form
\begin{equation} A_m^{(N+1)} \approx \tilde{A}_{T}\delta_{0m}, \end{equation}
where 
$\tilde {A}_{T}$ is unknown. For 
$x\in (nW,(n+1)W)$ and 
$0\leq n < N$, the wave is approximated as the superposition of Bloch waves that correspond to the infinite array where the barriers have submergence 
$d^{(n)}$ if solutions to the Bloch–Floquet problem exist. That is, we assume that the potential is a superposition of right- and left-travelling Bloch waves with complex amplitudes 
$\alpha ^{(n)}$ and 
$\beta ^{(n)}$, respectively. This implies
\begin{gather} \boldsymbol{A}^{(n)} \approx\alpha^{(n)}\exp(\mathrm{i} q^{(n-1)} n W)\boldsymbol{A}^{(n-1)+}+\beta^{(n)}\exp(-\mathrm{i} q^{(n-1)} n W)\boldsymbol{A}^{(n-1)-}, \end{gather}
\begin{gather}\boldsymbol{B}^{(n)} \approx\alpha^{(n)}\exp(\mathrm{i} q^{(n-1)} n W)\boldsymbol{B}^{(n-1)+}+\beta^{(n)}\exp(-\mathrm{i} q^{(n-1)} n W)\boldsymbol{B}^{(n-1)-}, \end{gather}
where 
$q^{(j)}$ is the Bloch wavenumber and 
$\boldsymbol {A}^{(j)\pm }$ and 
$\boldsymbol {B}^{(j)\pm }$ are the Bloch wave coefficients for the infinite array of barriers with submergence depth 
$d^{(j)}$. These quantities are computed using the method given in § 2.3. If the Bloch–Floquet problem has no solutions for the infinite array where the barriers have submergence 
$d^{(n)}$, then we assume 
$A_m^{(n)}=B_m^{(n)}=0$.
To compute the LBWA, we must consider three cases.
4.1.1. Case 1: total reflection by the first barrier
This case governs the situation where the region 
$0< x< W$ does not support propagating Bloch waves, that is, the frequency is in the stopband of the infinite array where the barriers have submergence 
$d^{(0)}$. The LBWA suggests that there is no wave amplitude anywhere in the array and that transmission by the array is zero. The interaction at the first barrier is modelled using the single semi-infinite array diffraction problem where the barriers occupy the positive half-line. Thus, we need only compute 
$\tilde {B}_{R}=\boldsymbol{\mathsf{S}}(0,d^{(0)})\tilde {A}_{Inc}$.
4.1.2. Case 2: total reflection within the graded array
This case describes the situation where the local Bloch wave is totally reflected inside the array. Let 
$p$ be the smallest value in 
$\{1,\ldots,N\}$ such that 
$\omega$ is in the stopband of the infinite array of barriers with submergence depth 
$d^{(p)}$. This implies that 
$x=pW$ is the turning point and no Bloch waves exist past this point. At the leftmost barrier of the array, we assume that the wave interaction can be locally modelled using the single semi-infinite array diffraction problem. Thus, we use (2.19) to relate the amplitudes of plane waves to the left and Bloch waves to the right of the leftmost barrier, i.e.
\begin{equation} \boldsymbol{\mathsf{S}}(0,d^{(0)})\begin{bmatrix} \tilde{A}_{Inc}\\ \beta^{(1)} \end{bmatrix}=\begin{bmatrix} \alpha^{(1)}\\ \tilde{B}_R \end{bmatrix}. \end{equation}
For 
$2\leq n\leq p$, we assume that the wave interaction can be locally modelled using the two coupled semi-infinite arrays problem. This problem must be transformed so that the first barrier of the right semi-infinite array occurs at 
$x=(n-1)W$ instead of at 
$x=0$. By combining (2.26) with the phase shift matrices 
${\mathsf{P}}_n=\mathrm {diag}[\exp (-\mathrm {i} q^{(n-1)}nW),\exp (\mathrm {i} q^{(n-2)}nW)]$ and 
${\mathsf{Q}}_n=\mathrm {diag}[\exp (-\mathrm {i} q^{(n-1)}nW),\exp (\mathrm {i} q^{(n-2)}nW)]$, we obtain
\begin{equation} {\mathsf{P}}_n\boldsymbol{\mathsf{S}}(d^{(n-2)},d^{(n-1)}){\mathsf{Q}}_n\begin{bmatrix} \alpha^{(n-1)}\\ \beta^{(n)} \end{bmatrix}=\begin{bmatrix} \alpha^{(n)}\\ \beta^{(n-1)} \end{bmatrix}. \end{equation}Lastly, the amplitudes of the Bloch waves in the region immediately preceding the turning-point barrier are again related by assuming that the wave interaction can be locally modelled using the two coupled semi-infinite arrays problem. However, in this problem no Bloch wave is defined in the right semi-infinite array. Thus, we use (2.28) to model the reflection of Bloch waves in the left semi-infinite array, which after transformation becomes
\begin{equation} \exp(2\mathrm{i} q^{(p-1)}pW)\boldsymbol{\mathsf{S}}(d^{(p-1)},d^{(p)})\alpha^{(p)}=\beta^{(p)}. \end{equation}
Equation (4.3a–c) form a system of 
$2p+1$ equations and 
$2p+1$ unknowns, which can be solved using standard methods. The approximation of the fluid potential can be recovered using (4.1a,b) and (4.2a,b).
4.1.3. Case 3: partial transmission
This last case governs the case where local Bloch waves are defined throughout the array and therefore the array can transmit energy. The solution proceeds similarly to case 2, where (4.3a) holds and (4.3b) holds for 
$2\leq n \leq N$. However, (4.3c) no longer applies. Instead, the interaction of the Bloch waves to the left of the 
$(N+1)$th barrier with the plane wave to the right is modelled using the semi-infinite array diffraction problem where the barriers occupy the negative half-line. This implies that
\begin{equation} {\mathsf{P}}_{N+1}\boldsymbol{\mathsf{S}}(d^{(N)},0){\mathsf{Q}}_{N+1}\begin{bmatrix} \alpha^{(N)}\\ 0 \end{bmatrix} = \begin{bmatrix} \tilde{A}_{T} \\ \beta^{(N)} \end{bmatrix} \end{equation}where the phase-shift matrices have been defined as
\begin{align} {\mathsf{P}}_{N+1}=\mathrm{diag}[1,\exp(\mathrm{i} q^{(N)}(N+1)W]\quad\mbox{and}\quad {\mathsf{Q}}_{N+1}=\mathrm{diag}[\exp(\mathrm{i} q^{(N)}(N+1)W,1]. \end{align}
Case 3 requires that 
$d^{(N-1)}=d^{(N)}$ so that Bloch waves can be consistently defined in the region 
$x\in ((N-1)W,NW)$.
4.2. Validation
The purpose of this subsection is to validate our implementation of the LBWA by comparing its predictions with the corresponding semianalytic solutions. In particular, we want to show that our implementation does not suffer from mathematical or coding mistakes, since later (in § 4.3) we will seek to identify and describe errors that are fundamental to the LBWA. With this in mind, we will not attempt to validate our implementation of the LBWA using a graded array, since this would not allow us to identify the origins of any discrepancies. Instead, we rule out the presence of mathematical or coding errors by first considering a case in which the LBWA should perform well, namely, that of a long array of regularly spaced and identical vertical barriers. Indeed, the LBWA reduces to a long array approximation analogous to that of Thompson et al. (Reference Thompson, Linton and Porter2008), as it assumes that only propagating Bloch waves excited at each of the array ends can interact at the opposite ends. In other words, the evanescent modes excited at each end do not interact. To evaluate the LBWA over a range of frequencies, we consider the absolute error of the reflection coefficient defined as 
$E_R = |\tilde {R}-R|$, where 
$\tilde {R}$ is computed from the LBWA (i.e. it is the value of 
$\tilde {B}$ when 
$\tilde {A}_{Inc}=1$) and 
$R$ is computed using the semianalytic method described in § 2.2. In figure 3, results are presented for a finite array of 
$N=50$ barriers with submergence depth 
$d^{(n)}=5$ m for all 
$0\leq n\leq N$. In figure 3(a), we observe that the absolute error of the reflection coefficient remains bounded below 
$5\times 10^{-7}$ below the cutoff frequency 
$\omega \approx 1.342$ s
$^{-1}$ and is numerically zero above the cutoff.
Figure 3. Validation of the LBWA for a finite array of vertical barriers for which 
$N=50$ and 
$d^{(n)}=5$ m for all 
$0\leq n\leq N$. Panel (a) displays the absolute error of the reflection coefficient. Panel (b) contrasts the absolute value of the free surface elevation computed using the LBWA (blue line) and using the semianalytic method (red line) at the angular frequency 
$\omega =1.2$ s
$^{-1}$. The 
$x$-coordinates of the vertical barriers are marked with dashed black lines.
We also assess how well the LBWA predicts local energy amplification within the array. To do this, we use the complex-valued free surface elevation given by
\begin{equation} \zeta(x) = \frac{\mathrm{i}\omega}{g}\phi(x,0). \end{equation}
In figure 3(b), the absolute value of the free surface elevation at 
$\omega =1.2$ s
$^{-1}$ is computed using both the LBWA (blue line) and the semianalytic solution (red line) for the array. We observe that there is excellent agreement between the two solutions across most of the horizontal domain. Small errors occur in the neighbourhood of the zeroth and 
$N$th barriers, which are due to evanescent modes excited at the array ends that are not considered by the LBWA. These results are completely analogous to those of Thompson et al. (Reference Thompson, Linton and Porter2008). We conclude that our method performs as expected for a problem that the LBWA should describe accurately. This suggests that any errors produced by our method in other problems are fundamental to the LBWA.
4.3. Results
We consider linearly graded arrays of 
$N$ vertical barriers described by 
${\rm \Delta} d = d^{(n)}-d^{(n-1)}$, where we restrict 
${\rm \Delta} d$ so that 
$N+1=d^{(N)}/{\rm \Delta} d\in \mathbb {N}$. Moreover, we fix 
$d^{(N)}=10$ m which implies that 
$d^{(0)}={\rm \Delta} d$. Three such arrays are considered, for which the grading parameters are 
${\rm \Delta} d=0.25$ m, 
${\rm \Delta} d=0.5$ m and 
${\rm \Delta} d = 1$ m. In this subsection, we restrict our investigation to case 2, i.e. we only consider frequencies at which the LBWA predicts total reflection within the array. Plots showing the absolute error of the reflection coefficient for each of these arrays are given in figure 4. We observe that the curves of 
$E_R$ are sawtooth-like (recall that 
$E_R$ is the absolute error of the reflection coefficient defined in § 4.2). The peaks of these error curves occur just above the angular frequencies 
$\omega ^{(n)}$, where 
$\omega ^{(n)}$ is the cutoff frequency of the infinite array of vertical barriers with submergence depth 
$d^{(n)}$. Away from these peaks, the error approaches a baseline that is lower for more weakly graded arrays. The error peaks for more weakly graded arrays are more densely packed because these arrays contain more different barrier submergence depths, i.e. they have more internal cutoff frequencies.
Figure 4. The absolute error of the reflection coefficient for arrays of vertical barriers parameterised by (a) 
${\rm \Delta} d=1$ m, (b) 
${\rm \Delta} d=0.5$ m and (c) 
${\rm \Delta} d=0.25$ m. The cutoff frequencies 
$\omega ^{(n)}$ are indicated with dashed vertical lines. Specifically, these are the cutoff frequencies for an infinite array of vertical barriers with spacing 
$W$ and submergence 
$d^{(n)}$, for 
$0\leq n \leq N$.
Next, we seek to explain why the cutoff frequencies induce peaks in the error curves. To do this, we use the absolute value of the free surface to compare the LBWA and semianalytic solutions for the graded array with 
${\rm \Delta} d=0.25$ m (figure 5) and 
${\rm \Delta} d=1$ m (figure 6). We compare the solutions at 
$\omega =1.341$ s
$^{-1}$ (figure 5a and figure 6a) and 
$\omega =1.343$ s
$^{-1}$ (figure 5b and figure 6b). The former angular frequency is just below the cutoff frequency 
$\omega ^{(j)}$ associated with the barrier with submergence depth 
$d^{(j)}=5$ m, whereas the latter angular frequency is just above this cutoff. The index 
$j$ of this barrier depends on 
${\rm \Delta} d$. Specifically, 
$j=19$ in figure 5 and 
$j=4$ in figure 6.
Figure 5. The free surface elevation in the graded array of vertical barriers with grading parameter 
${\rm \Delta} d=0.25$ m and angular frequencies (a) 
$\omega =1.341$ s
$^{-1}$ and (b) 
$\omega =1.343$ s
$^{-1}$. A complete description of the array geometry is provided in § 4.3. The elevation curves were computed using the semianalytic method (red line) and our implementation of the LBWA (blue line). The positions of the vertical barriers (i.e. 
$x=nW$ for 
$0\leq n \leq N$) are marked with dashed black lines and the position of the turning point (i.e. 
$x=pW$) is marked with a solid black line.
Figure 6. As for figure 5 but with 
${\rm \Delta} d = 1$ m.
In figures 5(a) and 6(a), which correspond to the case where 
$\omega$ is just below 
$\omega ^{(j)}$, we observe excellent agreement between the solutions for 
$x< pW$, where 
$p=j+1$ is the turning point. In results not shown here, similarly good agreement was found across a range of different array configurations and frequencies, provided that the frequency is away from the peaks in 
$E_R$. This suggests that decaying Bloch waves are negligible in these cases.
We briefly discuss the small discrepancy in figure 5(a) in the region 
$pW< x<(p+1)W$. As this region does not support propagating Bloch modes, our implementation of the LBWA defines the amplitude in this region to be zero. However, this region does support decaying Bloch waves. To explain why this error does not significantly affect the accuracy of the LBWA for 
$x< pW$, we note that the reflection coefficient used by the LBWA at 
$x=pW$ was computed from the version of the two coupled semi-infinite arrays problem in which there is no propagating mode in the right array. This means that the decaying Bloch wave in 
$pW< x<(p+1)W$ is internally accounted for by the LBWA, despite the fact that we later assume that there is zero amplitude in this region. Because 
$\omega$ is relatively far above 
$\omega ^{(p)}$, the decaying Bloch mode in this region attenuates rapidly (see Appendix for details). The rapid attenuation means that the effect of this mode on the solution is negligible.
In figures 5(b) and 6(b), which correspond to the case where 
$\omega$ is just above 
$\omega ^{(j)}$ (i.e. 
$\omega$ is a peak of 
$E_R$), we observe that the agreement between the semianalytic and LBWA solutions is significantly worse. As was the case in figure 5(a), the most notable difference occurs in the region 
$pW< x<(p+1)W$. However, the increased value of 
$\omega$ means that 
$p$ has shifted from 
$j+1$ to 
$j$. Since 
$\omega$ is now very close to 
$\omega ^{(p)}$, the decaying Bloch mode in the region 
$pW< x<(p+1)W$ attenuates slowly. We hypothesise that this slow attenuation allows the decaying wave to interact with the barriers at 
$x=nW$ for 
$n>p$ in a non-negligible way, which ultimately affects the solution.
To support this hypothesis, we consider a partially graded array of vertical barriers. The partial grading is described by
\begin{equation} d^{(n)}=\begin{cases} (n+1){\rm \Delta} d & \mbox{for }0\leq n< 20 ,\\ 5\ \mbox{m} & \mbox{for }20\leq n\leq 40 , \end{cases} \end{equation}
where 
${\rm \Delta} d=0.25$ m. That is, the submergence depth of the barriers increases linearly up to 
$5$ m, after which it remains constant. At 
$\omega =1.343$ s
$^{-1}$, the LBWA assumes that waves cannot propagate beyond the the 
$19$th barrier, since 
$\omega$ is in the stopband of this region. Importantly, this means that at 
$\omega$, our implementation of the LBWA is unable to distinguish between this partially graded array and the fully graded array with 
${\rm \Delta} d=0.25$ m. On the other hand, the semianalytic method can distinguish between these two cases. The free surface elevation in the partially graded array computed using both methods is given in figure 7. We observe excellent agreement between the two methods to the left of the turning point, despite the LBWA failing to predict the free surface elevation accurately for the fully graded array at the same frequency. This highlights the importance of the part of the array to the right of the turning point when the decaying Bloch wave in 
$pW< x<(p+1)W$ attenuates slowly. Indeed, if the part of the array to the right of the turning point is not well approximated by a semi-infinite array (as the LBWA assumes), these decaying Bloch waves can be reflected into the region to the left of the turning point and affect the solution there. This suggests that the peaks of 
$E_R$ could be reduced by replacing (4.3c) with an equation that accounts for the grading of the barriers beyond the turning point, the derivation of which is not considered here.
Figure 7. The free surface elevation in the graded array of vertical barriers described by (4.7). As in figure 5, the elevation curves were computed using the semianalytic method (red line) and our implementation of the LBWA (blue line). The positions of the vertical barriers are marked with dashed black lines and the position of the turning point is marked with a solid black line.
5. Concluding remarks
The LBWA has previously been established as a powerful tool for qualitatively understanding wave propagation through graded arrays. This paper has examined the LBWA in the context of graded arrays of vertical barriers by exploring the coupling mechanism between Bloch waves in adjacent unit cells. This mechanism was assumed to behave like the scattering of Bloch waves across the interface between two semi-infinite arrays of vertical barriers. We found that reflection by these interfaces is non-negligible when the differences between the submergence depths of the barriers of the two arrays is not small. This suggests that these reflections should be considered whenever the grading parameter of the metamaterial is significant.
We then described a numerical implementation of the LBWA which accounts for these small reflections yet continues to omit decaying Bloch modes. The method accurately predicts the amplitude of the free surface at a wide range of frequencies. Exceptions occur in frequency regions which lie just above the cutoff frequencies of individual barriers in the array, where error peaks occur. We argued that these errors result from the LBWA failing to account for decaying Bloch waves in the region just beyond the turning point, which can have a non-negligible effect if they attenuate slowly. It may be possible to create a more accurate implementation of the LBWA by accounting for these attenuating waves.
We note that the large amplitude oscillations reported in this paper may not be realistic, since linear water wave theory is only accurate for low amplitude waves (Mei et al. Reference Mei, Stiassnie and Yue2005) and the boundary conditions on the barriers are also greatly idealised when compared with real-world materials. The large responses in the array are related to the fundamental resonance of a pair of surface-piercing vertical barriers in isolation, which is known to become stronger as the ratio of the barrier submergence depth to the spacing increases (Newman Reference Newman1974; McIver Reference McIver1985; Wilks et al. Reference Wilks, Montiel and Wakes2022). Experimental or numerical tests (i.e. using wave flumes or computational fluid dynamics, respectively) would be required to evaluate the accuracy of the model.
The method presented in this paper may extend to graded arrays in which the scatterers are not vertical barriers, although further investigation will be required to confirm this. For example, it may be possible to extend our investigation of the LBWA to study the coupling of Rayleigh–Bloch waves in a finite, graded line array of scatterers in an unbounded, two-dimensional medium (in the sense of Bennetts et al. (Reference Bennetts, Peter and Craster2018)). However, it is likely that mathematical difficulties due to the underlying spectral theory would be encountered because the periodicity cells would be unbounded. It may also be possible to implement the LBWA to study the interaction of waves with two-dimensional graded arrays of scatterers in an unbounded, two-dimensional medium (in the sense of Bennetts et al. (Reference Bennetts, Peter and Craster2019)). The LBWA suggests that in the neighbourhood of each scatterer, the wave field is governed by the Bloch waves of the cognate doubly periodic lattices. Although the periodicity cells would be bounded in this case, difficulties are likely to arise because in general there is a continuous directional spectrum of the lattice Bloch waves, whereas only forward and backward-propagating Bloch waves needed to be considered in this paper. Despite the potential difficulties, both of these extensions could lead to better understanding of wave interaction with graded metamaterials.
Declaration of interests
The authors report no conflict of interest.
Appendix. Bloch waves above the cutoff
The following discussion is heavily based on an article by Bennetts & Peter (Reference Bennetts and Peter2022), who studied Rayleigh–Bloch waves above the cutoff in a medium governed by the Helmholtz equation that contains an infinite line arrays of circular cylinders (also see Bennetts, Peter & Montiel Reference Bennetts, Peter and Montiel2017; Bennetts et al. Reference Bennetts, Peter and Craster2018, Reference Bennetts, Peter and Craster2019). Our task is to show that in an infinite array of vertical barriers, the attenuation rate of Bloch waves above the cutoff is an increasing function of frequency.
Recall that 
$\exp (\pm \mathrm {i} q W)$ are eigenvalues of (2.10). Below the cutoff, 
$q$ is real and therefore these eigenvalues lie on the unit circle. The eigenvalue in the upper (lower) half-plane corresponds to the right-propagating (left-propagating) Bloch mode. As 
$\omega$ approaches the cutoff, 
$q$ approaches 
${\rm \pi} /W$ and both eigenvalues move around the unit circle, eventually merging at 
$-1$. Above the cutoff, the merged eigenvalues bifurcate as a real-valued reciprocal pair, with one moving inside the unit circle and the other moving outside the unit circle. The eigenmode corresponding to the eigenvalue on the inside (outside) of the unit circle describes a Bloch wave which decays towards the right (left). The decay rate of the Bloch wave above the cutoff is described by the imaginary part of 
$q$. In figure 8, we observe that above the cutoff, the imaginary part of 
$q$ is an increasing function of 
$\omega$, which is what we wanted to show. Although 
$d=5$ m and 
$W=2$ m have been fixed in figures 9 and 8, we have verified that the result holds for a range of other values of the parameters.
Figure 8. Plots showing (a) the real part of 
$q$ and (b) the imaginary part of 
$q$ as functions of the angular frequency.
Figure 9. Generalised eigenvalues of (2.10) for the infinite array of vertical barriers with 
$d=5$ m and 
$W=2$ m at (a) 
$\omega =1.2$ s
$^{-1}$, (b) 
$\omega =1.34$ s
$^{-1}$ and (c) 
$\omega =1.344$ s
$^{-1}$. The values 
$\exp (\pm \mathrm {i} q W)$ are marked with red dots and all other eigenvalues are marked with blue dots. The unit circle is marked with a dashed black line.
Open access
