1. Introduction
Problems involving the reflection and transmission of waves by thin vertical surface-piercing barriers under classical linearised theory have been the subject of research over many decades. For a fluid of infinite depth, Ursell (Reference Ursell1947) obtained an explicit solution for monochromatic waves normally incident upon a single barrier. When the fluid has a constant finite depth, Porter & Evans (Reference Porter and Evans1995) showed how to obtain accurate numerical solutions which provide upper and lower bounds on reflected and transmitted amplitudes by formulating complementary integral equations. They also produced numerical results for a pair of identical surface-piercing barriers, reproducing and extending results of Evans & Morris (Reference Evans and Morris1972), Newman (Reference Newman1974) and McIver (Reference McIver1985). Notably, Newman (Reference Newman1974) had shown, using matched asymptotic methods, that a pair of closely spaced barriers can totally transmit or reflect incoming wave energy at angular frequencies in the vicinity of a critical value of 
$\omega _c =\sqrt {g/d}$ related to the vertical fluid resonance in the narrow column between the two vertical plates, where 
$d$ is the depth of submergence of the plates and 
$g$ is the acceleration due to gravity. Evans (Reference Evans1978) later used this idea to model the operation of a narrow oscillating water column wave energy device. Non-identical barriers and arrays of more than two barriers have been considered by a number of authors including Evans & Porter (Reference Evans and Porter1997) and Roy, De & Mandal (Reference Roy, De and Mandal2019).
More recently, Wilks, Montiel & Wakes (Reference Wilks, Montiel and Wakes2022) have considered larger arrays of closely spaced barriers whose submergence increases gradually in the direction of the incident wave, i.e. so-called ‘graded array’. The graded array is a specific type of metamaterial that refers to a material designed to have specific properties not found in naturally occurring materials and typically consists of repeating sub-wavelength structures. The graded array in Wilks et al. (Reference Wilks, Montiel and Wakes2022) is designed to be resonant at multiple frequencies associated with the variable fluid column lengths across the array such that high reflection is sustained across a broad range of frequencies, known as rainbow reflection. Similar ideas based on local internal resonance provided by the displacement of the surface have been implemented in water waves using C-ring cylinders in channels, rather than vertical barriers, by Dupont et al. (Reference Dupont, Remy, Kimmoun, Molin, Guenneau and Enoch2017), Bennetts, Peter & Craster (Reference Bennetts, Peter and Craster2018), Archer et al. (Reference Archer, Wolgamot, Orszaghova, Bennetts, Peter and Craster2020), etc. Large arrays of floating buoys which resonate on the surface of the fluid with elements that either possess constant properties or vary in space have been considered for wave energy harvesting applications by, for example, Garnaud & Mei (Reference Garnaud and Mei2009) and Porter (Reference Porter2021), as well as by Wilks et al. (Reference Wilks, Montiel and Wakes2022). Arrays of resonators have been used to produce similar rainbow reflection on acoustic wave transmission in waveguides in the form of cavities attached to sidewalls (e.g. Tang Reference Tang2012; Jiménez et al. Reference Jiménez, Romero-Garcia, Pagneux and Groby2017; Jan & Porter Reference Jan and Porter2018) and on surface wave propagation in elasticity (e.g. Colquitt et al. Reference Colquitt, Colombi, Craster, Roux and Guenneau2017) in the form of mechanical oscillators attached to the surface of an elastic half-space.
In this paper, we return to consider large periodic arrays of surface-piercing vertical barriers which are equally spaced and submerged to the same constant depth, being simpler than the graded array problem. Part of the motivation for looking at this problem is to characterise the effects of resonance on models which are based on low-frequency homogenisation of plate array structures. This approximation replaces the discrete structure of the array by a continuous effective medium under an assumed contrast in scales between the wavelength and the spacing between plates. This approximation allows problems of a large array of closely spaced plates to be solved more easily and has been used to consider the interaction of waves with metamaterial structures (e.g. Jan & Porter Reference Jan and Porter2018; Liu et al. Reference Liu, Liang, Chen and Zhua2018; Zheng, Porter & Greaves Reference Zheng, Porter and Greaves2020). Away from resonant conditions, results from homogenisation of plate array structures compare favourably with those derived from the application of direct numerical methods for discrete arrays (e.g. Zheng et al. (Reference Zheng, Porter and Greaves2020) and Porter, Zheng & Liang (Reference Porter, Zheng and Liang2022), where comparison is made with boundary integral methods). However, low-frequency homogenisation fails when undamped fluid motion in narrow channels is resonant, occurring at the critical frequency 
$\omega _c$ indicated above. This is an example of the well-known phenomenon of local resonance in metamaterials (e.g. Ma & Sheng Reference Ma and Sheng2016).
By assuming a periodic array with barriers having an equal depth of submergence, 
$d$, we allow ourselves the opportunity of making an analytic comparison between a direct solution of the discrete array of 
$N+1$ plates separated by a non-zero distance 
$\delta$ and homogenisation methods based on 
$\delta /d \ll 1$ in order to understand issues relating to resonance. In particular, periodicity allows us to exploit Bloch–Floquet methods which are associated with the corresponding infinite periodic arrays. When considering the scattering by a finite number of identical elements arranged periodically in an array, a naive approach is to apply direct multiple scattering methods, which leads to a large coupled system (examples are given in Linton & McIver (Reference Linton and McIver2001)). For quasi-one-dimensional scattering systems, such as the one we consider here, transfer and scattering matrices are often used to reduce the computational effort, in which wave information is propagated left and right (see Porter & Porter Reference Porter and Porter2003; Wilks et al. Reference Wilks, Montiel and Wakes2022). For the transfer matrix method, reflection and transmission across the array are expressed as products of matrices which encode the scattering characteristics of a single element by converting incoming modes consisting of both propagating and a truncated set of evanescent waves into outgoing modes. This becomes impractical when considering closely spaced arrays since the size of the matrix must increase to capture a greater number of mode interactions as the value of 
$\delta /d$ is decreased. Results showed that for an infinite periodic array there exist ranges of frequencies for which wave propagation through the structure is prohibited. For the constrained mass system in Wilks et al. (Reference Wilks, Montiel and Wakes2022), they argued that they must replace scattering matrices in favour of a formulation involving a system of integral equations which link certain functions at the 
$n$th barrier to corresponding functions at the 
$(n\pm 1)$th barriers since the wave scattering problem and the equations of motion are needed to be solved simultaneously.
The frequency ranges indicated above for which wave propagation through the periodic structures is absolutely forbidden are well known as stop bands and the corresponding structure is called the band-gap structure. It is shown that the existence of stop bands is possible for one-dimensional periodically varying topography (e.g. Porter & Porter Reference Porter and Porter2003; An & Ye Reference An and Ye2004). Further, Chen et al. (Reference Chen, Kuo, Ye and Sun2004) and Yang et al. (Reference Yang, Wu, Zhong and Zhong2006) investigated the band-gap structures of liquid surface waves propagating through an infinite two-dimensional periodic topography of circle and square hollows, respectively. McIver (Reference McIver2000) and Linton (Reference Linton2011) also established the existence of a band-gap structure associated with water waves propagating over infinite periodic arrays of submerged vertical and horizontal cylinders. If an incident wave is subjected to a finite section of this infinite array at a certain frequency within such a stop band, most of energy is expected to be reflected.
In this paper we take a new approach to determining the scattering by periodic arrays of finite extent and develop concepts introduced in the work of Porter & Porter (Reference Porter and Porter2003) who highlighted connections between finite- and infinite-array problems. Specifically, they showed that the transfer matrices referred to above could be expressed in terms of generalised eigenvalues and eigenfunctions of the corresponding periodic Bloch–Floquet problem representing both propagating and decaying modes. After specifying the infinite periodic array problem in § 2 we introduce a novel orthogonality relation satisfied by the Bloch eigenfunctions, which is key to developing a solution in the interior of the array. Then we present two different forms of solution for the Bloch–Floquet problem by expanding the velocity potential with different orthogonal functions, each of which has its own advantages in numerical calculation. Information is now able to propagate from the left to the right of the array of 
$N+1$ barriers via a simple product of Bloch wavenumbers/eigenvalues; this is described in § 4 of the paper.
In § 3 we formally develop the homogenisation approximation for a continuum model, which relies on a separation of horizontal length scales based on the wavelength and the array separation, 
$\delta$. It turns out this assumption is violated as the resonance of fluid in the narrow channels is approached since solutions are predicted with a propagating wavelength which tends to zero. The wave scattering problem of this continuum model is solved in § 5. Numerical results are produced to show the comparison between the exact description of the finite array and the approximation based on homogenisation in order to demonstrate how the Bloch–Floquet solution produces accurate and efficient results and to demonstrate when homogenisation is a reliable approach to take. The work is summarised in § 6.
It is worth pointing out that as the spacing between plates tends to zero, viscous effects become important in the physical setting. Since our primary interest is to understand the resonance occurring in metamaterial structures and to evaluate the validity of homogenisation, the effect of viscosity is not taken into account in the present paper. For a consideration of the viscous effects the reader is referred to Mei, Stiassnie & Yue (Reference Mei, Stiassnie and Yue2005), whose idea has been implemented to consider the energy dissipation in the metamaterial (see Zheng et al. Reference Zheng, Porter and Greaves2020).
2. The periodic barrier problem
Consider an infinite periodic array of thin barriers with equal spacing 
$\delta$, which extend vertically downwards to a depth 
$d$ below the free surface of a fluid with constant depth 
$h$. Two-dimensional Cartesian coordinates are defined with the origin 
$O$ in the mean free surface and 
$z$ directed vertically upwards, such that barriers occupy 
$\{x=x_n = n \delta (n \in \mathbb {Z}),\ -d< z<0\}$.
Under the assumption of an incompressible and inviscid fluid and irrotational flow, the fluid motion of a single frequency 
$\omega$ can be described by 
$\mathrm {Re} \{ \phi (x,z) \exp ({-\mathrm {i} \omega t}) \}$ (in which 
$t$ is time), where the spatial velocity potential 
$\phi (x,z)$ satisfies
\begin{equation} \nabla^2 \phi = 0,\quad \mbox{in the fluid.}\end{equation}On the free surface, the combined linearised kinematic and dynamic boundary conditions result in
\begin{equation} K \phi - \phi_z = 0,\quad \mbox{on}\ z=0, \end{equation}
where 
$K = \omega ^2/g$. No-flow conditions apply on fixed rigid boundaries meaning
\begin{equation} \phi_x = 0,\quad \mbox{on}\ x=n \delta^\pm \quad (n \in \mathbb{Z}) \quad \mbox{for}\ -d < z< 0, \end{equation}where the positive and negative signs denote two sides of the plate, and
\begin{equation} \phi_z = 0, \quad \mbox{on } z=-h. \end{equation}
Besides, the analysis of the flow close to the edge of the barrier reveals that the fluid velocity should possess inverse square root behaviour. Finally, since the geometry is periodic we may invoke the Bloch–Floquet theory which allows us to consider just one periodic element of the array (here we choose 
$D = \{0 < x < \delta,\ -h < z< 0\}$) provided we introduce quasi-periodic boundary conditions on those fluid interfaces which connect one cell to the next. For our choice of 
$D$, we require
\begin{equation} \left.\begin{gathered} \phi(\delta,z)= {\rm e}^{\mathrm{i} \beta \delta} \phi(0,z)\\ \phi_x(\delta,z)= {\rm e}^{\mathrm{i} \beta \delta} \phi_x(0,z) \end{gathered}\right\} \quad \mbox{for } -h < z <-d, \end{equation}
where 
$\beta$ is the Bloch wavenumber needed to be determined.
Momentarily it helps to imagine that 
$\lambda$ replaces 
$\textrm {e}^{\mathrm {i} \beta \delta }$ in (2.5). Since the problem is unchanged by the mapping 
$x \to \delta - x$, if 
$\lambda$ is an eigenvalue then 
$\lambda ^{-1}$ is another eigenvalue. Also, if 
$\phi$ is a solution corresponding to 
$\lambda$ then 
$\bar {\phi }$, the complex conjugate of 
$\phi$, is also a solution with eigenvalue 
$\bar {\lambda }$. This implies that eigenvalues 
$\lambda$ must lie either on the real axis in reciprocal pairs or on the unit circle in complex conjugate pairs. Returning to 
$\beta$, this implies that 
$\beta$ is a real number or can be expressed as 
$n{\rm \pi} /\delta +\mathrm {i} \gamma _n$ (where 
$n \in \mathbb {Z}$ and 
$\gamma \in \mathbb {R}$) and that for every 
$\beta$ it is paired with 
$-\beta$.
When 
$\beta$ is real, it represents a wavenumber which encodes the phase shift of the fluid motion as waves propagate across one cell of the infinite periodic array. We note that we need only consider real values of 
$\beta \in (0,{\rm \pi} /\delta ]$ since 
$\beta ' = \beta + 2{\rm \pi} m/\delta$ for 
$m \in \mathbb {Z}$ leaves the problem unchanged and 
$\beta ' = 2 {\rm \pi}/\delta - \beta$ results in the same problem with 
$x$ mapped to 
$\delta - x$. That is to say, it reverses the wave direction, but not the solution. When 
$\beta$ becomes a complex number, for similar reasons we only need to consider the case of 
$\beta$ being 
$n{\rm \pi} /\delta +\mathrm {i} \gamma _n$ (
$n=0,1$ and 
$\gamma _n \in \mathbb {R}^+$). These values represent a solution with local wave decay from one cell to the next (i.e. evanescent waves) although the extension to the infinite periodic array is unphysical.
It can be shown that there is only one real value of 
$\beta$ which exists below one certain frequency and one complex value of 
$\beta ={\rm \pi} /\delta +\mathrm {i} \gamma _1$ that exists above the certain frequency, which means they will not appear at the same frequency. Meanwhile, there exist a number of pure imaginary values of 
$\beta =\mathrm {i}\gamma _0^{(k)}$ 
$(k=1,2,\ldots )$ over the whole frequency range. Thus, we can label different values of 
$\beta$ as 
$\beta = \pm \beta ^{(k)}$ (
$k = 0,1,2,\ldots$), where 
$k=0$ is reserved for an eigenvalue either on the positive real axis or in the complex plane and 
$k=1,2,\ldots$ are used for the pure imaginary values which are ordered with increasing magnitude along the imaginary axes.
The boundary-value problem described above is homogeneous (that is, free of forcing) and we can think of 
$\beta$ as playing the part of the eigenvalue and 
$\phi \neq 0$ the corresponding eigenfunction. Thus, each eigenvalue 
$\pm \beta ^{(k)}$ will be associated with a corresponding eigenfunction which is labelled as 
$\phi = \phi ^{(\pm k)}(x,z)$ such that
\begin{equation} \phi^{(-k)}(x,z) = \phi^{({+}k)}(\delta-x,z), \end{equation}
resulting in 
$\phi ^{(-0)}(x,z)$ being different from 
$\phi ^{(+0)}(x,z)$. In particular, this implies that
\begin{equation} \phi^{(-k)}(0,z) = \exp({\mathrm{i} \beta^{(k)} \delta}) \phi^{({+}k)}(0,z) \end{equation}and
\begin{equation} \phi^{(-k)}_x(0,z) =-\exp({\mathrm{i} \beta^{(k)} \delta}) \phi^{({+}k)}_x (0,z), \end{equation}which will be used extensively later.
2.1. An orthogonality relation
Before we set about solving the eigenproblem, we introduce a useful orthogonality relation. Consider two eigenfunctions 
$\phi ^{(+ k)}(x,z)$ and 
$\phi ^{(\pm j)}(x,z)$ satisfying all of the conditions of the problem described above. Then from Green's identity we have
\begin{align} 0 &= \iint_D \left[ \phi^{(+ k)} \nabla^2 \phi^{(+ j)} - \phi^{(+ j)} \nabla^2 \phi^{(+ k)} \right] \mathrm{d}\kern0.06em x\,\mathrm{d} z \nonumber\\ &= \int_S \left[ \phi^{(+ k)} \frac{\partial \phi^{(+ j)}}{\partial n} - \phi^{(+ j)} \frac{\partial \phi^{(+ k)}}{\partial n} \right] \mathrm{d} s \nonumber\\ &= \left[\exp({\mathrm{i} ( \beta^{(k)} + \beta^{(\,j)})\delta}) - 1 \right] \int_{-h}^{-d}\left[ \phi^{(+ k)} \frac{\partial \phi^{(+ j)}}{\partial x} - \phi^{(+ j)} \frac{\partial \phi^{(+ k)}}{\partial x}\right] _{x=0} \mathrm{d} z \end{align}and
\begin{align} 0 &= \iint_D \left[ \phi^{(+ k)} \nabla^2 \phi^{(- j)} - \phi^{(- j)} \nabla^2 \phi^{(+ k)} \right] \mathrm{d}\kern0.06em x\,\mathrm{d} z \nonumber\\ &= \int_S \left[ \phi^{(+ k)} \frac{\partial \phi^{(- j)}}{\partial n} - \phi^{(- j)} \frac{\partial \phi^{(+ k)}}{\partial n} \right] \mathrm{d} s \nonumber\\ &= \left[\exp({\mathrm{i} ( \beta^{(k)} - \beta^{(\,j)})\delta}) - 1 \right] \int_{-h}^{-d}\left[ \phi^{(+ k)} \frac{\partial \phi^{(- j)}}{\partial x} - \phi^{(- j)} \frac{\partial \phi^{(+ k)}}{\partial x}\right] _{x=0} \mathrm{d} z, \end{align}
after using the conditions on the boundary, 
$S$, of 
$D$ having elemental arclength 
$\mathrm {d} s$ and outward normal derivative 
$\partial /\partial n$. The factor in front of the integral in (2.10) is zero if 
$\beta ^{(k)} =\beta ^{(\,j)}$, while the factor in front of the integral in (2.9) cannot be zero except for 
$\beta ^{(k)} =\beta ^{(\,j)} = {\rm \pi}/\delta$, where 
$k=j=0$. Assuming for now that the eigenvalues 
$\beta ^{(k)}$ are distinct and 
$\beta ^{( 0)} \neq {\rm \pi}/\delta$, it follows the following orthogonality relation:
$$\begin{gather} \int_{-h}^{-d}\left[ \phi^{({+}k)}(0,z) \frac{\partial \phi^{(+ j)}}{\partial x}(0,z) - \phi^{(+ j)}(0,z) \frac{\partial \phi^{(+ k)}}{\partial x}(0,z) \right]\mathrm{d} z =0, \end{gather}$$
$$\begin{gather}\int_{-h}^{-d}\left[ \phi^{({+}k)}(0,z) \frac{\partial \phi^{(- j)}}{\partial x}(0,z) - \phi^{(- j)}(0,z) \frac{\partial \phi^{(+ k)}}{\partial x}(0,z) \right]\mathrm{d} z= E^{(+ k)} \delta_{kj}, \end{gather}$$
where 
$E^{(+ k)}$ is a scaling factor defined by
\begin{align} E^{(+ k)} &= \int_{-h}^{-d}\left[ \phi^{(+ k)}(0,z) \frac{\partial \phi^{(- k)}}{\partial x}(0,z) - \phi^{(- k)}(0,z) \frac{\partial \phi^{(+ k)}}{\partial x}(0,z) \right] \mathrm{d} z \nonumber\\ &=-2 \int_{-h}^{-d} \phi^{(-k)}(0,z) \frac{\partial \phi^{(+ k)}}{\partial x}(0,z)\,\mathrm{d} z=-E^{(- k)}, \end{align} The special case 
$\beta ^{(0)} = {\rm \pi}/\delta$ relates to standing waves in the cell and (2.6) no longer defines an independent second function since 
$\phi ^{(0)}(x,z)=\phi ^{(+ 0)}(x,z)=\phi ^{(- 0)}(x,z)$. Instead, we define 
$\phi ^{(0)}(x,z)$ satisfying (2.5) by imposing supplementary constraints:
\begin{equation} \phi^{(0)}(0,z) =\phi^{(0)}(\delta,z) = 0 \end{equation}or
\begin{equation} \phi^{(0)}_x(0,z) = \phi^{(0)}_x(\delta,z) = 0 , \end{equation}
for 
$-h < z < -d$. It follows that 
$\phi ^{(0)}_x(\delta /2,z) = 0$ or 
$\phi ^{( 0)}(\delta /2,z) = 0$ and the solutions here relate to sloshing modes in closed rectangular domains of width 
$\delta$ either with or without a vertical baffle along the centreline. For the latter case, it can be deduced that it happens at 
$K=(n{\rm \pi} /\delta )\tanh (n{\rm \pi} h/\delta )$ for 
$n=1,2,\ldots$ (see Mei et al. Reference Mei, Stiassnie and Yue2005). Sloshing modes of similar character, but of more complex geometry, were shown to emerge in Porter & Porter (Reference Porter and Porter2003). Besides, it should be noted that the orthogonality relation (2.11) no longer applies for 
$\beta ^{( 0)} = {\rm \pi}/\delta$.
2.2. Solution of two independent forms
There are many different approaches one could adopt to develop solutions to the single-cell problem which partly depend upon how the fundamental cell is defined. Since our choice of the fundamental cell, 
$D$, is rectangular with conditions on the boundary of 
$D$ it makes sense to use separation of variables. In the following, we describe two different forms of solution for the Bloch–Floquet problem given above.
In the first form, the cell is divided into two subdomains which are above and below the level 
$z=-d$, and the solution is expanded by eigenfunctions in 
$x$. In 
$-d < z < 0$, we write the general solution satisfying (2.1), (2.2) and (2.3) as
\begin{equation} \phi(x,z) = a_{1,0} (1+ Kz) + \sum_{n=1}^\infty a_{1,n} \cos (n {\rm \pi}x/\delta) \zeta_n(z), \end{equation}where
\begin{equation} \zeta_n(z) = \frac{\cosh (n {\rm \pi}z/\delta) + (K \delta/n{\rm \pi}) \sinh ( n {\rm \pi}z/\delta)}{\cosh (n {\rm \pi}d/\delta)},\quad n \geq 1 \end{equation}
and 
$a_{1,n}$ for 
$n=0,1,\ldots$ are coefficients to be determined. In 
$-h < z < -d$, the general solution of (2.1) satisfying (2.4) and (2.5) is
\begin{equation} \phi(x,z) = \sum_{n=-\infty}^\infty b_{1,n} \frac{\cosh \beta_n(h+z)} {\cosh \beta_n(h-d)} \exp({\mathrm{i} \beta_n x}), \end{equation}where
\begin{equation} \beta_n = \beta + 2 n {\rm \pi}/\delta \end{equation}
and 
$b_{1,n}$ for 
$n \in \mathbb {Z}$ are also undetermined coefficients.
The pressure and vertical component of velocity must coincide across the common fluid interface 
$z=-d$ for 
$0 < x < \delta$. We first define
\begin{equation} w(x) = \phi_z(x,-d), \quad 0 < x < \delta, \end{equation}
which represents the vertical velocity across 
$z=-d$. From (2.15) and the orthogonality of the cosine functions over 
$0 < x < \delta$, continuity of velocity allows us to write
\begin{equation} a_{1,n} = \left\{\begin{array}{@{}ll} \dfrac{1}{K \delta} \displaystyle \int_{0}^\delta w(x)\,\mathrm{d}\kern0.06em x, & n=0, \\ \dfrac{2}{n{\rm \pi}[K\delta-n{\rm \pi} \tanh (n{\rm \pi} d/\delta)]} \int_{0}^\delta w(x) \cos(n {\rm \pi}x/ \delta)\,\mathrm{d}\kern0.06em x, & n \geq 1. \end{array}\right. \end{equation}
Also, from (2.17) we have, using orthogonality of the functions 
$\textrm {e}^{\mathrm {i} \beta _n x}$ over 
$0 < x < \delta$,
\begin{equation} b_{1,n} = \frac{1}{\beta_n \delta \tanh\beta_n(h-d)} \int_0^\delta w(x) \exp({-\mathrm{i} \beta_n x})\,\mathrm{d}\kern0.06em x,\quad n \in \mathbb{Z}. \end{equation}
We now match pressure at 
$z = -d$ and substitute (2.19) and (2.20) for 
$a_{1,n}$ and 
$b_{1,n}$ to give the scalar homogeneous integral equation
\begin{equation} \int_0^\delta w(x') {\mathcal{L}}(x,x')\,\mathrm{d}\kern0.06em x' = 0,\quad 0 < x < \delta, \end{equation}where
\begin{align} {\mathcal{L}}(x,x') = \frac{Kd-1}{K\delta} + 2 \sum_{n=1}^\infty D_n \cos (n {\rm \pi}x/\delta) \cos (n {\rm \pi}x'/\delta)+ \sum_{n=-\infty}^\infty \frac{\exp({\mathrm{i} \beta_n (x-x')})}{\beta_n \delta \tanh \beta_n(h-d)} \end{align}and
\begin{equation} D_n =\left(\frac{1}{n {\rm \pi}} \right) \frac{(K\delta/n {\rm \pi}) \tanh (n {\rm \pi}d/\delta)-1}{(K \delta/n{\rm \pi})- \tanh (n {\rm \pi}d/\delta)} \sim \frac{1}{n {\rm \pi}}, \end{equation}
as 
$n \to \infty$ (or as 
$\delta /d \to 0$).
The second form is to use eigenfunctions in the depth coordinate to expand the solution which is a natural approach to solving water wave problems (Linton & McIver Reference Linton and McIver2001). This leads to the solution being posed in terms of integral equations over finite intervals of 
$x=0$: either for the unknown 
$\phi _x(0,z)$ between 
$-h < z < -d$ or for the unknown 
$\phi (\delta ^-,z) - \textrm {e}^{\mathrm {i} \beta \delta } \phi (0^+,z)$ between 
$-d < z < 0$. Given the relation derived in § 2.1, the first of these two options is particularly attractive. Thus, the velocity potential is first written as
\begin{equation} \phi(x,z)=\sum_{n=0}^\infty \left( a_{2,n} \,{\rm e}^{k_n x} +b_{2,n} \exp({-k_n x})\right ) \psi_n(z). \end{equation}
Here 
$k_n$ are the roots of the dispersion equation
\begin{equation} \omega^2/g =-k_n \tan k_nh, \end{equation}
where 
$k_n$ 
$(n\geq 1)$ is real and positive while 
$k_0 = -\mathrm {i} k$ and 
$k$ is the real positive wavenumber, and
\begin{equation} \psi_n(z) = N_n^{-1/2} \frac{\cos k_n (z+h)}{\cos k_n h}, \end{equation}with
\begin{equation} N_n = \frac{1}{2\cos^2 k_n h} \left( 1+ \frac{\sin 2 k_n h}{2 k_n h} \right), \end{equation}which satisfy the orthogonality relation
\begin{equation} \frac{1}{h} \int_{-h}^0 \psi_n(z) \psi_m(z)\,\mathrm{d} z = \delta_{mn}. \end{equation}We now define
\begin{equation} u(z) = \phi_x(0,z), \quad -h < z <-d, \end{equation}
which represents the horizontal velocity across 
$x=0$. From the velocity periodic condition in (2.5) and the orthogonality relation in (2.29), we have
\begin{equation} a_{2,n}=\frac{{\rm e}^{\mathrm{i} \beta \delta}-\exp({-k_n \delta})}{2 k_n h\sinh k_n \delta} \int_{-h}^{-d} u(z) \psi_n (z)\,\mathrm{d} z \end{equation}and
\begin{equation} b_{2,n}=\frac{{\rm e}^{\mathrm{i} \beta \delta}-{\rm e}^{k_n \delta}}{2 k_n h\sinh k_n \delta } \int_{-h}^{-d} u(z) \psi_n (z)\,\mathrm{d} z.\end{equation}Applying the pressure periodic condition in (2.5) with (2.31) and (2.32) results in another scalar homogeneous integral equation:
\begin{equation} \int_{-h}^{-d} u(z') {\mathcal{K}}(z,z')\,\mathrm{d} z' = 0,\quad -h < z <-d, \end{equation}where
\begin{equation} {\mathcal{K}}(z,z')=\sum_{n=0}^\infty \frac{\cos \beta \delta-\cosh k_n \delta}{k_n h \sinh k_n \delta}\psi_n (z) \psi_n (z'). \end{equation}2.3. Numerical approximation
The numerical approximation of the integral equations is based on methods described in Porter & Evans (Reference Porter and Evans1995) in which it is recognised that the end points, 
$(0,-d)$ and 
$(\delta,-d)$, of the intervals involved in the integral equations coincide with the sharp edges of the barriers and the fluid velocity behaves as the inverse square root of distance to the edge. Thus, for the first form given in the previous section we choose
\begin{equation} w(x) \approx \sum_{m=0}^{M_1} \alpha_{1,m} w_{m}(x), \end{equation}where
\begin{equation} w_m(x) = \frac{T_m(2 x/\delta - 1)}{{\rm \pi} \sqrt{ (\delta/2)^2 -(x-\delta/2)^2}}, \end{equation}
in which 
$T_m({\cdot })$ is a Chebyshev polynomial and 
$M_1$ is the designated truncation parameter. In what follows we use the results that follow from, for example, Erdélyi et al. (Reference Erdélyi, Magnus, Oberhettinger and Tricomi1954):
\begin{equation} \int_{0}^\delta w_m(x) \cos (n {\rm \pi}x/\delta) \, \mathrm{d}\kern0.06em x = \cos [(m+n){\rm \pi}/2] {\rm J}_m(n {\rm \pi}/2) \end{equation}and
\begin{equation} \int_{0}^\delta w_m(x) \exp({-\mathrm{i} \beta_n x})\,\mathrm{d}\kern0.06em x = \exp({-\mathrm{i} \beta_n \delta/2}) \exp({-\mathrm{i} m {\rm \pi}/2}) {\rm J}_m(\beta_n \delta/2), \end{equation}
where 
$\textrm {J}_m({\cdot })$ are the 
$m$th-order Bessel functions of the first kind. Substituting the approximation (2.35) into (2.22), multiplying through by 
$w_n(x)$ and integrating over 
$0 < x < \delta$ result in the following system of equations for the expansion coefficients 
$\alpha _{1,n}$:
\begin{equation} \sum_{n=0}^{M_1} \alpha_{1,n} L_{mn} = 0,\quad m = 0,1,\ldots,M_1, \end{equation}where
\begin{align} L_{mn} &= \frac{Kd-1}{K \delta} \delta_{m0} \delta_{n0} + 2 \sum_{r=1}^\infty D_r \cos\left[{\frac12} (m+r){\rm \pi}\right] \nonumber\\ &\quad \times\cos\left[{\frac12} (n+r){\rm \pi}\right] {\rm J}_m\left({\frac12} r {\rm \pi}\right) {\rm J}_n\left({\frac12} r {\rm \pi}\right) \nonumber\\ &\quad + \exp({-\mathrm{i} (m-n){\rm \pi}/2}) \sum_{r=-\infty}^\infty \frac{{\rm J}_m\left(\dfrac12 \beta_r \delta\right) {\rm J}_n\left(\dfrac12 \beta_r \delta\right)}{\beta_r \delta \tanh \beta_r (h-d)}. \end{align} Similarly, for the second form we expand the unknown horizontal velocity 
$u(z)$ as
\begin{equation} u(z) \approx \sum_{m=0}^{M_2} \alpha_{2,m} u_{m}(z) \end{equation}
in a series of 
$M_2+1$ prescribed functions:
\begin{equation} u_m(z) = \frac{2(-1)^m T_{2m}[(h+z)/(h-d)]}{{\rm \pi} \sqrt{ (h-d)^2 -(h+z)^2}}, \end{equation}
which satisfy the condition (2.4) on the sea bed. Substituting the approximation (2.41) into (2.33), multiplying through by 
$u_{n}(z)$ and integrating over 
$-h < z < -d$ lead to the following system of equations for the expansion coefficients 
$\alpha _{2,n}$:
\begin{equation} \sum_{n=0}^{M_2} \alpha_{2,n} K_{mn} = 0, \quad m = 0,1,\ldots,M_2, \end{equation}where
\begin{equation} K_{mn} = \sum_{r=0}^\infty \frac{(\cos \beta \delta-\cosh k_r \delta)}{k_r h \sinh k_r \delta} F_{mr} F_{nr}, \end{equation}in which we have defined
\begin{equation} F_{mr} = \int_{-h}^{-d} u_m(z) \psi_r(z)\,\mathrm{d} z = N_r^{-1/2} {\rm J}_{2m}[k_r(h-d)]. \end{equation} Numerically we fix 
$Kd$ and look for real values of 
$\beta \in (0,{\rm \pi} /\delta ]$, 
$\gamma _1 >0$ with 
$\beta ={\rm \pi} /\delta +\mathrm {i} \gamma$ and 
$\gamma _0^{(k)} >0$ with 
$\beta =\mathrm {i} \gamma _0^{(k)}$ for which the system of (2.39) or (2.43) has non-trivial solutions. When 
$\beta$ is real the matrix formed by 
$L_{mn}$ is Hermitian and when 
$\beta$ is complex the matrix is real, while the matrix formed by 
$K_{mn}$ is a real matrix no matter whether 
$\beta$ is real or complex. These result in that the determinants of the two matrices are always real, so zero eigenvalue of the corresponding matrix can be found numerically using standard root-finding methods. In particular, when 
$\beta$ becomes a pure imaginary number, from the last term in (2.40) we can see that when 
$\gamma _0 (h-d)$ passes across 
$k {\rm \pi}$ each element in 
$L_{mn}$ will tend to positive or negative infinity. Thus, it can be deduced that the overall behaviour of 
$L_{mn}$ is similar to that of the function 
$\cot \gamma _0 (h-d)$ and 
$\gamma _0^{(k)} (h-d) \in ((k-1){\rm \pi},k{\rm \pi} )$ though the same conclusion is not easily drawn from 
$K_{mn}$. Furthermore, all of the series in 
$L_{mn}$ and 
$K_{mn}$ are convergent as of 
$O(1/r^2)$. In order to accelerate the numerical computation, an effective treatment is performed for these series (see Appendix A for details).
As is shown later, either of the two methods of solution presented above can be used to find accurate values of eigenvalues, 
$\beta$. However, accurately determining the eigenfunctions 
$\phi (x,z)$ is more problematic. Even though both methods use functions which accurately capture the inverse square root singularity at the lower edges of the barriers as part of the solution method, the expressions for 
$\phi$ are formed by separation solutions which do not explicitly include these singularities. Consequently, the expansion coefficients associated with the separation solutions are slowly convergent for both methods (like 
$O(1/n^{3/2})$). In order to produce plots of the eigenfunctions so that they may be compared with the results of homogenisation (in § 3.1) we nevertheless find that the first method works well. This is because when 
$\delta /d$ is small, which is the primary interest of the present study, only one or two terms are required in the expansion of the unknown velocity across the level 
$z=-d$ to obtain very accurate solutions whilst these separation solutions are well suited to close spacing with the first term in the expansion above and below 
$z=-d$ being dominant. This is not true, however, for the second method since evanescent modes play an important role when 
$\delta /d$ is small. Thus it is hard to plot eigenfunctions across the whole domain accurately by using the second method for small 
$\delta /d$.
Also later, when we consider scattering by finite arrays using a discrete barrier method we are required to compute integrals over the intervals 
$-h < z < -d$ beneath the barriers involving the eigenfunctions and their derivative with respect to 
$x$. Although separation solutions for determining 
$\phi$ from the first method can do this accurately, the series in which derivatives are taken term-by-term are no longer convergent, with terms decreasing like 
$O(1/n^{1/2})$. For this reason, the second method is useful since the solution method provides highly accurate representations for the derivative which explicitly include the singularity at 
$z=-d$ as shown in (2.30) and (2.41).
2.4. Small 
$\delta /d$ and $\beta \delta$
We now assume that 
$\epsilon =\delta /d \ll 1$ and 
$\beta \delta \ll 1$, that is to say, 
$\beta d \ll 1/\epsilon$. Based on the solution in the first form presented in the previous section, the dominant entry in (2.40) is
\begin{equation} L_{00} \approx \frac{Kd-1}{K\delta} + \frac{1}{\beta \delta \tanh \beta (h-d)}, \end{equation}
after using 
$\textrm {J}_m(\beta \delta /2) \approx \delta _{m0}$ and assuming that 
$K\delta /|Kd-1|$ has the same order as 
$\beta \delta$, i.e. 
$K\delta /|Kd-1| \ll 1$. Thus, the leading-order approximation to values of 
$\beta$ for small 
$\delta /d$ is determined from solving 
$L_{00} = 0$, or
\begin{equation} \beta \tanh \beta(h-d) = \frac{K}{1-Kd}, \end{equation}
provided 
$|1-Kd|/Kd \gg \epsilon$ which implies that 
$|1-Kd|\gg \epsilon$. As for the velocity potential, if we normalise 
$\alpha _n$ by setting 
$\alpha _0=K\delta /(1-Kd)$, after using (2.47) the velocity potential can be written as
\begin{equation} \phi(x,z) \approx \left\{\begin{array}{@{}ll} (1+Kz)/(1-Kd), & -d< z<0, \\ {\rm e}^{\mathrm{i} \beta x}\cosh \beta(h+z)/\cosh \beta(h-d), & -h< z<-d. \end{array}\right. \end{equation}
From (2.48), we can see that when 
$\delta / d \ll 1$ and 
$\beta \delta \ll 1$, the expressions only including the first term 
$(n=0)$ in (2.15) and (2.17) are good at simulating the flat oscillation in the cell. Besides, (2.47) is similar to the dispersion equation (2.26) with 
$n=0$. It can be proved that here 
$\beta$ only can be real or pure imaginary but could not be a complex number and the real value exists only when 
$Kd<1$.
3. A continuum model
In this section, we develop an approximation to wave propagation through the infinite periodic array by directly applying asymptotic methods to the underlying boundary-value problem. The principal assumption is that 
$\epsilon = \delta /d \ll 1$ (close spacing between adjacent barriers), which is the same as in § 2.4. In 
$-d< z<0$, we make a multiple scales approximation, i.e. 
$x \to d \hat {x} + \delta X$, where 
$\hat {x}$ is the macroscale variable and 
$X$ operates on the scale of a single cell. We also scale 
$z \to d \hat {z}$ and write
\begin{equation} \phi(x,z) \approx \phi^{(0)}(\hat{x},X,\hat{z}) + \epsilon \phi^{(1)}(\hat{x},X,\hat{z}) + \epsilon^2 \phi^{(2)}(\hat{x},X,\hat{z}) + \cdots. \end{equation} In 
$0< X<1$ and 
$-1 < \hat {z} < 0$, from (2.1) we have
\begin{equation} \left[\frac{\partial^2 }{\partial X^2} + 2\epsilon \frac{\partial}{\partial \hat{x} \partial X}+ \epsilon^2 \left(\frac{\partial^2 }{\partial \hat{x}^2}+\frac{\partial^2 }{\partial \hat{z}^2} \right) \right] \left(\phi^{(0)}+ \epsilon \phi^{(1)}+ \epsilon^2 \phi^{(2)} + \cdots \right) = 0, \end{equation}with
\begin{equation} \left(\frac{\partial}{\partial X}+\epsilon \frac{\partial}{\partial \hat{x}} \right) \left( \phi^{(0)}+ \epsilon \phi^{(1)} + \epsilon^2 \phi^{(2)} + \cdots \right) = 0,\quad \mbox{on } X=0,1 \end{equation}and
\begin{equation} \left( \frac{\partial}{\partial \hat{z}} - Kd \right) \left( \phi^{(0)}+ \epsilon \phi^{(1)} + \epsilon^2 \phi^{(2)} + \cdots \right) = 0,\quad \mbox{on } \hat{z} = 0. \end{equation}Using (3.2) with (3.3) for the zero order gives
\begin{equation} \phi^{(0)}(\hat{x},X,\hat{z}) \equiv \phi^{(0)}(\hat{x},\hat{z}), \end{equation}
which is independent of 
$X$. At the first order, (3.2) is
\begin{equation} \frac{\partial^2 \phi^{(1)}}{\partial X^2} =- 2\frac{\partial^2 \phi^{(0)}}{\partial \hat{x} \partial X} = 0, \end{equation}
where (3.5) has been applied. Integrating (3.6) over 
$0 < X < 1$ and using the boundary conditions implied by (3.3) for 
$\phi ^{(1)}$ on 
$X=0,1$ give
\begin{equation} \phi^{(1)} =-\frac{\partial \phi^{(0)}}{\partial \hat{x}}X+f(\hat{x},\hat{z}), \end{equation}
where 
$f(\hat {x},\hat {z})$ is a function independent of microscale variable 
$X$. Further, if we consider the second-order term in (3.2) and (3.3), we can obtain
\begin{equation} \frac{\partial^2 }{\partial X^2}\phi^{(2)} + 2\frac{\partial^2 }{\partial \hat{x} \partial X} \phi^{(1)}+ \left(\frac{\partial^2}{\partial \hat{x}^2} +\frac{\partial^2}{\partial \hat{z}^2} \right) \phi^{(0)}= 0, \end{equation}with
\begin{equation} \frac{\partial \phi^{(2)}}{\partial X}+\frac{\partial \phi^{(1)}}{\partial \hat{x}}=0,\quad \mbox{on } X=0,1. \end{equation}
Again, integrating (3.8) over 
$0< X<1$ and applying (3.5), (3.7) and (3.9) result in
\begin{equation} \frac{\partial^2 \phi^{(0)}}{\partial \hat{z}^2}=0 \end{equation}
as the leading-order governing equation, with (3.4) applying at zeroth order. The general solution, satisfying (3.10), is 
$\phi ^{(0)} = X(\hat {x}) (1+Kd \hat {z})$ for an arbitrary function 
$X(\hat {x})$. Since we are concerned with wave propagation, we write
\begin{equation} \phi^{(0)}(\hat{x},\hat{z}) = A \exp({\mathrm{i} \mu d \hat{x}}) (1+Kd \hat{z}), \end{equation}
where 
$\mu$ is a coefficient to be determined and the assumption is that 
$\mu d \ll 1/\epsilon$ otherwise the horizontal variation is not on the macroscale.
In 
$-h < z < -d$, since the fluid is not bounded by barriers, we can drop the microscale, resulting in that we rescale with 
$x \to d \hat {x}$ and 
$z \to d \hat {z}$. After expanding the velocity potential with respect to 
$\epsilon$, we can find that 
$\phi ^{(0)}(\hat {x},\hat {z})$ still satisfies Laplace's equation:
\begin{equation} \left( \frac{\partial^2 }{\partial \hat{x}^2} + \frac{\partial^2 }{\partial \hat{z}^2} \right) \phi^{(0)} = 0. \end{equation}After applying the separation of variables, the solution of (3.12) satisfying the zeroth-order condition on the sea bed is
\begin{equation} \phi^{(0)}(x,z) = B \exp({\mathrm{i} \mu' d \hat{x}}) \cosh \mu' d(\hat{z}+h/d), \end{equation}
where 
$\mu '$ is also an undetermined coefficient.
Applying the continuity of pressure and vertical component of velocity for (3.11) and (3.13) on the common fluid interface 
$\hat {z}=-1$ results in
$$\begin{gather} \mu=\mu', \end{gather}$$
$$\begin{gather}B \cosh \mu (h-d) = A(1 -Kd) \end{gather}$$and
\begin{equation} B \mu \sinh \mu(h-d) = A K. \end{equation}
That is, 
$\mu$ satisfies
\begin{equation} \mu \tanh \mu (h-d) = K/(1-Kd), \end{equation}and the corresponding mode, written in terms of the original coordinates, is
\begin{equation} \phi(x,z) = {\rm e}^{\mathrm{i} \mu x} \left\{ \begin{array}{@{}ll} (1+Kz)/(1-Kd), & -d < z < 0, \\ \cosh \mu(h+z)/\cosh \mu(h-d), & -h < z <-d. \end{array}\right. \end{equation}
When 
$Kd \to 1$ we can see from (3.17) that 
$\mu d \to \infty$ and thus the assumption made in (3.11) that 
$\mu d \ll 1/\epsilon$ is violated.
In this section we have implemented a ‘low-frequency homogenisation’ and it can be expected to be valid if 
$|1-Kd| \gg \epsilon$ which is also aligned with the assumption made in § 2.4. The result (3.17) and (3.18) shows that the homogenisation of the boundary-value problem coincides with the small-
$\delta /d$ limit (2.47) and (2.48) of the discrete barrier array problem with the association that 
$\beta \to \mu$ as 
$\delta /d \to 0$. Also, it can be proved that there exists one real value of 
$\mu _0$ and a number of pure imaginary values of 
$\mu _n=\mathrm {i} \hat {\mu }_n$ (where 
$n=1,2,\ldots$ and 
$\hat {\mu }_n \in ((n-1){\rm \pi},n{\rm \pi} )$) satisfying (3.17).
It should be possible to perform a ‘high-frequency homogenisation’ (see Craster, Kaplunov & Pichugin Reference Craster, Kaplunov and Pichugin2010) by expanding about the state 
$\beta \delta = {\rm \pi}$ where a local standing mode exists and gives rise to the critical value of 
$K_c$ below which wave propagation exists and above which wave propagation is prohibited.
3.1. Results
First, we determine the accuracy of the numerical scheme of § 2 by varying the truncation parameter, 
$M_{k}$ 
$(k=1,2)$, to assess the convergence of the two schemes given in § 2.2. Fundamental cells with the same submergence 
$d/h=0.2$ but different spacings are considered in tables 1 and 2 which catalogue the numerical estimates of the first 
$\beta ^{(0)}$ and second 
$\beta ^{(1)}$ Bloch wavenumbers by solving each scheme. For the first Bloch wavenumber, when 
$K$ exceeds a value of 
$K_c$ corresponding to the critical frequency 
$\omega _c$ (where 
$K_c d < 1$), there no longer exists a real-valued Bloch wavenumber. Instead, following the system introduced in § 2, 
$\beta ^{(0)}$ records a complex Bloch wavenumber and represents decay rather than wave propagation through the array. For the real value of 
$\beta ^{(0)}$, when 
$M_{k}=6$ the non-dimensional wavenumber 
$\beta ^{(0)}\delta$ is determined to have nearly five-decimal-place accuracy except for the second scheme at relatively high frequencies. When the frequency exceeds the critical frequency, the larger truncation parameters are required for obtaining the first Bloch wavenumber with the same precision. As for the second Bloch wavenumber, the convergent results can be reached with very few terms. Generally, the first scheme tends to converge faster for small 
$\delta /d$ and the second scheme does better when 
$\delta /d$ takes larger values, the reason for which has been outlined earlier.
The convergence of first non-dimensional Bloch wavenumber 
$\beta ^{(0)}\delta$ against the truncation parameter, 
$M_{k}$, for the two schemes given in § 2.2 with 
$d/h=0.2$. A dash (—) indicates cannot determine a value of 
$\beta$ that makes the determinant zero.
The convergence of second non-dimensional Bloch wavenumber 
$\beta ^{(1)}\delta$ against the truncation parameter, 
$M_{k}$, for the two schemes given in § 2.2 with 
$d/h=0.2$. A dash (—) indicates cannot determine a value of 
$\beta$ that makes the determinant zero.
Next, we compare Bloch wavenumber 
$\beta ^{(k)}$ obtained by (2.39) or (2.43) with numerical roots 
$\mu _k$ in the homogenisation method obtained by (3.17). Figure 1(a) shows the variation of the first value 
$(k=0)$ against the non-dimensional wavenumber 
$Kd$ for submergence 
$d/h=0.2$. As mentioned in § 2, real 
$\beta$ are determined in the range of 
$\beta \in (0,{\rm \pi} /\delta ]$ and for the other ranges the problem is unchanged. Thus, in figure 1(a) the curves describing real Bloch wavenumbers terminate at 
$\beta ={\rm \pi} /\delta$ and the value of 
$K_cd$ corresponding to the critical frequencies in the figure shown are, respectively, 
$0.9891$, 
$0.9477$ and 
$0.9006$. Above these critical frequencies, a complex value of 
$\beta ={\rm \pi} /\delta +\mathrm {i} \gamma _1$ emerges from the real axis, the imaginary part of which is also shown in figure 1(a). Thus, it can be inferred that when the frequency exceeds the critical frequency, the real Bloch–Floquet wavenumber will move off the real axis and go along the semi-infinite line 
$\beta ={\rm \pi} /\delta +\mathrm {i} \gamma _1$ for 
$\gamma _1 = [0,\infty )$. Besides, we can see that for small spacing, the complex solution increases extremely fast with frequency. Since the propagating mode no longer exists, we also can conclude that the first stop band during which wave propagation is prohibited is the interval 
$Kd \in (K_{c} d, ({\rm \pi} d/\delta )\tanh ({\rm \pi} h/\delta ))$. The end points of this interval both correspond to 
$\beta \delta = {\rm \pi}$ and standing waves occurring in the cell; the lower value corresponds to the case which is equivalent to a vertical baffle placed along the centreline of the cell, i.e. the solution of (2.13), while the upper value corresponds the case in which barriers extend through the depth, i.e. the solution of (2.14). The real root, 
$\mu _0$, of the dispersion equation (3.17) tends to infinity as 
$Kd \to 1$. The stop band under the homogenisation approximation is 
$Kd \in (1, \infty )$, i.e. the critical frequency 
$\omega _c=\sqrt {g/d}$ coinciding with Newman (Reference Newman1974) and representing resonance in narrow channels. Figure 1(b) plots the variation of the first five pure imaginary values. Generally, as the dimensionless spacing, 
$\delta /d$, decreases 
$\beta$ tends to 
$\mu$ as we have anticipated.
The variation of the roots 
$\mu _k$ from homogenisation and the Bloch–Floquet wavenumbers 
$\beta ^{(k)}$ against the non-dimensional wavenumber 
$Kd$ for a submergence 
$d/h=0.2$: (a) the first value 
$(k=0)$; (b) the imaginary part of the first five pure imaginary values.
In figure 2, the velocity potential fields have been plotted in the range 
$0 < x/h < 0.3$ for a barrier submergence 
$d/h=0.2$ at a non-dimensional wavenumber 
$Kd=0.8$. In figure 2(a–c), the results correspond to 
$\beta ^{(0)}$ and the barrier spacings are reduced from 
$\delta /d = 0.5$, to 
$0.25$ and then to 
$0.05$ so that we see 
$3$, 
$6$ and 
$30$ cells respectively. It can be seen that oscillation within each channel is dominated by vertical fluid motion and as the channels decrease in width, the results tend towards the potential field obtained under homogenisation, shown in figure 2(d). Only the real part of the velocity potential is shown, but there is a similar agreement for the imaginary part.
The fields of the real part of the velocity potential in the cell with the submergence 
$d/h=0.2$ at 
$Kd=0.8$: (a) 
$\delta /d=0.50$; (b) 
$\delta /d=0.25$; (c) 
$\delta /d=0.05$; (d) homogenisation.
In figure 3, we plot the fields of the velocity potential at 
$Kd=0.9891$ where 
$\beta ^{(0)} \delta ={\rm \pi}$ for the case of 
$\delta /d=0.05$ (the smallest spacing used in the previous plot). This is the case in which standing waves occur in the cell and the imaginary part of the velocity potential vanishes according to the Bloch–Floquet theory. The potentials are not normalised giving rise to large values in the plots. Unlike in figure 2, there is no longer good agreement between the Bloch–Floquet approach shown in figure 3(a) and homogenisation (the real and imaginary parts of which are shown in figures 3b and 3c). This illustrates how homogenisation breaks down as an approximation to closely spaced discrete arrays as standing wave resonance is approached.
Velocity potential field plots across 
$0 < x < 0.03$ with the submergence 
$d/h=0.2$ at 
$Kd=0.9891$: (a) 
$\delta /d=0.05$ and 
$\beta ^{(0)} \delta = {\rm \pi}$; (b) homogenisation (real part); (c) homogenisation (imaginary part).
4. Scattering of incident waves by a finite periodic array of barriers
In this section, we consider the problem of 
$N+1$ identical barriers each submerged to the same depth 
$d$ and located at 
$x_n = n \delta$ for 
$n=0,1,\ldots,N$ as shown in figure 4, which is a finite section of the assumed infinite array given in § 2. Thus, the general solution in 
$(n-1) \delta < x < n \delta$ can be expressed as a combination of the eigenfunctions 
$\phi ^{(\pm k)}$ of the periodic Bloch–Floquet problem associated with 
$\beta = \pm \beta ^{(k)}$:
\begin{equation} \phi_n(x,z) = \sum_{k=0}^\infty \left[ c_n^{(k)} \phi^{({+}k)}(x-(n-1)\delta,z) + d_n^{(k)} \phi^{(-k)}(x-(n-1)\delta,z) \right], \end{equation}
for 
$n=1,2,\ldots,N$. This representation of the solution was established in Porter & Porter (Reference Porter and Porter2003). It should be noted that since the local wave with a slow decay should be given priority for considering the oscillation in the barrier array, in this section a new labelling rule for 
$\beta ^{(k)}$ is applied that 
$k=0$ is prepared for the real eigenvalue if exists otherwise the complex eigenvalues of 
$\beta =n {\rm \pi}/\delta +\mathrm {i} \gamma _n$ 
$(n=0,1)$ are ordered with their imaginary parts increasing, which is different from the labelling rule in § 2.
An illustration of scattering by a finite periodic array of 
$N+1$ identical surface-piercing barriers.
Continuity of pressure and velocity across the fluid interface under the barrier along 
$x=n\delta$ requires
\begin{equation} \phi_{n+1}(n \delta,z) = \phi_n(n\delta,z)\quad \mbox{and} \quad \frac{\partial \phi_{n+1}}{\partial x}(n \delta,z) = \frac{\partial \phi_n}{\partial x}(n \delta,z), \end{equation}
for 
$-h < z < -d$ and 
$n=1,2,\ldots,N-1$. Using the orthogonality relation (2.11) derived earlier in § 2.1, we find that
\begin{equation} c_{n+1}^{(\,j)} E^{({+}j)} = \int_{-h}^d \left[ \phi_{n+1}(n\delta,z) \frac{\partial \phi^{(-j)}}{\partial x}(0,z) - \phi^{(-j)}(0,z) \frac{\partial \phi_{n+1}}{\partial x}(n\delta,z) \right] \mathrm{d} z \end{equation}and
\begin{equation} d_{n+1}^{(\,j)} E^{(-j)} = \int_{-h}^d \left[ \phi_{n+1}(n\delta,z) \frac{\partial \phi^{({+}j)}}{\partial x}(0,z) - \phi^{({+}j)}(0,z) \frac{\partial \phi_{n+1}}{\partial x}(n\delta,z)\right] \mathrm{d} z. \end{equation}Using the matching conditions (4.2a,b) and the definition (4.3) gives
\begin{align} c_{n+1}^{(\,j)} E^{({+}j)} &= \int_{-h}^d \left[\phi_{n}(n \delta,z) \frac{\partial \phi^{(-j)}}{\partial x}(0,z) - \phi^{(-j)}(0,z) \frac{\partial \phi_n}{\partial x}(n \delta,z)\right]_{x=0} \mathrm{d} z \nonumber\\ &= c_{n}^{(\,j)} E^{({+}j)} \exp({+\mathrm{i} \beta^{(\,j)} \delta}), \end{align}
after using the phase relations (2.5) for the 
$j$th eigenfunction to transfer information from 
$x=\delta$ to 
$x=0$ and the orthogonality relation (2.11) again. We do the same for 
$d_n^{(\,j)}$ allowing us to deduce that
\begin{equation} c_{n+1}^{(\,j)} = \exp({+\mathrm{i} \beta^{(\,j)} \delta}) c_{n}^{(\,j)}\quad \mbox{and} \quad d_{n+1}^{(\,j)} = \exp({-\mathrm{i} \beta^{(\,j)} \delta}) d_{n}^{(\,j)}. \end{equation}That is, there is no coupling between eigenmodes as waves propagate through the periodic array having the consequence that
\begin{equation} c_{N}^{(\,j)} = \exp({\mathrm{i} \beta^{(\,j)} (N-1) \delta}) c_{1}^{(\,j)} \quad \mbox{and} \quad d_{N}^{(\,j)} = \exp({-\mathrm{i} \beta^{(\,j)} (N-1) \delta}) d_{1}^{(\,j)}. \end{equation}
In other words, if the solution in 
$0 < x < \delta$ is expressed as
\begin{equation} \phi_1(x,z) = \sum_{k=0}^\infty \left[ c_1^{(k)} \phi^{({+}k)}(x,z) + d_1^{(k)} \phi^{(-k)}(x,z) \right], \end{equation}
then the general solution across the whole barrier array domain 
$0 < x < N \delta$ is represented in terms of just one set of expansion coefficients, 
$c_1^{(k)}$ and 
$d_1^{(k)}$. In particular, the solution in 
$(N-1) \delta < x < N \delta$ is
\begin{align} \phi_N(x,z) &= \sum_{k=0}^\infty [c_1^{(k)} \exp({\mathrm{i} \beta^{(k)} (N-1) \delta}) \phi^{({+}k)}(x-(N-1)\delta,z) \nonumber\\ &\quad +d_1^{(k)} \exp({-\mathrm{i} \beta^{(k)} (N-1) \delta}) \phi^{(-k)}(x-(N-1)\delta,z)]. \end{align} The scattering problem involves waves incident from and reflected into the domain 
$\{x < 0$, 
$-h < z < 0\}$ in which the general solution is represented by the standard expansion (e.g. Linton & McIver Reference Linton and McIver2001)
\begin{equation} \phi_0(x,z) = ({\rm e}^{\mathrm{i} k x} + R_N \exp({-\mathrm{i} k x})) \psi_0(z) + \sum_{n=1}^\infty a_n \,{\rm e}^{k_n x} \psi_n(z), \end{equation}
where 
$R_N$ is the reflection coefficient, 
$a_n$ (and, later, 
$b_n$) are expansion coefficients and 
$\psi _n (z)$ is the vertical eigenfunction which has been given in (2.27). It helps to write (4.10) as
\begin{equation} \phi_0(x,z) = 2 \cos kx \psi_0(z) + \sum_{n=0}^\infty a_n \,{\rm e}^{k_n x} \psi_n(z), \end{equation}
where 
$R_N = 1+ a_0$. In 
$x > N \delta$, waves are transmitted and the general solution is represented by
\begin{equation} \phi_{N+1}(x,z) = \sum_{n=0}^\infty b_n \exp({-k_n (x-N \delta)}) \psi_n(z), \end{equation}
with transmission coefficient 
$T_N = b_0\exp ({-\mathrm {i} k N \delta })$.
With (4.8) holding in the region adjoining 
$x=0$ and (4.9) holding in the region adjoining 
$x=N\delta$, the remaining conditions that need to be enforced in order to determine the values of 
$a_n$, 
$b_n$ for 
$n=0,1,\ldots$ and 
$c_1^{(k)}$, 
$d_1^{(k)}$ for 
$k = 0,1,\ldots$ are
\begin{equation} \frac{\partial \phi_0}{\partial x}(0,z) = 0 \quad \mbox{and}\quad \frac{\partial \phi_{N+1}}{\partial x}(N\delta,z) = 0, \end{equation}
for 
$-d < z < 0$, and
$$\begin{gather} \phi_0(0,z) = \phi_1(0,z) \quad \mbox{and} \quad \phi_N(N \delta,z) = \phi_{N+1}(N \delta,z), \end{gather}$$
$$\begin{gather}\frac{\partial \phi_0}{\partial x}(0,z) = \frac{\partial \phi_1}{\partial x}(0,z) \equiv U(z) \quad \mbox{and} \quad \frac{\partial \phi_N}{\partial x}(N \delta,z) = \frac{\partial \phi_{N+1}}{\partial x}(N \delta,z) \equiv V(z), \end{gather}$$
for 
$-h < z < -d$. Applying the condition (4.15a,b) to (4.11) and (4.12) and the orthogonality of the vertical eigenfunctions (2.29) results in
\begin{equation} \phi_0(x,z) = 2 \cos kx \psi_0(z) + \sum_{n=0}^\infty \frac{\psi_n(z) \,{\rm e}^{k_n x}}{k_n h} \int_{-h}^{-d} U(z') \psi_n(z')\,\mathrm{d} z' \end{equation}
in 
$x < 0$ and
\begin{equation} \phi_{N+1}(x,z) =-\sum_{n=0}^\infty \frac{\psi_n(z) \exp({-k_n (x-N \delta)})}{k_n h} \int_{-h}^{-d} V(z') \psi_n(z')\,\mathrm{d} z' \end{equation}
in 
$x > N \delta$. It also follows that
\begin{align} R_N = 1 + \frac{\mathrm{i}}{kh} \int_{-h}^{-d} U(z) \psi_0(z)\,\mathrm{d} z \quad \mbox{and} \quad T_N =-\frac{\mathrm{i} \exp({-\mathrm{i} kN\delta})}{kh} \int_{-h}^{-d} V(z) \psi_0(z)\,\mathrm{d} z. \end{align}
The matching across 
$x=0$ and 
$x = N \delta$ requires some work since the representation of the solution in 
$x < 0$ and 
$x > N \delta$ in terms of eigenfunctions in 
$z$ is fundamentally different from the representation of the solution in 
$0 < x < N \delta$ which is based on Bloch–Floquet eigenmodes. We start by using (4.3) and (4.4) with 
$n=0$ to give
\begin{align} c_{1}^{(\,j)} E^{({+}j)} &= \int_{-h}^{-d} \left[ \phi_{1}(0,z) \frac{\partial \phi^{(-j)}}{\partial x}(0,z) - \phi^{(-j)}(0,z) \frac{\partial \phi_{1}}{\partial x}(0,z)\right] \mathrm{d} z \nonumber\\ &= \int_{-h}^{-d} \left[ \phi_{0}(0,z) \frac{\partial \phi^{(-j)}}{\partial x}(0,z)- \phi^{(-j)}(0,z) U(z)\right] \mathrm{d} z \end{align}and
\begin{align} d_{1}^{(\,j)} E^{(-j)} &= \int_{-h}^{-d} \left[ \phi_{1}(0,z) \frac{\partial \phi^{({+}j)}}{\partial x}(0,z) - \phi^{({+}j)}(0,z) \frac{\partial \phi_{1}}{\partial x}(0,z)\right] \mathrm{d} z \nonumber\\ &= \int_{-h}^{-d} \left[ \phi_{0}(0,z) \frac{\partial \phi^{({+}j)}}{\partial x}(0,z) - \phi^{({+}j)}(0,z) U(z)\right] \mathrm{d} z, \end{align}where the matching conditions (4.14a,b) and (4.15a,b) have been applied. On account of the relations (2.7) and (2.8) and using (4.16), we can rewrite (4.19) and (4.20) as
\begin{align} c_{1}^{(\,j)} \exp({-\mathrm{i} \beta^{(\,j)} \delta}) E^{({+}j)} &=-2H_{j0}-\sum_{n=0}^\infty\frac{ H_{jn}}{k_n h} \int_{-h}^{-d} U(z') \psi_n(z')\,\mathrm{d} z' \nonumber\\ &\quad -\int_{-h}^{-d} \phi^{({+}j)}(0,z) U(z)\,\mathrm{d} z \end{align}and
\begin{equation} d_{1}^{(\,j)} E^{(-j)}=2H_{j0} +\sum_{n=0}^\infty \frac{ H_{jn}}{k_n h}\int_{-h}^{-d} U(z') \psi_n(z')\,\mathrm{d} z'-\int_{-h}^{-d} \phi^{({+}j)}(0,z) U(z)\,\mathrm{d} z, \end{equation}where
\begin{equation} H_{jn}=\int_{-h}^{-d}\frac{\partial \phi^{({+}j)}}{\partial x}(0,z) \psi_n(z)\,\mathrm{d} z. \end{equation}
Using (4.5) with 
$n=N$ and following the same procedure give
\begin{align} c_1^{(\,j)} \exp({\mathrm{i} \beta^{(\,j)} (N-1) \delta}) E^{({+}j)} =\sum_{n=0}^\infty \frac{H_{jn}}{k_n h} \int_{-h}^{-d} V(z') \psi_n(z')\,\mathrm{d} z'- \int_{-h}^{-d} \phi^{({+}j)}(0,z) V(z)\,\mathrm{d} z \end{align}and
\begin{align} d_1^{(\,j)} \exp({-\mathrm{i} \beta^{(\,j)} N \delta}) E^{(-j)} =-\sum_{n=0}^\infty\frac{H_{jn}}{k_n h} \int_{-h}^{-d} V(z') \psi_n(z')\,\mathrm{d} z' -\int_{-h}^{-d} \phi^{({+}j)}(0,z) V(z)\,\mathrm{d} z. \end{align}
The algebraic manipulations above allow us to express the coefficients 
$c_1^{(\,j)}$ and 
$d_1^{(\,j)}$, and hence the solutions 
$\phi _1$ and 
$\phi _N$ in 
$0 < x < \delta$ and 
$(N-1) \delta < x< N \delta$ (respectively) are in terms of the unknown functions 
$U(z)$ and 
$V(z)$. This replicates what we had already achieved in (4.16) and (4.17) for 
$\phi _0$ and 
$\phi _{N+1}$ in 
$x < 0$ and 
$x > N \delta$.
Eliminating 
$c_1^{(\,j)}$ in (4.21) and (4.24) eventually results in
\begin{equation} \int_{-h}^{-d} U(z) {\mathcal{K}}^{(1)}_j(z)\,\mathrm{d} z+ \int_{-h}^{-d} V(z) {\mathcal{K}}^{(2)}_j(z)\,\mathrm{d} z =-2\exp({\mathrm{i} \beta^{(\,j)} N \delta})H_{j0}, \end{equation}where
\begin{equation} {\mathcal{K}}^{(1)}_j(z) = \exp({\mathrm{i} \beta^{(\,j)} N \delta}) \left[\sum_{n=0}^\infty\frac{ H_{jn} }{k_n h} \psi_n(z) + \phi^{({+}j)}(0,z) \right] \end{equation}and
\begin{equation} {\mathcal{K}}^{(2)}_j(z) =\sum_{n=0}^\infty\frac{H_{jn} }{k_n h} \psi_n(z) -\phi^{({+}j)}(0,z). \end{equation}
Also eliminating 
$d_1^{(\,j)}$ in (4.22) and (4.25) gives
\begin{equation} \int_{-h}^{-d} U(z) {\mathcal{K}}^{(2)}_j(z)\,\mathrm{d} z +\int_{-h}^{-d} V(z) {\mathcal{K}}^{(1)}_j(z)\,\mathrm{d} z =-2 H_{j0}. \end{equation}We note that if we write
\begin{equation} U^s(z) = U(z) + V(z) \quad \mbox{and} \quad U^a(z) = U(z) - V(z), \end{equation}then (4.26) and (4.29) decouple into a pair of scalar integral equations:
\begin{equation} \int_{-h}^{-d} U^{s,a}(z) {\mathcal{K}}^{s,a}_j(z)\,\mathrm{d} z =-2 H_{j0}, \end{equation}where
\begin{align} {\mathcal{K}}^{s}_j(z) &= \frac{1}{\exp({\mathrm{i} \beta^{(\,j)} N \delta}) +1} \left[{\mathcal{K}}^{(1)}_j(z) + {\mathcal{K}}^{(2)}_j(z) \right] \nonumber\\ &= \sum_{n=0}^\infty\frac{ H_{jn} }{k_n h} \psi_n(z) +\mathrm{i} \tan \left(\beta^{(\,j)} N \delta/2 \right) \phi^{({+}j)}(0,z) \end{align}and
\begin{align} {\mathcal{K}}^{a}_j(z) &= \frac{1}{\exp({\mathrm{i} \beta^{(\,j)} N \delta}) -1} \left[{\mathcal{K}}^{(1)}_j(z) - {\mathcal{K}}^{(2)}_j(z) \right] \nonumber\\ &= \sum_{n=0}^\infty\frac{ H_{jn} }{k_n h} \psi_n(z) -\mathrm{i} \cot \left(\beta^{(\,j)} N \delta/2 \right) \phi^{({+}j)}(0,z). \end{align}
The use of superscripts 
$s$ and 
$a$ indicates that the two integral equations can be thought of as representing the components of the scattering by the barrier array which are symmetric and antisymmetric about the mid-plane, 
$x=N \delta /2$, of geometric symmetry.
4.1. Numerical approximation
The pair of integral equations (4.31) is approximated numerically using the identical method in § 2.3. Thus the two functions 
$U^s(z)$ and 
$U^a(z)$ have the same form as (2.41):
\begin{equation} U^{s,a}(z) \approx \sum_{n=0}^M \alpha_n^{s,a} u_n(z), \end{equation}
where 
$u_n(z)$ are defined in (2.42).
After substituting (4.34) into (4.31) we can obtain the following systems of equations:
\begin{equation} \sum_{n=0}^M \alpha_n^{s,a} K_{mn}^{s,a} =-2 H_{m0}, \end{equation}
for 
$m=0,1,\ldots,M$, where
\begin{equation} K_{mn}^s = \sum_{r=0}^\infty \frac{H_{mr}F_{nr}}{k_r h} +\mathrm{i} \tan \left(\beta^{(m)} N \delta/2 \right) G_{mn} \end{equation}and
\begin{equation} K_{mn}^a = \sum_{r=0}^\infty \frac{H_{mr}F_{nr}}{k_r h} -\mathrm{i} \cot \left(\beta^{(m)} N \delta/2 \right) G_{mn}, \end{equation}
in which 
$F_{mr}$ have already been defined in (2.45) and
\begin{equation} G_{mn} = \int_{-h}^{-d} u_n(z)\phi^{({+}m)}(0,z)\,\mathrm{d} z. \end{equation}
The system of equations (4.35) has been truncated with the parameter 
$M$ which need not be the same as 
$M_{k}$ in § 2.3. From (4.34) we can see that 
$M$ denotes the number of the vertical eigenfunctions used to approximate the horizontal velocity on the interface; on the other hand, from (4.36), (4.37) and (4.38) we can see that 
$M$ also represents the number of the evanescent modes applied to simulate the oscillation in the barrier array.
In the Bloch–Floquet problem, two solutions of 
$\phi ^{(+m)}(0,z)$ in different forms are presented in § 2 such that the eigenfunctions can be approximated by either (2.17) or (2.25) resulting in that 
$G_{mn}$ can be written as
\begin{equation} G_{mn}=\sum_{l=-\infty}^{\infty} \frac{(-1)^n b_{1,l}^{(m)}}{\cosh \beta_{l}^{(m)}(h-d)}I_{2m} \left[\beta_{l}^{(m)}(h-d) \right] \end{equation}or
\begin{equation} G_{mn}=\sum_{l=0}^{\infty} \left(a_{2,l}^{(m)}+b_{2,l}^{(m)} \right) F_{nl}. \end{equation}
The series in 
$G_{mn}$ decays like 
$O(1/l^2)$ as 
$l$ tends to infinity, which can be treated with the same procedure shown in Appendix A. However, for 
$H_{mr}$ in (4.23) we find that either solution will produce a slowly convergent series decreasing like 
$O(1/l^{3/2})$ if the expressions of 
$\phi ^{(+m)}(0,z)$ are applied directly. In order to calculate 
$H_{mr}$ efficiently, we use the approximation in (2.30) which is based on the second form in § 2. Then, 
$H_{mr}$ can be written as
\begin{equation} H_{mr}=\sum_{k=0}^{M_2} \alpha_{2,k}^{(m)} F_{kr} \end{equation}
after using (2.41), where 
$\alpha _{2,k}^{(m)}$ are eigenvectors corresponding to the eigenvalue 
$\beta ^{(m)}$. Now a slowly convergent infinite series is replaced with a truncated series such that the series in (4.36) and (4.37) decaying like 
$O(1/r^2)$ also can be computed efficiently if the method in Appendix A is used.
Once 
$\alpha _n^{s,a}$ have been determined from (4.35), we can recover the reflection and transmission coefficients from the use of (4.30a,b) in (4.18a,b) with (4.34) and (2.45) to give
\begin{align} R_N = 1 + \frac{\mathrm{i}}{2kh} \sum_{n=0}^M (\alpha_n^s+\alpha_n^a) F_{m0}\quad \mbox{and} \quad T_N =-\frac{\mathrm{i} \exp({-\mathrm{i} kN\delta})}{2kh} \sum_{n=0}^M (\alpha_n^s-\alpha_n^a) F_{m0}. \end{align} Combined with the description in § 2.3, it can be seen that the two different forms of eigenfunction 
$\phi ^{(\pm m)}(x,z)$ have their own advantages. The first form is better for approximating solutions to the Bloch–Floquet problem, especially for closely spaced barriers and the explicit limit of vanishing spacing can be taken. On the other hand, with the help of the second form, the slowly convergent series appearing in the scattering problem for a finite periodic array can be treated efficiently. In addition, from (4.36) and (4.37) we also can see that the present method has the same advantages as the recursive transfer matrix method (e.g. Porter & Porter Reference Porter and Porter2003) in that the dimension of the equation system is independent of the size of the array.
Furthermore, it should be noted that the reflection and transmission coefficients (4.42a,b) cannot be applied for the case of the critical frequency (i.e. 
$\beta ^{(0)}={\rm \pi} / \delta$) since the eigenfunctions 
$\phi ^{(\pm 0)}$ no longer satisfy the orthogonality relation (2.11) which has been widely used in the above derivation.
5. Scattering using the continuum model
Consider that the region 
$0 < x < N \delta$ is governed by the continuum model described in § 3, so that the potential in this region may be written as
\begin{equation} \phi_h(x,z) = \sum_{n=0}^\infty \left( c_n \exp({\mathrm{i} \mu_n x}) + d_n \exp({\mathrm{i} \mu_n (N \delta- x)}) \right) Z_n(z), \end{equation}where
\begin{equation} Z_n(z) = \varepsilon_n^{-1/2} \left\{\begin{array}{@{}ll} (1+Kz)/(1-Kd), & -d < z < 0, \\ \cosh \mu_n(h+z)/\cosh \mu_n(h-d), & -h < z <-d, \end{array}\right. \end{equation}
and 
$\mu _n$ are the roots of (3.17). Also, in this section a new labelling rule is used for 
$\mu _n$ that the real value takes precedence over increasing pure imaginary values and if there does not exist the real value the smallest value on the imaginary axis takes 
$\mu ^{(0)}$. In (5.2), 
$\varepsilon _n$ are normalisation factors given by
\begin{equation} \varepsilon_n = \frac{1}{2\cosh^2 \mu_n(h-d)} \left[ 1 + \frac{\sinh 2 \mu_n (h-d)}{2 \mu_n (h-d)} \right], \end{equation}so that
\begin{equation} \frac{1}{h-d} \int_{-h}^{-d} Z_n(z) Z_m(z)\,\mathrm{d} z = \delta_{mn}. \end{equation}
This orthogonality relation for the functions 
$Z_n(z)$ follows since 
$\mu _n$ are distinct and
\begin{align} (\mu_n^2 - \mu_m^2) \int_{-h}^{-d} Z_n(z) Z_m(z) \mathrm{d} z &= \int_{-h}^{-d} \left[Z_n''(z) Z_m(z) - Z_n(z) Z_m''(z) \right] \mathrm{d} z \nonumber\\ &= \left[ Z_n'(z) Z_m(z) - Z_n(z) Z_m'(z) \right]_{-h}^{-d} = 0, \end{align}where the dispersion relation (3.17) has been used.
The matching conditions shown in (4.13a,b), (4.14a,b) and (4.15a,b) still hold but 
$\phi _h$ will replace 
$\phi _1$ and 
$\phi _N$. Continuity of the horizontal velocity at 
$x=0$ and 
$x=N \delta$ results in
\begin{equation} c_n - d_n \exp({\mathrm{i} \mu_n N \delta}) = \frac{1}{\mathrm{i} \mu_n (h-d)} \int_{-h}^{-d} U(z) Z_n(z)\,\mathrm{d} z \end{equation}and
\begin{equation} c_n\exp({\mathrm{i} \mu_n N \delta}) - d_n = \frac{1}{\mathrm{i} \mu_n (h-d)} \int_{-h}^{-d} V(z) Z_n(z)\,\mathrm{d} z . \end{equation}This gives
\begin{equation} c_n = \frac{1}{2 \mu_n(h-d) \sin \mu_n N \delta} \int_{-h}^{-d} \left[ U(z) \exp({-\mathrm{i} \mu_n N \delta}) - V(z) \right] Z_n(z)\,\mathrm{d} z \end{equation}and
\begin{equation} d_n = \frac{1}{2 \mu_n(h-d) \sin \mu_n N \delta} \int_{-h}^{-d} \left[ U(z) - V(z)\exp({-\mathrm{i} \mu_n N \delta}) \right] Z_n(z)\,\mathrm{d} z. \end{equation}
Matching (5.1) to (4.16) across 
$x=0$ gives
\begin{equation} \int_{-h}^{-d} U(z') {\mathcal{L}}^{(1)}(z,z')\,\mathrm{d} z' + \int_{-h}^{-d} V(z') {\mathcal{L}}^{(2)}(z,z')\,\mathrm{d} z' =-2 \psi_0(z) \end{equation}
and to (4.17) across 
$x = N \delta$ gives
\begin{equation} \int_{-h}^{-d} U(z'){\mathcal{L}}^{(2)}(z,z')\,\mathrm{d} z' +\int_{-h}^{-d} V(z'){\mathcal{L}}^{(1)}(z,z')\,\mathrm{d} z' = 0, \end{equation}where
\begin{equation} {\mathcal{L}}^{(1)}(z,z')=\sum_{n=0}^\infty \left[\frac{\psi_n(z) \psi_n(z')}{k_n h} - \frac{Z_n(z) Z_n(z')}{\mu_n (h-d) \tan(\mu_n N \delta)} \right] \end{equation}and
\begin{equation} {\mathcal{L}}^{(2)}(z,z')=\sum_{n=0}^\infty \frac{Z_n(z) Z_n(z')}{\mu_n (h-d) \sin (\mu_n N \delta)}. \end{equation}For (5.10) and (5.11), we also can decouple the pair of integral equations into their symmetric and antisymmetric components like (4.30a,b) which satisfy
\begin{equation} \int_{-h}^{-d} U^{s,a}(z'){\mathcal{L}}^{(s,a)}(z,z')\,\mathrm{d} z' =-2 \psi_0(z) , \end{equation}where
\begin{align} {\mathcal{L}}^{(s)}(z,z') &= {\mathcal{L}}^{(1)}(z,z')+{\mathcal{L}}^{(2)}(z,z') \nonumber\\ &= \sum_{n=0}^\infty \left[\frac{\psi_n(z) \psi_n(z')}{k_n h} + \frac{\tan (\mu_n N \delta/2)}{\mu_n (h-d) } Z_n(z) Z_n(z')\right] \end{align}and
\begin{align} {\mathcal{L}}^{(a)}(z,z') &= {\mathcal{L}}^{(1)}(z,z')-{\mathcal{L}}^{(2)}(z,z') \nonumber\\ &= \sum_{n=0}^\infty \left[\frac{\psi_n(z) \psi_n(z')}{k_n h} - \frac{\cot (\mu_n N \delta/2)}{\mu_n (h-d) } Z_n(z) Z_n(z')\right]. \end{align}
These are the equations that would be derived had the original problem been decomposed into the sum of problems symmetric and antisymmetric about the mid-plane 
$x=N\delta /2$.
5.1. Numerical approximation
The approximation (4.34) will be used again. We substitute (4.34) into (5.12), multiply through by 
$u_m(z)$ and integrate over 
$-h < z < -d$, a process which characterises the Galerkin method and results in the following systems of equations:
\begin{equation} \sum_{n=0}^M \alpha_n^{s,a} L_{mn}^{s,a} =-2 F_{m0}, \end{equation}
for 
$m=0,1,\ldots,M$, where
\begin{equation} L_{mn}^s = \sum_{r=0}^\infty \frac{F_{mr} F_{nr}}{k_r h} + \sum_{r=0}^\infty \frac{\tan (\mu_r N \delta/2)}{\mu_r(h-d) } P_{mr} P_{nr} \end{equation}and
\begin{equation} L_{mn}^a = \sum_{r=0}^\infty \frac{F_{mr} F_{nr}}{k_r h} - \sum_{r=0}^\infty \frac{\cot (\mu_r N \delta/2)}{\mu_r(h-d) }P_{mr} P_{nr}, \end{equation}in which
\begin{equation} P_{mr} = \int_{-h}^{-d} u_m(z) Z_r(z)\,\mathrm{d} z = \varepsilon_m^{-1/2} (-1)^m I_{2m} [\mu_r(h-d)]. \end{equation}
All of the series in (5.18) and (5.19) have the order of 
$O(1/r^2)$ when 
$r \to \infty$, so the treatment shown in Appendix A can be applied to accelerate the series convergence. After the systems of equations are solved numerically, the reflection and transmission coefficients can be determined also by (4.42a,b).
It can be seen that the derivation for the reflection and transmission coefficients between the discrete model and the continuum model is different. For the discrete model, the number of the truncated evanescent mode 
$M$ is equal to the dimension of the equation system (see (4.35)). As shown later, 
$M=12$ is usually sufficient to obtain convergent results. For the continuum model, the evanescent mode is included in the series shown in (5.18) and (5.19). We generally need to consider roughly 1000 terms to guarantee the accuracy of the series calculation. Actually, for the continuum model, we also can develop a system of equations similar to (4.35). However, we find that this would present a troublesome series when close to the critical frequency.
5.2. Results
We first examine the convergence of the scheme for the discrete model and the continuum model. In both settings, the barriers are submerged to a depth 
$d/h=0.2$ and the total distribution length of the barriers is 
$N\delta =h$. Table 3 shows how the modulus of the reflection coefficient, 
$|R_N|$, converges with the truncation parameter, 
$M$. At low frequencies, results can be seen to converge quickly requiring only a small system of equations, but when the frequency approaches the critical frequency results tend to converge slowly with 
$M$ since the amplitude of the fluid oscillation in the barrier array structure becomes increasingly severe. When the frequency is in the stop band, the wave motion decays through the array and there is practically no transmission. As mentioned at the end of § 4.1, the second method in § 2, which tends to converge fastest when 
$\delta /d$ is larger, is used in this scattering problem for determining the slowly convergent series. Thus, for the discrete model, the convergent results for distribution with large spacing can be obtained by a small truncation parameter used. In general, 
$M=12$ is sufficient to produce results with the accuracy of roughly four significant figures although computations are more demanding when close to the critical frequency.
The convergence of the modulus of reflection coefficient 
$|R_N|$ computed using the discrete model and the continuum model against the truncation parameter, 
$M$, in the case of 
$d/h=0.2$ and 
$N\delta =h$.
Next, some cases have been chosen to allow the comparison between the discrete model based on an expansion in terms of Bloch–Floquet eigenfunctions and existing results. A pair of barriers 
$(N=1)$ with the submergence 
$d/h=0.2$ is first examined for which Porter & Evans (Reference Porter and Evans1995) previously provided accurate computations using the Galerkin approximation method. As shown in figure 5, the results of the discrete model compare favourably with these existing results, accurately replicating total reflection and transmission. Notice that the heavily suppressed transmission beyond 
$Kd \approx 1$ can now be understood as being associated with the stop band for the periodic barrier array despite there only being two barriers and one cell in the present example.
Comparison of the modulus of the reflection coefficient 
$|R_1|$ for the discrete model with existing results (Porter & Evans Reference Porter and Evans1995) for the case of 
$N=1$ and 
$d/h=0.2$: (a) 
$\delta /d=0.50$; (b) 
$\delta /d=0.05$.
When the number of barriers 
$N+1$ is large, direct solution methods such as those used by Porter & Evans (Reference Porter and Evans1995) are algebraically cumbersome and lead to 
$N+1$ coupled equations in terms of 
$N+1$ unknown functions eventually implying that numerical computations are of 
$O(N^3)$. To mitigate against this, previous authors (e.g. Porter & Porter Reference Porter and Porter2003) have used transfer matrices in which the scattering by 
$N+1$ elements of the array is accounted for by the multiplication of 
$N+1$ matrices whose size depends on the number of evanescent wave interactions retained in the exchange of information between adjacent elements in the array. Superficially the computational effort is of 
$O(N)$. In figure 6, we fix the number of barriers 
$N=10$ and the submergence 
$d/h=0.2$ and compare the modulus of the reflection coefficient 
$|R_{10}|$ computed using the present Bloch–Floquet discrete model (which does not scale with 
$N$) with those computed using transfer matrices. The results agree well, apart from when very close to the critical frequency where resonance occurs and when 
$N$ is large. On account of the high-frequency oscillations in 
$|R_{N}|$ close to the critical frequency, even small errors in either the transfer matrix method or the present approach can lead to large changes in 
$|R_{N}|$ and it is not easy to determine which is more accurate. In particular, as the spacing decreases an increasing number of evanescent modes is required to maintain accurate computations resulting in larger transfer matrices and, in turn, this leads to numerical instability caused by rounding errors even though a treatment for avoiding these rounding errors devised by Porter & Porter (Reference Porter and Porter2003) has been applied. Thus, for the case of 
$\delta /d=0.05$, the calculation by the method of the transfer matrix fails and figure 6(c) only includes the results from the present discrete model. During the review of this paper, one of the reviewers pointed out that the method devised by Ko & Sambles (Reference Ko and Sambles1988) may be used to overcome the numerical instability issues caused by using transfer matrices.
Comparison of the modulus of the reflection coefficient 
$|R_{10}|$ for the discrete model with the scattering matrix for the case of 
$N=10$ and 
$d/h=0.2$: (a) 
$\delta /d=0.50$; (b) 
$\delta /d=0.25$; (c) 
$\delta /d=0.05$.
As mentioned in the Introduction, one aim of the present study is to assess the validity of the homogenisation method for wave interaction with plate array structures. We first investigate the validity of the homogenisation method by varying the number of barriers (or the total length of the barrier array). In figure 7, the case of 
$\delta /d=0.05$ and 
$d/h=0.2$ is investigated with 
$N=1$, 
$N=5$ and 
$N=10$ and results from the exact discrete model are plotted against the results from the homogenisation approximation. We recall that homogenisation is not expected to work for 
$Kd$ sufficiently close to a value of 
$K_c d = 1$ corresponding to the critical frequency and curves will oscillate infinitely quickly as 
$Kd = 1$ is approached. On the other hand, the curves of 
$|R_N|$ computed under the discrete model oscillate and the total transmission will happen 
$N$ times before the critical frequency is reached. Thus, we see in figure 7 overall good agreement between the exact and approximate models apart from close to 
$Kd = 1$.
The variation of the modulus of reflection coefficient 
$|R_N|$ for the discrete model and the continuum model against the non-dimensional frequency 
$Kd$ for the case of 
$\delta /d=0.05$ and 
$d/h=0.2$: (a) 
$N=1$; (b) 
$N=5$; (c) 
$N=10$.
In figure 8 we present the modulus of the reflection coefficient 
$|R_N|$ for the discrete model and the continuum model. The total length of the barrier array is fixed (i.e. 
$N\delta =h$), but the spacing in different arrays varies. This plot allows us to see how scattering computed from the discrete model converges to the results predicted from homogenisation. Again the submergence is 
$d/h=0.2$ (the depth of the fluid is relatively unimportant for the effects we are observing). The spacing 
$\delta /d = 0.5$ corresponds to 
$N=10$ whilst 
$\delta /d = 0.05$ corresponds to 
$N=50$. The two curves corresponding to these two cases hit the horizontal axis (i.e. 
$|R_N| = 0$) 10 and 50 times respectively (although this cannot be captured by the resolution in the plots). We can see the agreement is good for low frequencies and gets better as 
$\delta /d$ decreases, although the rapid oscillations in 
$|R_N|$, which occur as the critical frequency is approached, mean that the models diverge for 
$Kd$ sufficiently close to one. The continuum model serves as a good approximation for the spacing 
$\delta /d<0.05$ when 
$Kd \lesssim 0.7$ or 
$\delta /d<0.5$ when 
$Kd \lesssim 0.4$.
The variation of the modulus of reflection coefficient 
$|R_N|$ for the discrete model and the continuum model against the non-dimensional wavenumber 
$Kd$ for the case of 
$N\delta =h$ and 
$d/h=0.2$.
Finally, in figure 9 we plot the velocity potential field for two cases of scattering of incident waves by a barrier array computed using the continuum model for 
$d/h=0.2$ and 
$N \delta = h$ at 
$Kd = 0.8$ and 
$Kd = 1.2$, that is, above and below the critical frequency respectively. In the first case, we can see wave propagation through the barrier array leading to transmission beyond the array. In the second case, we see rapid decay of the wave field through the array and near perfect reflection of incident waves. We have been unable to show a field plot from the direct numerical approach since the solution to the scattering problem requires that we used the second method presented in § 2 to determine eigenfunctions. As explained at the end of § 2.3, this method is poor at producing convergent representations for the field.
The fields of the real part of the velocity potential for the continuum model with 
$d/h=0.2$ and 
$N\delta =h$: (a) 
$Kd=0.8$; (b) 
$Kd=1.2$.
6. Conclusions
The main focus of this paper has been on describing and comparing two approaches to solving the problem of two-dimensional wave propagation through periodic arrays of surface-piercing barriers with a particular focus on the small spacing between adjacent barriers. A continuum model is described which is derived formally for small barrier spacing using homogenisation methods. This model is shown to be valid away from resonance occurring at 
$\omega = \sqrt {g/d}$; the propagating wavelength is predicted to become vanishingly small as resonance is approached, signalling a breakdown in the multiple-scale assumption underpinning homogenisation. This conclusion sheds light on previously unexplained ill-posed behaviour associated with the use of a continuum description for wave interaction with resonance in plate arrays in problems encountered by, for example, Jan & Porter (Reference Jan and Porter2018) and Zheng et al. (Reference Zheng, Porter and Greaves2020). A more complicated approach is based on an exact description of the barrier array for non-zero spacing, 
$\delta$. First, by considering propagation in infinite periodic barrier arrays we have been able to show that there is a critical frequency 
$\omega _c$ which lies below 
$\sqrt {g/d}$ and acts to divide wave propagation from wave decay. At the critical frequency, standing waves exist between the barriers in the array and carry no energy. As 
$\delta \to 0$, this model tends to the continuum model provided frequencies are sufficiently far away from resonance. The approach shows that, for a non-zero spacing, there is a well-behaved transition from passing to stopping associated with wavelengths on the scale of the separation 
$\delta$.
The exact description in terms of non-zero 
$\delta$ for a finite number of 
$N+1$ barriers is also considered. Use is made of an orthogonality relation which applies to eigenmodes derived from the infinite periodic array to express the solution through the entire finite barrier array region in terms of the solution in just one period. The results show that the limit of the discrete problem is the continuum problem provided resonance is avoided.
Although this problem has been developed for a simple geometry where extensive use of separation of variables has been made to develop a semi-analytical approach to the problem in terms of solutions to a pair of scalar integral equations, it is clear that there will be other problems in water waves, linearised acoustics, etc. involving scattering by finite periodic arrays which can be analysed by the same method. In particular, the orthogonality condition satisfied by the eigenfunctions for the periodic Bloch problem is key to connecting solutions from one edge of a periodic array to the other. Its construction is in § 2.1 of the paper, which is not dependent on the geometry of the structure and can be applied to fields governed by the Helmholtz equation, for example.
Following this work, it would be interesting to consider wave propagation through an array of barriers in which the barrier length is a slowly varying function of space. For example, the methods described by Porter (Reference Porter2020) could be used to transform the scattering process into an ordinary differential equation in the continuum limit 
$\delta \to 0$ as a means of qualitatively understanding the onset of resonance/rainbow reflection in graded arrays such as those described by Wilks et al. (Reference Wilks, Montiel and Wakes2022).
Funding
This work was supported by the EPSRC, grant number EP/V04740X/1.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Numerical treatment for the slowly convergent series in $L_{mn}$ and $K_{mn}$
In § 2.3, the matrices formed by 
$L_{mn}$ and 
$K_{mn}$ are presented for determining the Bloch–Floquet wavenumber 
$\beta$. It can be found that all of the series in 
$L_{mn}$ and 
$K_{mn}$ are convergent as of 
$O(1/r^2)$. In order to speed up the series convergence, the last term in (2.40) (taken as an example) can be written as
\begin{align} \sum_{r=-\infty}^\infty l_{mn,r} &= l_{mn,0}+\frac{\left[1+(-1)^{m+n} \right]^2}{24{\rm \pi}} \nonumber\\ &\quad + \sum_{r=1}^\infty \left[ l_{mn,r}-\frac{\cos\left[{\tfrac12} (m-n){\rm \pi}\right]+ \sin\left[\beta \delta- {\tfrac12} (m+n){\rm \pi}\right]}{2r^2{\rm \pi}^3} \right] \nonumber\\ &\quad +\sum_{r=-\infty}^{-1} \left[ l_{mn,r}-(-1)^{m+n}\frac{ \cos\left[{\tfrac12} (m-n){\rm \pi}\right]- \sin\left[\beta \delta+ {\tfrac12} (m+n){\rm \pi}\right]}{2r^2{\rm \pi}^3} \right], \end{align}with
\begin{equation} l_{mn,r}=\frac{{\rm J}_m\left(\tfrac12 \beta_r \delta\right) {\rm J}_n\left(\tfrac12 \beta_r \delta\right)}{\beta_r \delta \tanh \beta_r (h-d)}, \end{equation}
where the asymptotic form of Bessel function has been used, i.e. 
$\textrm {J}_n(z)\thicksim \sqrt {2/({\rm \pi} z)} \cos (z-n{\rm \pi} /2-{\rm \pi} /4)$ when 
$|z| \to \infty$ and 
$|\arg z|<{\rm \pi}$ (see Abramowitz & Stegun Reference Abramowitz and Stegun1972). The infinite series in (A1) now decay like 
$O(1/r^4)$ such that results can be efficiently computed to high accuracy (we aim for an error of less than 
$10^{-8}$). For other slowly convergent series, we apply the same treatment.
Open access
