A wide-spacing approximation model for the reflection and transmission of water waves over an array of vertical obstacles

Abstract With a view to numerical modelling and optimisation of wave energy farms, a simple recursive formulation is employed to solve for the reflection and transmission of plane water waves by a number of rows of vertical obstacles, under the wide-spacing approximation. The proposed recursive formulation relies on the ‘concatenation’ of any two sets of obstacles, for which the reflection–transmission problem is already resolved. Furthermore, the obstacles are allowed to move in any combination of pitch and surge. The proposed recursive model is validated by means of physical experiments in a small-scale wave flume, whereby waves are reflected and transmitted by one, two and three rows of vertical, flexible blades, taking into account dissipation within the fluid along the wave propagation direction. For the special case of identical, regularly spaced rows, under the adopted formalism, distinct theoretical behaviours are highlighted, depending on whether or not individual obstacles absorb (or dissipate) energy as they interact with incoming waves. In a ‘non-dissipative’ case, the well known fact that discrete values of the row-to-row distance $L$ completely cancel reflection is retrieved, as well as the existence of ‘band-gap’ intervals, i.e. intervals for $L$ where reflection is high, with maximum reflection occurring away from the Bragg condition. In contrast, when the obstacles dissipate or absorb energy as they interact with the fluid, reflection is always non-zero, and, as the number of rows tends to infinity, forms marked Bragg peaks, reaching unity when $L$ is a multiple of half a wavelength.


Introduction
Ocean waves transport massive amounts of energy. In terms of power per metre of wavefront, the monthly median power from wind waves above 30 • N for instance, ranges from 17 to 130 kW m −1 (Arinaga & Cheung 2012). Wave energy has the potential to be transformed into useful energy (Callaway 2007). In the case of Europe, it could be a significant contributor to the electricity supply, with an estimated 300-400 GW potential along European Atlantic coastlines alone (Babarit 2017). In addition to being a source of clean, renewable power, large arrays of wave energy absorbing structures may also serve the objective of mitigating coastal erosion along the shoreline (Nové-Josserand et al. 2018), mimicking the wave reduction effect of natural coastal defences such as mangroves, salt marshes or seagrass and kelp beds (Abanades, Greaves & Iglesias 2015;Narayan et al. 2016).
Among wave energy converter concepts, oscillating wave surge converters (OWSCs), which primarily exploit the horizontal fluid motion, are identified as one of the most promising and mature technological options (Babarit 2017). Traditional OWSC designs consist of a rigid flap, oscillating around a rotation point, which can be fixed to the sea bed or floating (Sarkar, Renzi & Dias 2014). Recent work (Nové-Josserand 2018;Nové-Josserand et al. 2018) has investigated an artificial canopy of bio-inspired, flexible OWSCs, with the two objectives of coastline protection and energy harvesting. These arrays of flexible structures may indeed present benefits in terms of survivability and absorption capabilities, in the same way that aquatic vegetation can withstand and dissipate surface wave energy (Koehl & Wainwright 1977;Denny & Gaylord 2002). From the theoretical point of view, on the one hand, each flexible element of the artificial canopy acts as an oscillator, with its intrinsic natural frequency or damping coefficient. On the other hand, the array as a whole can be viewed as a metamaterial, with properties related to its internal structure, such as the existence of crystallographic effects analogue to those observed in solid-state physics or acoustics, e.g. Bragg resonances (Garnaud & Mei 2009;Arnaud et al. 2017;Rey et al. 2018). Several studies have explored these ideas in different systems designed in view of engineering water wave propagation, ranging from the refraction phenomenon of water waves propagating through an array of bottom-mounted structures (Hu & Chan 2005;Arnaud et al. 2017), to the tuning of the sea bed topography (Davies & Heathershaw 1984;Berraquero et al. 2013), or the deployment of floating membranes with a crystalline array of defects that confer unique propagation features to the hydroelastic waves that ensue (Domino, Fermigier & Eddi 2020). Rey et al. (2018) study the propagation of water waves through dense and sparse arrays of vertical, bottom-mounted cylinders, including the reflection of oblique flume modes in the analysis. For dense arrays (where cylinder spacing in the wave direction is much smaller than the wavelength), the cylinder array is better treated as a porous medium, while Bragg reflection patterns become governed by the array length, as opposed to the cylinder spacing (Arnaud et al. 2017;Rey et al. 2018).
In order to assess the potential, and improve the design of such flexible OWSC arrays, it is necessary to articulate individual OWSC parameters, together with the array layout, within a reasonably detailed dynamical model. Ideally, an appropriate dynamical model should be able to predict global quantities, such as the reflected, transmitted, absorbed and dissipated energy, as well as more detailed information on the dynamics of individual flexible OWSCs, such as wave loads, motion and energy absorption. The present work paves the way towards those objectives.
In the work of Nové-Josserand, Godoy-Diana & Thiria (2019), assuming that flexible OWSCs are arranged in several rows parallel to the wavefront, a simple interaction model is proposed, which allows for assessing the effects of various parameters upon the global reflection and transmission coefficients of the array (see also Renzi & Dias (2012)). Analytical formulations are derived, which relate the transmission and reflection coefficients of a single row, as well as the spacing between rows, to the transmission and reflection coefficients of the whole array. Those formulae, however, neglect successive reflections between non-neighbouring rows within the array. Furthermore, since reflection and transmission coefficients are represented as real numbers, the model does not account for the possible phase shift, induced by transmission and reflection.
In this work, building upon the ideas of Nové-Josserand et al. (2019), an improved interaction model, similar to the one in Evans (1990), is proposed, which does not neglect interactions between any pair of rows, and considers complex reflection and transmission coefficients. Using the wide-spacing approximation, the proposed model is based on an iterative approach, whereby the transmission-reflection problem is solved recursively as new rows are successively added to the array. The recursion proposed in the present study determines the reflection and transmission coefficients for the concatenation of any two arrays, of which the reflection and transmission coefficients are already known. Furthermore, the obstacles are allowed to move, and absorb energy, in any combination of surge and pitch. Finally, wave dissipation within the fluid along the propagation direction is taken into account within the model.
The interaction model is validated against the results of physical experiments carried out by Nové-Josserand et al. (2019), whereby several rows of vertical, flexible blades are subject to regular plane waves in a small-scale wave flume. Experimental results support the model accuracy.
Taking a step back from the specific problem of flexible OWSCs, the paper then focuses on the case of arrays made of identical, regularly spaced rows of obstacles. Exact formulae are retrieved for the array reflection and transmission coefficients, as functions of the properties of individual rows and array spacing, as well as asymptotic expressions when the number of rows tends to infinity. Although the derivation of those formulae does not differ, in principle, from existing work such as Linton (2011), in the present study the individual transmission and reflection coefficients are expressed in a formalism which emphasises the role of energy dissipation or absorption by the obstacles. As a result, the formulae obtained for the global array reflection and transmission highlight the existence of two significantly different behaviours, depending on whether or not the interaction between the wave and individual rows involves energy absorption (or dissipation).
The remainder of this paper is organised as follows. The problem and associated assumptions are defined in § 2. Section 3 presents a recursion relation to solve the reflection and transmission problem. In § 4, the predictions of the proposed model are compared with experimental results. In § 5, the special case of N identical rows is given special attention. Conclusions are provided in § 6.

