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:

(2.1a)\begin{gather} \bigtriangleup \phi =0 \quad (x,z)\in\varOmega, \end{gather}

(2.1b)\begin{gather}\partial_z\phi=\frac{\omega^2}{g}\phi \quad z=0, \end{gather}

(2.1c)\begin{gather}\partial_z \phi=0 \quad z={-}H, \end{gather}

(2.1d)\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}

(2.1e)\begin{gather}\partial_x\phi=0 \quad (x,z)\in\varGamma, \end{gather}

(2.1f)\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

(2.2)\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

(2.3)\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

(2.4)\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

(2.5)\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

(2.6)\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

(2.7)\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

(2.8)\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

(2.9a)\begin{gather} A_m^{(n)}=\exp({\pm}\mathrm{i} q n W)A_m^{{\pm}}, \end{gather}

(2.9b)\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:

(2.10)\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

(2.11)\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

(2.12)\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

(2.13a)\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}

(2.13b)\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

(2.14a)\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}

(2.14b)\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

(2.15a)\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}

(2.15b)\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:

(2.16a)\begin{gather} \boldsymbol{\tilde{A}}^{(n)}= \sum_{j=0}^n \exp(\mathrm{i} q (n-j)W)\boldsymbol{\tilde{C}}^{(j)}, \end{gather}

(2.16b)\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

(2.17a)\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}

(2.17b)\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

(2.17c)\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}

(2.17d)\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

(2.18a)\begin{gather} {\mathsf{T}}_f^+\boldsymbol{A}^{(0)}+{\mathsf{R}}_b^+\beta = \alpha , \end{gather}

(2.18b)\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

(2.19)\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:

(2.20a)\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}

(2.20b)\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

(2.21)\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

(2.22)\begin{gather} {\mathsf{T}}_f^-\gamma + {\mathsf{R}}_b^-\boldsymbol{B}^{(0)}=\boldsymbol{A}^{(0)}, \end{gather}

(2.23)\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

(2.24)\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:

(2.25)\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:

(2.26)\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

(2.27)\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:

(2.28)\begin{equation} \boldsymbol{\mathsf{S}}(d^-,d^+)\begin{bmatrix} \alpha^- \end{bmatrix}=\begin{bmatrix} \beta^- \end{bmatrix}. \end{equation}

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.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

(4.6)\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.