Problem definition and assumptions
2.1. Overview of the problem Consider an array of thin, vertical structures arranged in parallel rows, such as that illustrated in figure 1, where the extent of each row is infinite along the y axis. The water depth is assumed constant. The undisturbed incident wave is a plane wave propagating in the x direction, orthogonal to the rows. The array is assumed periodic in the y direction, with periodicity smaller than the wavelength; in other words, λ > l + W. Therefore, transverse modes can be neglected (Dalrymple & Martin 1990). The incoming wave is described by means of the free-surface elevation (FSE), η 0 , given by the following equation: ( 2.1) where k = 2π/λ is the wavenumber, ω is the wave frequency andη 0 is the complex wave amplitude. As the incoming plane wave travels through a given row, the wave-row interaction results in a transmitted plane wave and a reflected plane wave. Transmitted and reflected wave amplitudes are deduced from the incoming plane wave through transmission and reflection coefficients, respectively. Such a simple representation omits the evanescent modes, which are essential to describe the flow in the close vicinity of each row of obstacles. This implies, in particular, that the rows are sufficiently far from each other (Evans 1990), hence the terminology 'wide-spacing approximation' (also known as 'plane wave approximation' in naval and offshore hydrodynamics). In theory, row-to-row distances greater than one wavelength are necessary for the wide-spacing approximation to be theoretically consistent (Martin 1984), although many studies (McIver & Bennett 1993;chapter 6 of Linton & McIver 2001;chapter 3 of Li 2006) find, through comparison with more accurate numerical techniques, that the wide-spacing approximation results remain accurate for distances significantly smaller than the wavelength. Also note that the wide spacing approximation does not assume that evanescent modes are non-existent: it merely assumes that they can be neglected when analysing row-to-row interaction. The wave field, as computed through the wide-spacing approximation, is only a partial description of the exact wave field, valid outside the close vicinity of the rows considered. Finally, the interaction theory presented in this work does not assume identical rows: the obstacle characteristics, as well as their spacing, may vary across rows, but not within a given row; furthermore, the distance between consecutive rows may also vary.

Diffraction problem for a single row of fixed obstacles
First consider a single row of fixed obstacles, which constitutes a diffraction problem. For example, concerning the obstacles illustrated in figure 1, this assumption implies that the blades are infinitely rigid. The fixed row, located at a position x 1 , is subject to an incoming wave described as in (2.1). The wave, transmitted by the obstacle, propagates in the positive x direction, and can be written as follows: (2. 2) The reflected wave takes the form In (2.2) and (2.3),t andr are complex-valued transmission and reflection coefficients, which apply a change in amplitude and a phase shift to the incoming wave as it reaches the obstacle position. The coefficientst andr depend on the obstacle geometry, and on the incoming wave frequency. Also note that, since the row is thin, it is considered to have a plane of symmetry transverse to the x axis; therefore, the reflection coefficient is identical for incoming waves propagating in the positive and in the negative x directions. It is possible to state a few general properties whichr andt must satisfy (Linton 2011). In linear wave theory, and considering no energy absorption or dissipation at the interface with the obstacles, preservation of energy implies that the row reflection and transmission coefficients satisfy the following equality: (2.4) Equation (2.4) can be verified, for example, in Ursell (1947), for the case where the rows of blades in figure 1 are replaced with infinitely long, submerged vertical barriers, extending along the whole y axis, or in Linton (2011) for the case of submerged, horizontal cylinder extending along the y direction.
In the more general case, where some of the incoming wave energy can be dissipated at the interface with the obstacle, energy conservation implies the following inequality (see e.g. Isaacson, Premasiri & Yang (1998) for a vertical slotted barrier): |t| 2 + |r| 2 ≤ 1. (2.5) Note that the energy-preservation property, whether it takes the form of (2.5) or (2.4), should also be satisfied by the array as a whole. Considering infinitely thin rows (or at least sufficiently thin with respect to the wavelength), the following relation must hold between the complex reflection and transmission coefficientsr andt:t +r = 1.
(2.6) Indeed, one may think of the row, excited by two incoming waves of identical amplitude, and propagating in opposite directions. The two incoming waves together form a standing wave pattern. At the antinodes of the standing waves (i.e. at the locations where the wave amplitude is the largest), the horizontal fluid velocity is zero along the whole water column. Therefore, by synchronising the two exciting waves in such a way that the obstacle is at one such antinode, the no-flow boundary conditions on the obstacle vertical boundary are naturally satisfied. The presence of the obstacle thus leaves the flow unchanged. Formulated in terms of the reflection and transmission coefficients, this is simply written ast +r = 1. This can be verified, for example, in Huang (2007), for the case of porous vertical barriers, or in Ursell (1947), for the case of impermeable vertical barriers covering some fraction of the distance between the surface and the sea bottom. If the assumption that the row is infinitely thin is removed, the relationr +t = 1 may be replaced with |r +t| = 1 (Linton 2011), with little changes in the calculations carried out in § 5. The properties represented by (2.6), (2.4) and (2.5) may be visualised geometrically, as shown in figures 2(a) (for a non-dissipative case) and 2(b) (in the general case). With a non-dissipative interaction,t andr can be visualised in the complex plane as forming two sides of a right-angled triangle, of which the hypotenuse is of unitary length, as seen in figure 2(a). In contrast, a dissipative case is shown in figure 2(b). Energy preservation implies thatt must lie inside the circle represented in the figure.

Pitching and surging obstacles
In the foregoing section, the obstacles were considered fixed. However, pitching and surging obstacles, such as rigid or flexible OWSCs, can be modelled in the same way. Let r d andt d be the reflection and transmission coefficient of an individual row, assumed fixed (d stands for diffraction). When the obstacles are oscillating, the reflected and transmitted waves are the sum of diffraction and radiation effects, the latter being the effect of the obstacle oscillating motion. Letĥ + rad andĥ − rad be the two coefficients such that for forced pitching or surging oscillations with complex velocity amplitudeΘ, the wave radiated backward has amplitudeĥ − radΘ , and the wave radiated forward has amplitudeĥ + radΘ . Because of the obstacle symmetry,ĥ − rad andĥ + rad verify thatĥ − rad = −ĥ + rad . Consider an incoming waveη 0 propagating forward from x = −∞. LetĤ + denote the transfer function that relates the incoming wave to the row oscillation velocity, i.e.Θ =Ĥ +η 0 . Under wave excitation, the reflected wave thus has complex amplitude (r d −ĥ + radĤ + )η 0 , while the transmitted wave has amplitude (r d +ĥ + radĤ + )η 0 .

A recursive numerical approach
In this section, the case of several rows is considered. As illustrated in figure 3, each row n is located at a position x n along the x axis, and is characterised by its transmission and reflection coefficients,t n andr n , respectively.
Define X 0 , X 1 , . . . , X N such that X 0 < x 1 , x N < X N , and, for every n such that 0 < n < N, x n < X n < x n+1 . In the fluid domain surrounding X n , and delimited by the previous and next rows, the wave field is specified by means of two coefficientsη + n andη − n , representing 'forward' and 'backward' complex wave amplitudes, such that For n = 0 (respectively, n = N), the domain where (3.1) is valid extends towards −∞ (respectively, +∞). Furthermore, if incoming waves propagate from −∞ in the positive x direction,η + 0 represents the forcing term of the whole system, whileη − N = 0, i.e. there is no incoming wave propagating from +∞.
Consider a section S of the array comprised between X and X , with one or more obstacles located inside the interval [X; X ]; S is represented in figure 4. Define the Figure 3. Wave transmission and reflection through N obstacles. 'forward' and 'backward' complex wave amplitudesη + andη − such that, in the neighbourhood of X, the FSE is written as follows: Similarly, define the complex wave amplitudesη + andη − such that, in the neighbourhood of X , the FSE is written as follows: Assume that the reflection and transmission problems have been solved for the domain S = [X; X ], i.e. that complex coefficients R − , R + , T + = T − = T have been found, such that, for incident wave componentsη + andη − propagating into the domain S, the wave componentsη − andη + , propagating away from S, are derived as follows: (3.4) In the above expression, the 'forward' and 'backward' transmission coefficients, T + and T − , are assumed identical, which will receive proper justification subsequently. Similarly, let S be the domain comprised between X and another position X > X , and represented in figure 4. Complex coefficients R − , R + , T have been found, such that, for incident wave componentsη + andη − , the wave componentsη − andη + are determined as follows:η (3.5) Now define S as the domain extending from X to X , as represented in figure 4. By combining (3.4) and (3.5), it is straightforward to obtain a linear relation between {η + ,η − }, on the one hand, and {η − ,η + } on the other hand, such that where the coefficients R + , R − and T are calculated as follows: Note that, in the above expression, if both S and S satisfy the condition that their forward and backward transmission coefficients are identical, the same holds for S . The recursion of (3.7) can be put into practice by defining the elementary domains S n for n = 1, . . . , N, only containing the nth row, and extending from X n−1 to X n . For one such domain, it is easy to find that the transmission and reflection coefficients are as follows: The recursion can be initiated using the elementary domain containing the first row, thus obtaining T 1 , R + 1 and R − 1 . Additional rows are then added sequentially through the use of (3.7). It can be seen that, for every elementary domain, the forward and backward transmission coefficients are identical, and that this property is preserved through successive iterations of (3.7), which justifies a posteriori the simplifications in (3.4) and (3.5).
The above recursion yields the global transmission and reflection coefficients of the array, denoted as T {1,...,N} , R + {1,...,N} and R − {1,...,N} . For a given incident waveη + 0 , the wave field between each pair of consecutive rows, i.e.η + n andη − n , can be reconstructed. To that end, the reflected wave coefficientη − 0 is first calculated asη Then, other 923 A2-8 coefficientsη + n andη − n for n > 0 are calculated recursively using the relationŝ This completes the description of the plane waves across the whole array.
In the above developments, it is assumed that, behind the last row, waves are free to propagate away from the array (η − N = 0). However, it can be of interest to investigate the case where waves are blocked by some reflecting structure (e.g. the coastline) behind the last row. In such a case, it is easy to simply add one more row, characterised by a transmission coefficient equal to zero, and follow the very same procedure.
An important extension of the proposed approach consists of modelling some dissipation within the fluid, by defining with κ defined as κ = k + iν, where ν > 0 is some dissipation factor. The inclusion of such a dissipation factor is crucial in the analysis of small-scale experimental data, as reported in § 4. It can easily be shown that, by choosing X n = x n+1 , the recursion of (3.12) coincides with the one proposed by Evans (1990). Finally, from an array optimisation perspective, the proposed 'concatenation' formulae of (3.7) are advantageous, since they allow the combination of arbitrary array sections.

Experimental procedure
The proposed recursive numerical approach is validated by means of an experimental data set, recorded by Nové-Josserand et al. (2019). The set-up consists of a small-scale wave flume (0.6 m wide, 1.8 m long), with a wave maker at one end and a sloped damping beach at the other end. The water depth is 8 cm. Flexible, vertical blades are placed in the wave flume, arranged in one, two and three rows extending along the width of the wave flume. In these experiments, the blades are surface piercing (although this is not necessary to employ the theory of § § 2 and 3). A three-dimensional map η(x, y, t) of the FSE is recorded by means of an optical measurement system, as reported briefly in Appendix A, and in more detail in chapter 3 of Nové-Josserand (2018). All three-dimensional mappings of the FSE are averaged along the y axis, i.e. along the flume width, so that each experiment yields a signal η(x, t). All experiments are carried out using the same excitation frequency (5 Hz In each longitudinal position x, the η(x, t) signal is Fourier filtered at the 5 Hz excitation frequency, in order to cancel parasitic frequencies due, for example, to undesired transient events. An example of experimentally recorded signal η(x, t) is shown in figure 6 (before and after Fourier filtering). Marked reflection patterns are visible in the up-wave zone.
The experimental validation procedure relies on the following steps.
(i) Estimate the reflection and transmission coefficientsr andt of a single row of obstacles, from the results of a single-row experiment. (ii) Estimate the reflection and transmission coefficients R 2 , T 2 , R 3 , T 3 of twoand three-row arrays, for different configurations in terms of spacing between consecutive rows. (iii) Compare R 2 , T 2 , R 3 , T 3 with their theoretical counterparts, computed using the wide-spacing recursive model proposed in § 3.
However, before the steps outlined above can be followed, several practical challenges must be overcome. First of all, because of the small size of the wave flume, wave dissipation is expected to have a significant impact on observed results. Therefore, the results of a blade-free experiment are used to calibrate the coefficient ν of an exponential dissipation law as in (3.13). The procedure to calculate ν is reported in more detail in Appendix A. Furthermore, the observed wavelength (λ = 6.90 cm) departs slightly from that expected from linear wave theory (λ = 6.67 cm) for capillary-gravity waves (see Appendix A). Accurate measurement of λ was found to be highly important, especially in the process of estimating the phase of reflection and transmission coefficients. To that end, the spatial autocorrelation a( x) = E[η(x, t)η(x + x, t)] was calculated in the blade-free experiment, and λ = 6.90 cm was obtained as the distance between successive maxima of the autocorrelation function (more specifically, taking the distance between five consecutive maxima allows for excellent accuracy, below 1 %.) In addition, bearing in mind that surface tension might vary slightly across experiments, depending on impurities at the water surface, it was verified that λ remained consistent across all four sets of experiments.
Finally, the recorded signal η(x, t) must be decomposed into two free wave componentŝ η + andη − propagating in the positive and negative x directions, respectively. Those two components are required in order to estimate the array reflection and transmission coefficients. In this work,η + andη − are estimated based on a least squares method, thoroughly detailed in Appendix A.

Experimental results
The magnitude and phase of the transmission coefficient of a single row of flexible blades are found to be |t| ≈ 0.73 and φ t ≈ 0.1 rad (see in Appendix A how uncertainty in λ and ν is accounted for in the estimate oft). Note that, given those values, the row of blades is a 'dissipating' obstacle, since clearly |t| 2 + |r| 2 = |t| 2 + |1 −t| 2 ≈ 0.61 < 1 (see Appendix A). The missing energy is dissipated either internally (through bending) or through viscous effects.
In the set of experiments with N = 2 rows, the first row (i.e. the farther up-wave) is held fixed, while the spacing L between both rows is gradually increased by changing the second row position; L thus covers a range of discrete values, between approximately 0.25 and one wavelength. parameter uncertainty is also visualised, by plotting the model reflection and transmission coefficients, obtained when |t| is changed by ±5 % (thin lines, below and above the nominal results). Taking into account the significant sensitivity of the model results to the identifiedt magnitude, the agreement between experimental and predicted coefficients can be considered excellent, although transmission seems to be slightly overestimated by the model. In the experimental data set with N = 3, the length L 12 between the first two rows is held fixed to approximately 0.25λ, while L 23 is gradually increased. As evidenced by figure 8, the agreement between the model and experimental results remains satisfactory, although the model slightly underestimates reflection.

The particular case of identical rows
In this section, the case of identical rows is considered, whereby every row has the same transmission and reflection coefficientst andr, and the distance between each pair of consecutive rows is L. Indeed, in such a case, analytical formulae can be found for the array reflection and transmission coefficients (Evans 1990;Linton 2011). Several important results arise, regarding the effect of L upon the array reflection properties, such as the existence of discrete values which completely cancel reflection, and the existence of 'band-gaps', i.e. intervals where reflection is high, and tends to unity as the number of rows increases. Furthermore, it has been pointed out (Linton 2011) that, for some cases such as rows consisting of submerged, horizontal cylinders, spacing values which yield maximum reflection (i.e. the centres of the reflection band-gap intervals) are shifted with respect to the Bragg condition. In Arnaud et al. (2017), which investigates a dense array of wave-dissipating vertical cylinders, treated as a porous medium, it can be noted that the Bragg-type reflection minima are not null (in contrast to cases, see e.g. Rey et al. (2018), where no dissipation is considered). In the present section, those results are revisited and clarified by adopting a formalism, which highlights the role of energy absorption or dissipation by the obstacles.
Define X 0 = x 1 − L/2 and, for n > 0, X n = x n + L/2. For every domain S n , the transmission and reflection coefficients are given as where the 2 × 2 matrix A has complex entries defined as follows: Therefore, η 0 and η N are related as follows: However, in the equation above, only η + 0 and η − N are known, where η + 0 corresponds to the system input, while η − N is zero. Manipulating (5.5), it is easy to find the following relations: Therefore, the global transmission and reflection coefficients are calculated as where (A N ) m,n denotes entry n, m of the matrix A N . Of course, the solution above may be computed numerically. However, it is also interesting to investigate the behaviour of R N and T N qualitatively, depending on the basic parameterst,r and L, in order to gain more general insight. In general, any pair of coefficientst,r satisfying (2.6) may be constructed from a 'non-dissipative' pair of coefficientst ,r , as follows:t =t − δ, where |δ| = 1/2 − |t − 1/2| and 2φ is the phase oft − 1/2 (and φ is the phase oft ). Coefficientst,r,t ,r , δ are illustrated in figure 9. Essentially,t is the transmission coefficient which is the closest tot, while representing a non-dissipative obstacle. It is easy to show thatt andr verifyt = cos φe iφ , r = −i sin φe iφ . (5.10) With the decomposition of (5.8), the entries of A may be rewritten as a function of L, |δ| and φ, using the following equations: as well as the new variable θ = kL + φ.
The matrix A may now be reformulated as follows: (5.13) In order to obtain an analytical solution for the global transmission and reflection coefficients, A N should be calculated explicitly as a function of the model parameters. Instead of the recursion in Evans (1990), we can also proceed through a diagonalisation of A, or use the direct result proposed by Williams (1992) for powers of 2 × 2 matrices. Let ρ 1 and ρ 2 be the two eigenvalues of A. Then (5.14) The two eigenvalues of A are found to be It is easily verified that ρ 1 ρ 2 = 1. Adopting the same idea as Evans (1990), define α ∈ C such that cosh(α) = f (θ)/f (φ), and β ∈ C such that cosh(β) = g(θ)/g(φ). Then, ρ 1 = e α and ρ 2 = e −α . Using (5.14), we finally use the fact that f (φ) sinh(α) = g(φ) sinh(β) to obtain relatively simple expressions (Evans 1990) for R N and T N as follows: , it is easy to see that α(θ + π) = α(θ) + iπ and β(θ + π) = β(θ) + iπ. Hence, R N (θ + π) = −R N (θ ) and T N (θ + π) = (−1) N T N (θ ). As far as the reflection and transmission magnitudes are concerned, it is thus sufficient to study R N and T N on any interval of length π.
The mathematical characterisation of |R N | can be appreciated in the different graphs of figure 10 as follows.
It is interesting to see that maximum reflection is obtained for θ within intervals of the form [nπ − φ; nπ + φ], which corresponds to kL ∈ [nπ − 2φ; nπ]. Those intervals are wider as φ increases, i.e. as the reflection coefficient of individual rows predominates over their transmission coefficient. Furthermore, as N increases, reflection tends to unity over the whole band-gap width, while it remains bounded from above by the function sin φ/ sin θ outside band-gap intervals. Finally, band-gaps are centred around  Figure 10. Reflection coefficient magnitude for a set of N equally spaced rows of obstacles with no energy dissipation (|δ| = 0), with N = 2, 3, 4, 5, 10 (a-e,g-k), for φ = π/8 ('transmission-dominated') and φ = 3π/8 ('reflection-dominated'). The two parameter options forr,t of individual rows are represented in panels ( f,l). See also figure 9 for the meaning of φ and |δ|. kL = nπ − φ, which is in contrast to 'classical' Bragg peaks which occur for kL = nπ, as pointed out by Linton (2011). Therefore, strictly speaking, those two conditions only coincide if the phase of the transmission coefficient is small. However, given the location of the band-gap interval as specified above, the Bragg condition kL = nπ is never outside, but always at the limit of the band-gap interval (even for large φ).

Dissipative obstacles
The case where dissipation occurs is more complicated, since the parameter |δ| is now different from zero. However, it is easy to find the asymptotic behaviour of R N and T N as N tends to infinity: while T N tends to zero, R N tends to R ∞ = e −β .
The behaviour of |R N | is illustrated in figure 11 for N = 2, 3, 4, 5, 10 as well as |R ∞ |, for 'transmission dominated' cases (panels (a-g)) and 'reflection dominated' cases (panels (h-n)). In both sets of graphs, two different values of |δ| are shown: |δ| = 0.1 (weak dissipation, lighter colour) and |δ| = 0.4 (strong dissipation, darker colour). As can be seen in figure 11, |R N | oscillates around |R ∞ |. There is a larger number of oscillations as N increases, and oscillations are less marked as |δ| approaches its maximum possible value of 0.5. The figure also suggests that, even for relatively small N, the asymptotic function |R ∞ | is an acceptable approximation for |R N |.
This time, there are no values of θ, for which |R N | vanishes. It is interesting to see that a similar result is obtained in the case of an array of cylinders, as studied by Arnaud et al. (2017). Therefore, it seems to hold rather generally that, in the case of obstacle arrays made of identical rows, the presence of absorption or dissipation by individual rows makes it impossible to obtain zero reflection. One of the consequences is in terms of wave energy absorption: absorbing the totality of the incoming wave energy would imply cancelling both reflection and transmission coefficients. Therefore, total absorption of the incoming wave energy is impossible with an array of identical rows.
In the non-dissipative scenarios previously studied, the transmitted wave amplitude |T N | could simply be calculated using the fact that |T N | 2 + |R N | 2 = 1. In the present case, some part of the wave energy is dissipated or absorbed by the obstacles, so that |T N | cannot be 923 A2-17  Figure 11. Transmission coefficient magnitude for a set of N equally spaced rows, with N = 2, 3, 4, 5, 10 and N → ∞, and four parameter options regarding φ and |δ|. In each graph, the lighter and darker lines correspond, respectively, to |δ| = 0.1 (weak energy absorption or dissipation) and |δ| = 0.4 (strong energy absorption or dissipation). The panels (af ) and (h-m) are for 'transmission-dominated' and 'reflection-dominated' interactions, respectively. The four parameter options are represented in panels (g,n). See also figure 9 for the meaning of φ and |δ|.  Figure 12. Transmission coefficient magnitude for a set of N equally spaced rows, with N = 2, 3, 4, 5, 10, and four parameter options regarding φ and |δ|. In each graph, the lighter and darker lines correspond, respectively, to |δ| = 0.1 (weak energy absorption or dissipation) and |δ| = 0.4 (strong energy absorption or dissipation). The panels (a-e) and (g-k) are for 'transmission-dominated' and 'reflection-dominated' interactions, respectively. The four parameter options are represented in panels ( f,l). See also figure 9 for the meaning of φ and |δ|.
immediately deduced from |R N |. Therefore, it is interesting, for the sake of completeness, to also illustrate |T N |; figure 12 shows |T N | for the same range of cases as figure 11 (except for N → ∞ since T ∞ is simply zero). As can be appreciated in the figure, |T N | also exhibits some oscillatory behaviour. Unsurprisingly, |T N | fades out to zero more rapidly for stronger dissipation (larger |δ|), and when reflection of individual rows is larger than their transmission. The asymptotic function |R ∞ | deserves particular attention. To that end, figure 13 shows |R ∞ | for a wide range of possible scenarios regarding φ and |δ| (illustrated in the circles on both sides of the figure). As can be appreciated in figure 13, for |δ| > 0, the asymptotic function |R ∞ | shows a marked peak, reaching unity for distances L such that kL = nπ. This corresponds to classical Bragg peaks. Unsurprisingly, the peak width is smaller for larger |δ| (i.e. individual rows tend to dissipate more energy) and for obstacles for which transmission is predominant. It is also interesting to note that, even for relatively small (but non-zero) |δ| values, a strong contrast can be observed with respect to the band-gaps associated with 'non-dissipative' obstacles. We now briefly turn back to the problem of wave energy absorption. Given that T ∞ = 0, the energy absorbed by the whole array is proportional to 1 − |R ∞ | 2 , and thus larger for smaller |R ∞ |. Figure 13 allows the visualisation of thoser,t parameters, which are more favourable from an energy absorption perspective. It appears that stronger transmission by individual rows (e.g. the top-left graph) is desirable, so thatt should be chosen close to unity; with such parameters, the array absorbs most of the incoming wave energy, if the inter-row spacing does not match the Bragg condition. However, the particular casê t = 1,r = 0 merely corresponds to an absence of obstacles (and hence yields no energy absorption at all).
Finally, the limiting case |δ| = 0.5 (i.e.t =r = 0.5) corresponds to the reflection/ transmission characteristics of a single row of wave energy devices operating under optimal power-maximising control (Evans 1981). Nevertheless, as evidenced by the top-left graph, the choice |δ| = 0.5 is suboptimal for power absorption by the array as a whole. In other words, parameter tuning may be optimal for individual rows, but suboptimal for the interacting array. The need for model-based optimisation, articulating individual row parameters together with array interaction, is thus confirmed.

Conclusion
With a view to the modelling and optimisation of wave energy farms, a simple recursive formulation has been employed, based on the wide-spacing approximation, to solve for the transmission and reflection of plane water waves through a number of rows of thin, vertical obstacles, where each row is characterised by complex reflection and transmission coefficients. The method allows the combination of any pair of adjacent subdomains of the fluid, for which the reflection and transmission coefficients have already been solved, thus lending itself to efficient numerical solution. Additional effects, such as linear dissipation along a wave flume, or partial reflection by the shore or an absorption beach, may be easily incorporated in the analysis. The theoretical limitations of the model, in terms of applicability, are related to the array periodicity in the direction, transverse to the incoming wave direction, which must be smaller than the wavelength, and to the distance between consecutive rows, which must not be too small relative to the wavelength.
The array reflection and transmission, predicted by the recursive model from the reflection and transmission coefficients of individual rows, are compared with the results of small-scale physical experiments comprising N = 2 and N = 3 rows of vertical, flexible blades. There is a close agreement between experimental and theoretical results, even though the distance between consecutive rows is at times significantly smaller than the wavelength, which corroborates the results of several numerical studies (McIver & Bennett 1993;chapter 6 of Linton & McIver 2001;chapter 3 of Li 2006), as well as experimental results (e.g. Liu, Li & Zhu (2016), dealing with submerged semicircular bars). Those conclusions suggest that the approach can be considered a suitable modelling option for further simulation and optimisation of OWSCs.
In the case of identical, regularly spaced rows, distinctive behaviours are identified, depending on whether individual rows predominantly reflect or transmit waves, and whether or not they dissipate energy. With no dissipation, the well known fact that discrete values of the row-to-row distance L cancel reflection is retrieved and clarified. In each interval of width λ/2, where λ is the wavelength, there are exactly N − 1 values of L which cancel reflection, and those values are all located within an interval, with width governed by the reflection/transmission balance of individual rows. The existence of band-gap intervals is also retrieved, centred away from the classical Bragg resonant condition. In contrast, when dissipation is present, the reflection magnitude |R N | is always non-zero and, as L grows, oscillates close to a reasonably simple asymptotic function R ∞ , which is interesting to investigate. In contrast to the non-dissipative case, reflection maxima take the form of the well known Bragg peaks, reaching unity when L is a multiple of λ/2. Such peaks are narrower as more dissipation takes place. Finally, from the point of view of power absorption, results illustrate the fact that design parameters may be optimal for individual rows, but suboptimal as soon as array interaction is considered, thus confirming the need for joint optimisation of device and array characteristics.
Regarding the objective of wave energy farm modelling and optimisation, further work is still required towards a sufficiently detailed mathematical model of flexible OWSC farms. In particular, it should be noted that reflection and transmission coefficients of individual rows are determined by both the obstacle geometry, and by their dynamical behaviour, including the possible use of control. Array optimisation of flexible OWSCs will require the distinct articulation of geometry and control effects, since those may have significantly different implications in terms of investment cost.
In addition, ocean waves are best described as Gaussian, random processes (Pierson 1952), as opposed to the monochromatic case treated in the present study. Such random waves are statistically characterised by their spectral density function, which shows how energy is distributed across frequencies. It will be interesting to assess how the array affects the wave field energy content, both up-wave and down-wave -the latter being particularly important for coastal protection purposes. To that end, the methodology employed in this paper can be readily generalised to polychromatic (random) wave excitation signals. It is sufficient to know the reflection and transmission coefficients of individual rows of obstacles over a range of frequencies covering the incident wave spectrum bandwidth. The array reflection and transmission coefficients then become complex functions of the frequency, and every variable of interest (reflected or transmitted spectrum, absorbed power, obstacle motions, and so on) can be analysed in terms of dynamical systems driven by Gaussian inputs. In particular, the tuning of row-to-row spacing and other design parameters will need to be chosen, taking into account typical wave spectra at the location considered, not just one frequency.
wave maker to the beach, and back to the observation zone where the blade array is located. Similarly, the analysis time window starts sufficiently late so that transient dynamics have faded out. In this work, the analysis time window of each experiment has a duration of 2 s, which covers 10 wave periods. The spatial zone of analysis, illustrated in figure 6, is restricted to a 40 cm long area, encompassing the Lego base.
Finally, the η(x, y, t) mapping is averaged along the y direction, thus obtaining a two-dimensional η(x, t) mapping, as illustrated in figure 6. Here η(x, t) is Fourier-filtered at the 5 Hz excitation frequency, to obtain a single position-dependent complex amplitude coefficientη(x), which characterises the wave field.

A.2. Wavenumber and dissipation coefficient estimation
For waves at a 5 Hz frequency, surface tension effects are not negligible anymore, so that the theoretical dispersion relation should follow that of capillary-gravity waves, where σ = 0.074 N m −1 and ρ = 1000 kg m −3 are the water surface tension and density, respectively. As mentioned in § 4, the observed wavelength (λ = 6.90 cm) departs slightly from that expected from the above theoretical relation (λ = 6.67 cm). Such a difference may be attributed to imperfect cleanliness of the water surface, which may have increased the surface tension value. Alternatively, the deviation may have been due to the presence of a relatively slow Stokes-induced current (less than 1 cm s −1 ), which would have resulted in a Doppler effect. Those possibilities are the subject of ongoing investigation by the authors.
Here |η(x)| is smoothed using a moving average filter, thus yielding the blue curve in figure 14 (shown in logarithmic scale). The best fit is found for a dissipation rate ν = 2.15 m −1 , also shown in figure 14, along with the corresponding lines obtained when ν is changed by ±5 %. Note that the assumption of a homogeneous dissipation coefficient along the wave flume is somewhat a simplification. In fact, the figure suggests that dissipation might be stronger in the central area (15-25 cm into the analysis zone), where the Lego base is located.
In all subsequent experiments, with blades in the small-scale flume, the equivalent wavenumber is considered to be κ = 2π/λ + iν, where λ and ν are the measured wavelength and dissipation rate.

A.3. Transmission and reflection coefficient estimation
The experimental identification of an array's reflection and transmission coefficients requires the separation of the forward-and backward-propagating wave amplitudes in a given observation zone, which is not a simple matter. In this work, a least squares method is employed, based on the idea introduced by Mansard & Funke (1980). the zone is the sum of a forward-propagating waveη + and a backward-propagating wavê η − , such that in every position x m ,η( (A2a-c) Findingη + andη − is achieved by solving the equation Ay = z in a least-square error sense, i.e.ỹ = (A * A) −1 A * z, where * denotes conjugate transpose. The same method is applied in the down-wave zone with a set x 1 , . . . , x M of down-wave positions and a reference position x 0 , to obtainη + andη − , the forward and backward wave coefficients in the down-wave zone. As pointed out by Mansard & Funke (1980), this procedure exhibits some sensitivity to the choice of measurement points x m and x m , but less so than by using only two measurement points and proceeding through a direct inversion of a 2 × 2 system.
Once forward and backward wave amplitudes have been estimated in the up-wave and down-wave zones, the reflection and transmission coefficients of any array located between positions x 0 and x 0 can be calculated as follows: However, as mentioned above, results are sensitive to the choice of measurement points x 1 , . . . , x M and x 1 , . . . , x M . Therefore, similarly to the procedure followed by Nové-Josserand et al. (2019), a set of 30 positions x 1 , . . . , x M=30 is randomly selected in the up-wave interval, as well as in the down-wave interval, andR andT are estimated based on those two sets of random points. The same procedure is iterated 100 times (i.e. with 100 different choices of 30 random locations), and the averageR andT coefficients are taken as the final estimates. The up-wave interval is between 0 and 100 mm inside the frame, while the down-wave interval is between 250 and 350 mm inside the frame. The identification ofr andt for a single row is of particular importance. This was done by employing the procedure outlined above, and choosing both x 0 and x 0 equal to the row position. The results of the identification procedure are illustrated in figure 15, where the 100 different pairs of coefficients are shown in the complex plane as small markers, along with their average as larger, circular markers. It can be seen that the sum of the estimated t andr identified is reasonably close to unity, as expected from (2.6) (see also figure 2).
However, the phase oft is found to exhibit high sensitivity to the exact wavelength estimate λ, while that ofr depends strongly on the exact referencing of the row position in the recorded frame x 0 . In contrast, the norm of both coefficient estimates is almost insensitive to those parameters. Therefore, it is preferable to confirm the estimates above through the use of an alternative method. This is done by estimating only |r| and |t| (which are almost insensitive to λ and x 0 ), using the approach outlined above, and using the fact thatr +t = 1. More specifically, in each random choice of M = 30 random points up-wave and down-wave, the norms |r| and |t| are calculated by taking the norm of (A3). The procedure is iterated 100 times, which yields reliable estimates of |r| and |t|. Finally, usingr +t = 1, it is found that Among the two possible values for the phase oft, the negative one is more consistent with the observations in figure 15, and is therefore chosen. The resulting estimatet is indicated by means of a square in figure 15, very close to the one obtained through the direct estimate (average of the green crosses in figure 15).
Finally, it is worth noting that the row of blades is a 'dissipative' obstacle, since clearly |t| 2 + |r| 2 < 1 (in other words,t andr are not located on the circle pictured in the figure). The energy dissipation splits into external dissipation (through friction and drag) and internal dissipation (through the blade internal damping). Only the latter may be converted into useful energy. However, in this work, we do not investigate the balance between external and internal dissipation.