Assessment of the Power Output of a Two-Array Clustered WEC Farm Using a BEM Solver Coupling and a Wave-Propagation Model

One of the key challenges in designing a Wave Energy Converter (WEC) farm is geometrical layout, as WECs hydrodynamically interact with one another. WEC positioning impacts both the power output of a given wave-energy project and any potential effects on the surrounding areas. The WEC farm developer must seek to optimize WEC positioning to maximize power output while minimizing capital cost and any potential deleterious effects on the surrounding area. A number of recent studies have shown that a potential solution is placing WECs in dense arrays of several WECs with space between individual arrays for navigation. This innovative arrangement can also be used to reduce mooring and cabling costs. In this paper, we apply a novel one-way coupling method between the NEMOH BEM model and the MILDwave wave-propagation model to investigate the influence of WEC array separation distance on the power output and the surrounding wave field between two densely packed WEC arrays in a farm. An iterative method of applying the presented one-way coupling to interacting WEC arrays is used to compute the wave field in a complete WEC farm and to calculate its power output. The notion of WEC array ‘independence’ in a farm from a hydrodynamic point of view is discussed. The farm is modeled for regular and irregular waves for a number of wave periods, wave incidence angles, and various WEC array separation distances. We found strong dependency of the power output on the wave period and the wave incidence angle for regular waves at short WEC array–array separation distances. For irregular wave operational conditions, a large majority of WEC array configurations within a WEC farm were found to be hydrodynamically ‘independent’.


Introduction
Ocean wave energy is a promising source of clean electricity that has the potential to make a significant contribution in reducing the world's dependence on fossil fuels. However, in order for it to follow the path of offshore wind and become a commercially viable power source, significant cost reductions must be made. Because of physical restrictions on the size of individual wave energy converters (WECs), it is the established view of the wave-energy community that WECs have to be deployed in farms to be economically viable. To benefit from developing offshore infrastructure and the maritime support industry, such farms need to have a power rating in the order of hundreds of megawatts. With the most promising current WEC technology, this corresponds to farms of hundreds of WECs. How these WECs are grouped and arranged within a WEC farm to maximize profitability while minimizing detrimental effects is still an open question. For a key group of WECs nearing commercial deployment, i.e., heaving axisymmetrical point absorbers, a number of recent studies have numerically and experimentally investigated the layout and spacing of WECs within WEC farms [1][2][3][4][5]. Although the terms "WEC farm" and "WEC array" are used interchangeably, we define a "WEC farm" as comparable in size to an offshore wind farm that may consist of a large number of sparsely separated WECs or clusters of densely packed WECs, which we hereby term "WEC array". All of the aforementioned investigations utilized potential flow theory, specifically the Boundary Element Method (BEM), to resolve intra-array effects, that is, those between the WECs in the array. While effective for arrays with a small number of bodies, BEM modeling becomes computationally demanding as the number of bodies and modeled frequencies increase.
We follow an alternative approach whereby a WEC farm comprising two WEC arrays is modeled using a one-way coupling technique between BEM model NEMOH [6] and wave propagation model MILDwave [7,8]. One-way coupling means that the perturbed wave field is only propagated from the inner domain to the outer domain, as evidenced in Figure 1. In our investigation, the inner domain is the WEC arrays' near field, while the outer domain spans the entire farm area and the far-field. We used the BEM model in the near-field area of the WEC arrays and the wave propagation model in the far-field WEC farm area external to the WEC arrays (see Figure 1). A key feature of the proposed one-way coupling technique is that waves are propagated from the near-field model domain (NEMOH) to the far-field model domain (MILDwave) via a transfer of information on a wave-generation circle at coupling radius r c . A schematic of these domains and the clustered layout is presented in Figure 1. Wave loading in NEMOH is determined by the wave conditions in the domain at the WEC array location. If the effect of one array on another is sufficiently small, then these disturbances in the wave field due to the interaction can be ignored; therefore the arrays can be simulated by using the same incident wave conditions. If the WEC arrays are sufficiently close, allowing for mutual hydrodynamic interaction, however, the effect of the perturbed (radiated and diffracted) waves from one array on another needs to be taken into account. Such an approach would, of course, require multiple simulations and would take longer to perform. The crucial question then is at what distance we can consider two arrays to be sufficiently hydrodynamically independent to model them as isolated. Depending on the wave type, wave incidence angle, and interarray separation distance, we explore the magnitude of the effect of the presence of one array in the proximity of another on the total WEC farm power output. Two closely spaced staggered arrays of nine-point absorber-type heaving WECs are modeled using the aforementioned coupling technique. We investigated various interarray separation distances for a range of wave-incidence angles for regular waves and various interarray separation distances in irregular waves. The power output for the different configurations of the WEC farm is calculated and compared to that of a WEC farm of hydrodynamically independent WEC arrays, i.e., those operating in isolation. The minimum interarray separation distance, D, for which two WEC arrays in a farm can be considered as hydrodynamically independent is defined for each simulated wave period. As our focus is on operational sea states, in this work we operate on the paradigm of linear potential theory, as detailed in Section 2.

Linear Potential Flow
This investigation assumes linear potential flow theory [9], a subset of linear wave theory that allows flow velocity v to be expressed as the gradient of the potential Φ (Equation (1) The assumptions underlying potential flow are the following: • the fluid is inviscid; • the fluid is incompressible; and • the flow is irrotational.
The standard assumption of linear theory that the bodies' motion amplitudes are much smaller than the wavelength also applies. Linear potential flow theory has hitherto been utilized in a majority of the investigations into WEC array modeling, for example, see References [2,5,10,11]. Due to the principle of superposition, linear potential theory allows for the separation of the total wave field into the following components (Equation (2)): where φ t is the total velocity potential, φ i is the incident wave potential, φ d is the diffracted wave potential, and ∑ 6 i φ r is the sum of radiated wave potentials for each Degree of Freedom (DoF) of the WEC. In our investigation, we only modeled the heave motion for simplicity and because heave is the primary operating DoF of the modeled WEC. We also introduced the term 'perturbed wave' to denote the wave resulting from the sum of the diffracted and radiated potentials.

Boundary Element Method Solver
In our coupling technique, the intra-array effects, induced by the hydrodynamic interaction between the WECs, are resolved by simulating the WEC motions using open-source potential flow BEM solver NEMOH [6]. Given Equation (1), NEMOH solves the Laplace Equation (3) for complex velocity potential φ: given a set of boundary conditions on the wetted body surface, the free surface, sea bottom, and far field. The motion equations of each WEC are solved in the frequency domain using the method of Green's functions, as explained in Reference [6]. NEMOH has been extensively used for modeling WECs and WEC arrays in recent years, and has been validated against other popular BEM solvers, such as WAMIT R [12]. An important restriction imposed by the method is the assumption that water depth h is constant throughout the near-field domain (NEMOH domain in Figure 1). Free surface elevation η is calculated by taking the real part of complex potentialη that is, in turn, obtained in NEMOH from free-surface boundary condition Equation (4). From the superposition principle of Equation (2), free-surface elevations η can be separately obtained for the WEC motions due to the diffracted and the radiated potentials.
where g is acceleration due to gravity, and z = 0 is the undisturbed water surface.

Mild-Slope Wave Propagation Model
For simulating the far-field effects, e.g., the 'shadow zone' in the lee of the array, wave propagation model MILDwave was employed [7,8] in the outer domain. MILDwave, developed at the Coastal Engineering Research Group of Ghent University, Belgium, is a phase-resolving model based on depth-integrated mild-slope equations (Equation (5a,b)) in the form proposed by Radder and Dingemans [13]. MILDwave has been used in modeling WEC arrays in a number of recent publications [8,[14][15][16][17]. The mild-slope equations (Equation (5a,b)) are solved using a finite-difference scheme that consists of a two-step space-centred, time-staggered computational grid, as detailed in Reference [18].
Here, η and φ t are, respectively, surface elevation and total velocity potential at the free water surface, g is the gravitational acceleration, C is the phase velocity, and C g the group velocity for a wave with wave number k and angular frequency ω.

Modeled WECs
The WEC type modeled in this study is a flat circular cylinder with a diameter of 10 m and a draft of 2 m. The shape was selected based on its overall dimensions being similar to several promising WEC technologies, namely, Seabased, Seatricity, and Carnegie Wave [19]. All three WECs are in the planning stages of a precommercial WEC array. The Power Take Off (PTO) of each WEC is modeled as a resistive damper with a B PTO value of 3.6 × 10 5 kg s −2 , which is representative for a resistive PTO of the WEC type we model [20] targeting a sea state with a peak period of 8 s. This would correspond to a mode of the wave-climate period encountered in parts of the North Atlantic where WEC array demonstration projects are in the planning stages. The natural or resonance period of the WEC, T r , is equal to 4.6 s and, for simplicity, the value of B PTO was set identical to each of the WECs in an array. Further details can be found in Reference [17].

WEC Array and WEC Farm Layout
To simulate a realistic array of WECs, we chose a staggered configuration that has been shown in a number of numerical studies of heaving WECs [21][22][23], to maximize power in both regular and irregular sea states. Similar results for staggered configurations were shown in experimental studies in References [3,4,15]. For each of the farm configurations, we simulated two nine-WEC arrays as shown in Figure 2 within the farm shown in Figure 1 at various interarray separating distances D 1 from each other. The array orientation was held constant, while the angle of the incoming waves relative to the x-axis, β, was set at 0 • , 45 • , and 90 • . A schematic of the farm layout is shown in Figure 1. In this investigation, water depth was held constant at 40 m.

Coupling of NEMOH to MILDwave
In order to model the far-field effects in an efficient manner and with reasonable accuracy, a one-way coupling method introduced in [8,15,17] is employed. In brief, the perturbed wave field is calculated in the BEM code NEMOH and is propagated into the depth-integrated wave model MILDwave along a circle large enough to enclose the near-field domain that surrounds the WECs. Based on the aforementioned analysis, we set the coupling radius r c at the smallest possible value which results in a discrepancy of less than 2% in |η| between NEMOH and MILDwave. For the present investigation the value of r c is set at 100 m. The MILDwave grid resolution is set at ∆x = ∆y = 1 m and the outside boundary conditions are Sponge Layers which are calibrated to minimize reflection [14]. Further details on the coupling are available in [17], which introduced the present NEMOH-MILDwave coupling method for WEC arrays.

Calculating the Total Wave Field of the Perturbed Sea State-Demonstration for a Regular Wave
To assess the effects of the two WEC arrays within a WEC farm on each other, and in order to evaluate the total power output of the WEC farm, we needed to calculate the total perturbed wave field in the MILDwave domain. As we assumed linear theory in our work, we could use the superposition principle to sum up the total wave field by combining an iterative approach with the coupling method presented in Section 3.3. The technique employed is illustrated in Figure 3. The initial step (Step 1) was to propagate the incident wave field in the empty numerical basin in MILDwave to obtain the undisturbed wave elevation. In Step 2, the incident wave field was used as input into NEMOH, whence the 1st order perturbed wave of WEC Array I, p 1i , was evaluated. In Step 3, the average wave amplitude at the location of the center of Array I, was used as input into NEMOH to calculate the 1st order perturbed wave of WEC Array II, p 1ii . In Step 4, the process in Step 2 was repeated, with p 1ii as the new input perturbed wave. Finally, in Step 5, the same process was performed for the 2nd perturbed wave of WEC Array I, p 2i . Since the input perturbed wave field in each subsequent step was reduced by approximately an order of magnitude, for all practical purposes this process could be terminated at Step 4 without any appreciable loss in accuracy, even for closely spaced cases where interaction is maximized. Therefore, Step 5 is only displayed for a complete description of the proposed coupling method.

Coupling Irregular Waves
In this paper, we modeled an irregular long-crested sea state using a nondirectional Pierson-Moskowitz spectrum S PM (H m0 , T p , ω) with N = 20 frequency components, which, according to analysis of existing work in Reference [24] is sufficient for WEC motion simulation. Wave amplitude A i of each irregular wave component is calculated as: where δω is the angular frequency increment. Total wave elevation η for an irregular wave field is then the sum of ζ i , and unit amplitude total wave η i obtained using the procedure in Figure 3 in Section 3.4 for each frequency component i:

Determining the Power Output of a Nine-WEC Array
To evaluate the influence of the interarray interaction effects on the performance of a WEC farm, we computed the total power produced by the two WEC arrays after having obtained the modified wave field in the WEC farm using the approach outlined in Section 3.4. For each WEC array, using the amplitude of the total modified wave field at the locations of the WECs as the input, we calculated the power output by simulating the WEC motions in NEMOH using Equation (8) for regular waves for each frequency component i.
Here, M is the number of bodies in the array, ω is the wave angular frequency, Z j is the complex heave Response Amplitude Operator (RAO) of WEC j, and B PTO is the PTO damping coefficient, set equal to 3.6 × 10 5 kg s −2 for each WEC. For modeling irregular wave cases, we modeled the power output as the sum of the power at each wave component frequency i calculated by Equation (8) weighted by the nondirectional Pierson-Moskowitz spectrum S PM ( f ): In Equation (9), ∆ω is the frequency bandwidth of the spectrum discretization and N = 20. The total power output of the WEC farm is the sum of the power produced by the two WEC arrays. For a WEC array in regular waves, the power output is modeled for three wave incidence angles β: 0 • , 22.5 • , and 45 • . In exploring the effect of irregular waves in a WEC farm, only the head on incidence angle β = 0 • was simulated. Although B PTO was set constant and at the same value for both regular and irregular wave cases, sensitivity analysis was performed with varying values of B PTO , which showed that the constant of 3.6 × 10 5 kg s −2 chosen for this paper results in the maximal power for one WEC for a variety of wave conditions.'

Wave Field around two WEC Arrays within a WEC Farm in Regular Waves
We begin our analysis by qualitatively looking at the coupled total wave fields of the two array WEC farms produced by the iteration method outlined in Section 3.4. Representative results shown in this section in Figures 4 and 5 are for a wave height of H = 2 m regular wave with head-on incidence angle β = 0 • for wave periods of T = 6, T = 8, and T = 10 s. Figure 4 shows |η| for an array separation distance of D 1 = 200 m, while Figure 5 shows the total |η| for D 1 = 600 m. We note that the wave field shown inside the circular region of radius r c of both WEC arrays is the wave field calculated in NEMOH, initialized by the average |η| given by MILDwave at the end of Step 4, as shown in Figure 3 in Section 3.4. We noticed a strong contrast between the results for T = 6 s and the other two simulated periods, namely, T = 8 and T = 10 s. For the former, both the magnitude and the extent of the disturbances in the wave field due to the presence of the array are quite notable. We also observed strong positive |η| anomalies on the y = 0 axis that were not present at T = 8 s and T = 10 s. The positive anomalies at T = 6 s were due to interference between radiation and the diffraction of the WEC array optimized for β = 0 • , as detailed in Section 8. The absence of the aforementioned anomalies for T = 8 and T = 10 s is explained by the fact that radiation is minimal and the observed |η| mainly represents diffraction. We observed that, in Figure 5 for the D = 600 m, behind the second array, there was a reduction in |η| (the so-called 'wake-effects') for T= 10 s and an increase for T = 6 s, with little change for T = 8 s. These results can be explained in part by the more favorable position of Array II with respect to Array I that enabled it to radiate more waves for the case of larger separation distance. Such disturbances in |η| can be correlated with the WEC farm's performance as is shown in Section 5.2. In terms of the magnitudes of the anomalies in the wave field, there was again a marked contrast between the case of T = 6 s and that of the other two periods. For the former, the positive anomalies reached a value of 1.35, meaning that, in places, the presence of the WEC arrays increased the undisturbed wave field by up to 35%. The same negative anomaly was not as strong, as the region of marked decrease in |η| was mitigated by the positive anomalies in |η| due to the radiation of WEC array II. In contrast, for T = 8 s and T = 10 s, the greatest positive effects were no more than 10%, and there was a symmetry in the values of the maximum positive and negative |η| anomalies that deviated around 10% from the mean value of |η|.

Power Output of a WEC Farm Composed of Two WEC Arrays in Regular Waves
In the next two subsections, we expanded on the qualitative observations made in Section 5.1 by quantifying the power output by a WEC farm composed of two WEC arrays separated by distance D 1 for incident waves of T = 6 s, T = 8 s, and T = 10 s. The procedure outlined in Section 4 was employed to calculate the power output of the two-array WEC farm for a range of separation distances D 1 . Total average power output is displayed in the graphs in Figure 6 for each period and wave incidence angle β. The thinner level lines are P isolated , the power output of a farm of hydrodynamically isolated WEC arrays, or 2× the power of a single nine-WEC array, while the thick lines represent P farm , the power output of the hydrodynamically coupled WEC farm. The results are also presented in a nondimensional manner in Figure 7, where P farm normalized by P isolated on the y-axis is plotted versus the nondimensional ratio of D 1 /d, where d is the WEC diameter. The ratio of P farm /P isolated is analogous to the q-value, a commonly used metric to assess array effects within individual WEC arrays [21].
We first noticed the oscillating nature of the power output, with values both above and below the line showing the power output of arrays operating in isolation. The oscillations decreased in magnitude as we moved the arrays away from each other. Observing the trend from Figure 6a to Figure 6c, we noted the absolute value of the power output decreasing with increasing period. This is an expected trend given the behavior of the disk-shaped buoy with resistive control in regular waves that maximizes the motion close to the resonance period, T r , of 4.6 s. Note also that, in addition to the decrease in P farm with wave period T, there is a slight decrease with increasing incidence angle β for T = 6 s and T = 10 s, especially in the case of the former. This is a consequence of the WEC arrays' shape, as seen in Figure 2, where an increasing intra-array shadowing on the second row of WECs for each WEC array was observed, as β increases toward 45 • , at which angle the WEC array effectively becomes aligned.

Wave Incidence at β = 0 •
In Figure 6, we plotted the power output for increasing separation distance D 1 between Arrays I and II for three incidence angles, β = 0 • (solid lines), β = 22.5 • (dash-dot lines), and β = 45 • (dashed lines). Figure 6a shows the result for T = 6 s, Figure 6b for T = 8 s, and Figure 6c for T = 10 s. Observe that the result for T = 6 s for β = 0 • shows the greatest power oscillations. This should come as no surprise, seeing that, in Figure 4a,b, there is a strong rapidly oscillating pattern of |η| in front of and in between the WEC arrays. Note also that, despite a single peak giving higher power output than the case of WEC arrays operating in isolation, the rest of the points fall well below the line of P isolated . This trend demonstrates that the optimized staggered WEC array configuration results in substantial power extraction from the incoming waves, and that when one WEC array shadows another, the effect is strongly negative. This deleterious effect on power of placing one WEC array in lee of another is mirrored in the results for T = 8 s (Figure 6b) and T = 10 s (Figure 6c).

Wave Incidence at β = 22.5 • and 45 •
In this subsection, we compared and contrasted the reposes of the WEC farm power output for WEC array off-axis wave incidence. The key message of the curves in Figure 4 is that, unlike the result for β = 0 • , the array off-axis graph shapes do not greatly vary across the three tested periods. In other words, the WEC farm exhibits similar behavior in power output across the three modeled wave peak periods. The power output of the WEC farm for β = 22.5 • is always higher than for β = 45 • , with the magnitude of oscillations about the P isolated also less for β = 45 • . We point out that, although for T = 6 s the power of the WEC farm at β = 22.5 • is generally lower than that for a head-on wave (β = 0 • ), this is not true for T = 8 s and T = 10 s. For T = 8 s, power output for β = 0 • and β = 22.5 • of both P isolated and (P farm ) is higher than for β = 0 • . For T = 10 s, P farm with β = 22.5 • is lower than the value of P isolated in head-on waves, but is higher than the result for P farm with β = 0 • . Again, we can link this result to the |η| plots in Section 5.1, where β = 22.5 • is generally in an area of positive interference and low variability compared to β = 0 • . Note that, for β = 45 • , P farm is consistently lower than for β = 0 • and β = 22.5 • . This outcome is explained in Section 8. Note that, by looking at Figure 7, the graphs converge toward unity as we increase the relative distance, a result expected from theory and presented in many studies, among them References [1,21,25,26]. However, from a practical point of view, it is important to remark that presenting the results in this manner hides the absolute difference in power. For example, for the closest separation distances for T = 6 s, the ratio P farm /P isolated is below 0.95 for β = 22.5 • and β = 45 • , and is greater than 0.95 for T = 8 s and T = 10 s. However, the improved relative performance comes at the cost of significant decrease in absolute power, as witnessed in Figure 6.

Quantifying Percent Difference between P farm and P isolated
While in Section 5.2 we explored the trends in WEC farm power production in absolute terms, in order to answer the question posed in the introduction of this paper, namely, that of the error introduced by assuming WEC array independence, we needed to quantify the percent difference between P farm and P isolated . We calculated the percent difference between P farm and P isolated for the three regular waves for the three wave incidence angles. As expected from analyzing power output in Section 5.2, the relatively large errors for T = 6 s β = 0 • stand out compared to the rest of the data. We see that the error was as large as 20% for D 1 = 700 m, and did not consistently decline below 10% until D 1 of 2000 m. This is a consequence of strong interference between the perturbed waves of the two arrays when they were aligned with the wave direction and with each other. For β = 22.5 • and β = 45 • , we noted that the percent error was below the 5% threshold for the former and 1% for the latter, meaning that the array effect played a minor role in modifying the behavior of P farm . Still, the trend was a decrease in the magnitude of the percentage difference, with an increase in array-array separation distance. To see this trend more clearly, we could plot the power in graphic format in Figure 8. Here, the decreasing asymptotic trend is obvious save for the anomalous result of β = 0 • for T = 6 s. The strong oscillations in the graph for D < 1000 m are a consequence of the resolution in x of our D 1 , where the strongly varying graph oscillating about P isolated is sampled frequently enough to capture the peaks and troughs of the perturbed waves of the two arrays. For D 1 > 1000 m, the lower resolution in x only shows the envelope of the trend. For the other eight cases, the frequency of the variability was not as strong; therefore, the curves look to be smoother over the entire range of separation distances.

Irregular Wave Results
In this section, we present the results for an irregular wave incident on the WEC farm in Figure 1. As mentioned in Section 3.5, the irregular waves in the study are modelled on a Pierson-Moskowitz spectrum with no directional spreading. The peak periods analyzed are T p = 6 s, T p = 8 s, and T p = 10 s, matching the period of the regular wave cases. Each irregular wave result is a weighted sum of the coupled wave field at each modeled frequency. As the dependency of the wave field on the incidence wave angle yields a similar pattern for irregular waves as for regular waves, we only present results for wave incidence β = 0 • . We begin, as in Section 5, by looking at the total coupled wave fields for the farm in Section 6.1 and then explore the WEC farm power output in Section 6.2.

Wave Field around Two WEC Arrays within a WEC Farm in Irregular Waves
The irregular wave results for T p = 6 s and T p = 8 s are plotted in Figures 9 and 10. The total wave field is obtained as the sum of the undisturbed (incident) and the perturbed wave field at 20 different frequency components. The chief difference we noted in comparing the irregular wave results in Figures 9 and 10 to the regular counterparts in Figures 4 and 5 was the decrease in the overall magnitude of the interaction, as would be expected for the case where wave energy was not concentrated at one frequency but was, instead, spread out. Moreover, we could observe an absence of significant areas of positive interactions, such as those encountered for a regular wave case in the top panels in Figures 4 and 5, in the bands surrounding the wake of the arrays, at approximately 20 • to 30 • off the y-axis. At the same time, the wake was quite strong, notably for T p = 8 s. For T = 6 s, we could explain the decreased wake by the ability of the WECs in the array to radiate, acting to immediately offset the decrease in the wave height in lee of the arrays. This contrast between the two modeled wave periods is starkest for the arrays separated by 200 m in Figure 9, where we observed a region of neutral or positive |η| from 200 to 800 m behind the arrays for T p = 6 s, and only negative anomalies in |η| for T p = 8 s. As is expected for multifrequency sea states, the interaction pattern witnessed in Figures 9 and 10 is more complex than that of the regular wave fields in Figures 4 and 5. Nonetheless, the majority of the array effect of our particular array configurations are in the region immediately on the array axis or slightly off it. This is the reason for the significant effects on power absorption for the arrays in a head-on sea state detailed in Section 5.2.1 and witnessed in Figures 6a and 8a. In the next section, we explore the power output for the analogous case in irregular waves for various array separation distances D 1 .

Power Output of a WEC Farm Composed of Two WEC Arrays in Irregular Waves
When we compare the result in Figure 11 to the corresponding regular wave case in Figure 6, we noted two large differences. The first was that for the case of an irregular wave, all of the array effects were detrimental to the power absorption of the WEC farm. This is further highlighted in Figure 12, where the P farm /P isolated ratio was below unity for all plotted nondimensional distances. Such was not the case for the regular wave for T = 6 s in Figure 6a, where there were certain distances for which the power output of the interacting WEC farm was greater than that of the isolated one. The second difference was that total power output for the three irregular wave peak periods was lower than for the corresponding regular wave periods. This is not surprising, as the efficiency of a heaving WEC decreased when the energy was spread out over many frequencies in the irregular wave case. We also observed that, unlike the regular wave case where the power output was highest for T = 6 s, for the irregular wave case power was highest for T p = 8 s. This difference in the behavior of the WEC farm could be due to the fact that for a Pierson-Moskowitz spectrum, which we used here to model the irregular waves, as the peak period increases, the frequency spectrum also narrows. Thus, the spectrum was widest for T p = 6 s and narrowest for T p = 10 s. Therefore, the difference between the energy bandwidth was greatest between T = 6 s and T p = 6 s. This effect is enough to decrease the performance of the heaving WECS by a factor of 3, and reverse relative power output vis-a-vis T = 8 s. Even if we did not explicitly model off-axis wave incidences β for irregular waves, we expected the results to mirror those for regular waves in that power output would be improved compared to the case of β = 0 • , but not much greater than unity. This can be easily seen in Figures 9 and 10, where the 'wake zone' extends out to ±15 degrees on either side of β = 0 • behind the WEC arrays.

Defining 'Hydrodynamic Independence' in a WEC Farm Composed of Two WEC Arrays
We have seen in Sections 5.1 and 6.1 that the various factors in play influencing the power output of a WEC farm lead to a very complicated pattern of interaction that can be hard to discern. It is natural, then, to ask how we can extract practical information from such data that can both serve to optimize the WEC farm layout for a specific goal, as well as to accurately calculate the wave fields around the WEC arrays. For this reason, we attempted to simplify the problem by quantifying the significance of the interactions by first setting the value of 5% as an 'independence' threshold. Consequently, we defined a WEC farm of two WEC arrays as hydrodynamically 'independent' if the power output was within ±5% of the power output by two independent WEC arrays that operate in isolation (the case of 2 × P array ). We recall here that in the hydrodynamically independent case, power output was computed for each WEC array in isolation, specifically that the undisturbed wave field is used as input for the motion equations of the WEC array. The power output for the case where there was interaction between the WEC arrays was determined by the iterative procedure outlined in Section 3.4. Here, the input wave field is the sum of the incident and perturbed waves from both arrays. For the case of irregular waves outlined in Section 3.5, wave field summation is performed over each frequency ω. In order to visualize this concept, we turn back to Figures 8 and 13, where we plotted the percent difference between P farm and P isolated for regular and irregular waves, respectively. Starting with the regular wave cases in Figure 8, we immediately observe that only for the case of T = 6 s and β = 0 • was the difference consistently greater than 10% for a range of separation distances D 1 . For the rest of the investigated regular cases, the difference was small, and, in fact, for T = 10 s only the β = 22.5 • waves resulted in a difference larger than 1% in power output. For T = 8 s, for all three wave incidence angles, the percent difference was below the 5% "hydrodynamic independence" threshold. We can therefore safely assume array 'independence' for an overwhelming majority of the regular wave cases presented in this study. For the irregular wave scenarios for β = 0 • in Figure 13, we see a slightly different but marked decrease trend in the difference between P farm and 2 × P array with increasing D 1 . While the case of T p = 6 s is still the 'worst' in terms of percent difference because of the frequency spread of the Pierson-Moskowitz waves, the percent difference for T p = 6 s was greater than that for T = 8. For T p = 10 s, the percent difference was less than 5% for all separation distances greater than 400 m. Although we did not model them in this investigation, based on Figures 9 and 10, we could surmise that, for the 'off-axis' wave incidence angles, the difference between 2 × P array and P farm would again be smaller. In summary, making the assumption of array 'independence' in a WEC farm, where the WEC arrays are modeled as isolated, is safe as long as one array is not directly in lee of another. Moreover, for D 1 greater than 1000 m, all modeled cases except for those of waves of T = 6 s and T p = 6 s were below the 5% threshold and could be deemed 'independent', allowing for a significant reduction in modeling complexity without loss of fidelity.

Discussion
In Section 7 we saw that the separation distance between the WEC arrays in a WEC farm is not the only factor that plays a role in determining the extent to which two WEC arrays are hydrodynamically linked. However, the asymptotic behavior of WEC array interaction with respect to array separation distance D 1 is evident in Figures 6, 8, 11, and 13. It should be noted that the extent of the separation distance that we have modeled is limited from a practical standpoint to 2800 m, and several studies [1,25] show that, in regular waves, two WECs can have an appreciable hydrodynamic influence on each other, even when they are separated by more than 5 km. Here, it is important to remark that, in those particular investigations and in others that showed similar strong interactions at large separation distances, the modeled WECs were optimally tuned to magnify WEC motions. Therefore, we should expect the perturbed waves from WECs tuned in such a manner to be greater than those of linearly resistively tuned WECs of the type that we modeled in this paper. Of the factors influencing the strength of both WEC farms, wave field modification and its power output, D 1 , had the largest influence. However, as we saw in Section 5.2, the period of the modeled wave has an appreciable influence on interarray iteration, especially close to the resonance period of the WECs. For irregular waves, as witnessed in Section 6.2, the difference between the three peak periods was not as strong given the frequency spreading inherent in the Pierson-Moskowitz spectrum we used in our model. We also demonstrated the linking of influence of the wave incidence angle and D 1 to the power output. Not only did the overall magnitude of the interaction effects decrease as the wave incidence changed from a β = 0 • heading to β = 45 • , but the variability over the range of D 1 decreased as well. This was a result of the relative position of the WEC arrays; when one array was not directly shadowing another, the likelihood of a decrease in performance of a WEC array located downwave was reduced. Consequently, for incidence angles |η| away from 0 • , the waves that interacted with WEC array II located downwave were closer to the undisturbed incident wave. Of note is the suboptimal performance of β = 45 • observed in Figure 6, specifically for the case of T = 6. The 'underperformance' of β = 45 • could be explained by the fact that, at this wave incidence angle, the staggered configuration became aligned, and the back row of the array was strongly shadowed by the front row, as can be witnessed in Figure 2. As the staggered configuration of the WECs became roughly aligned for the waves with β = 45 • , there was a significant 'wake effect' inside the WEC array, but not at the WEC-farm level. This is why there was also less oscillation in power output over the WEC farm separation distances D 1 for β = 45 • . These results remind us that, in order to construct an optimized WEC farm, both the design of the micro elements, i.e., the layout of the individual WECs in a clustered array, and the macro elements, that is the WEC farm layout composed of larger units, should be considered.
We should remark an important point about the trends seen in Figures 6 and 11. In particular, WEC farm interaction is beneficial to only a small subset of the regular wave cases modeled, and is never beneficial for irregular waves. This outcome is largely due to WEC type and the limitations of the linear resistive PTO modeled in this investigation. As was shown in a number of previous studies [5,21,26,27], one needs to implement active frequency-dependent control in order to fully take advantage of WEC motions to induce beneficial hydrodynamic interactions between WECs and, by extension, between WEC arrays. While we observed an overall decrease in the magnitude of interarray interactions as we increased the array separation distance consistent with the 1/ √ 2 asymptotic trend defined in Reference [25], there was significant difference in the smoothness of the power-output curve between the various tested wave periods and incidence angles β. It should be pointed out that the result was mainly due to the configuration choice of individual WEC arrays that were optimized for a certain incident wave direction, specifically, β = 0 • . Thus, when placed one behind another, the power output of the WEC farm substantially decreased.
The observed response of the WEC farm exactly mirrors the trend that was demonstrated for individual WECs placed at increasing intra-array distances from each other, such as in References [11,25,28]. In these papers, the net power in a WEC array trends to the sum of the power of isolated WECs as the separation distance becomes larger. In our investigation, we were able extend said observation to WEC farms composed of multiple arrays. Note that a similar conclusion was reached in Reference [28], where the authors separated a WEC farm into two clusters of WECs, concluding that offsetting array clusters so that one is not directly behind another is the best array layout design strategy. However, in Reference [28], the authors employed a BEM solver to simultaneously calculate all interactions, an approach that has limits as the number of simulated WECs increases. In contrast, our coupling method permitted us to model arbitrary large numbers of WECs, provided they were split into individual clustered arrays. As was noted in the introduction, this constraint could almost certainly be applied to WEC farms from practical and economical considerations. As we have shown in Section 7, unless WECs are closely spaced and are directly aligned with the incoming wave direction, such clusters can be assumed 'independent' with only a small error in the WEC farm power output estimate.

Conclusions
In this paper, an iterative coupling method between the near-field BEM solver NEMOH and far-field wave propagation model was applied to examine the WEC array interaction effects in a WEC farm composed of heaving resistive WECs. The method provides a robust and efficient means of calculating the wave field around compact WEC arrays and, in turn, allowed us to estimate the total power output of a WEC farm. Although the coupling gives accurate results to an arbitrary degree of precision, even a few orders of interactions require a complicated web of iterations as explained in Section 3.4 in Figure 3. Hence, it is natural to seek further simplification of hydrodynamic calculations. If we can assume that two WEC arrays (I and II) in a farm are hydrodynamically independent, i.e., they behave as isolated, then the power absorbed by each WEC array can simply be computed in one iteration. The total wave field in a farm could then be calculated as the sum of two perturbed wave fields generated by WEC Arrays I and II, where the motion of both arrays is forced only by the incident wave. We saw that the primary determinant for the power output of a WEC farm composed of linearly resistive heaving WECs for a given regular or irregular wave period is interarray separation distance D 1 . Nonetheless, wave incidence angle β plays a significant role in determining not only the total power output of a WEC farm for a given wave period but also the attenuation of the array effects with D 1 . It should be mentioned that, in this investigation, we focused on a narrow subset of modeling scenarios, namely, that the study was performed for heaving WECS with a linearly resistive PTO. Although we expect the same overall trends to hold for various classes of WECs, it is evident that, for actively controlled WECs that are able to be tuned for a particular sea state, WEC motion and, by extension, the perturbations in the wave field would be increased in magnitude and felt over a larger distance away from the array. It will be the topic of a future investigation to model a more realistic type of WEC, where each WEC's PTO is subject to active control. Finally, we should note that, although we have demonstrated the coupling technique of Section 3.4 and the power-output trends in Sections 5.2 and 6.2 for a WEC farm composed of only two WEC arrays, the method could easily be extended to WEC farms composed of multiple arrays. In this paper, we also showed that, for both regular and irregular waves, for a large majority of cases, two WEC arrays in a farm could be considered hydrodynamically independent for the purposes of assessing the power output of a WEC farm. In this case, a simple and fast coupling method consisting of only one summation for each array could estimate power production with high accuracy. In Section 7, we investigated the error magnitude that is introduced into the calculation by making the assumption of hydrodynamic independence of the WEC arrays. We observed that the error introduced by the array independence assumption was within 5% for all cases except for only the closest separation distances D 1 for T = 6 and T p = 6 for β = 0 • . As was noted in Section 6, in this work we did not explicitly model off-axis wave incident angles β. This work on short-crested irregular wave modeling of WEC arrays has been accepted for publication as of the date of this paper. Additionally, as we alluded in the discussion in Section 8, our results depend on the type of modeled PTO. To this end, our group is integrating various realistic PTO types into our coupled modeling via a PTO module. These results will be published in an upcoming series of articles. If we extend our scope beyond purely hydrodynamic considerations and consider the economic constraints of a commercial WEC farm, we can see that separating WEC farms into hydrodynamically independent clusters of WECs can have advantages beyond simplifying power-output calculations. By concentrating many WECs in close proximity, a WEC farm developer can save on the cost of marine cables that are known to be a significant expense item for offshore energy projects. Furthermore, spacing constraints, such as leaving navigation channels for operations and maintenance navigation and other sea users in between the WEC arrays, naturally result in a clustered layout for the WEC array. Based on the results obtained in this paper, we can consider such a WEC farm to be composed of 'independent' arrays and apply the present coupling method to the problem of technoeconomic optimization of WEC farms that is necessary for the commercial viability of wave energy.
Author Contributions: P.B. generated the idea of the paper, set up the numerical simulations, and wrote the manuscript. G.V.F. performed numerical simulations and assisted in the developing the conceptual model. V.S. and P.T. proofread the text and provided the funding resources.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

DoF
Degree of Freedom PTO Power Take-Off RAO Response Amplitude Operator WEC Wave Energy Converter β angle of incidence of the incoming wave to the x-axis ( • ) d x , d y intra-array WEC separation distances in the x and y direction (m) B PTO linear power-take-off damping coefficient (kg/s 2 ) D 1 interarray centre-to-centre separation distance (m) M number of bodies in the WEC array N number of frequencies in irregular sea state discretization η free-surface elevation (m) |η| absolute value of the wave amplitude η (m) p ij perturbed wave of order j for array i (-) P i (ω, β) mechanical power produced by the WEC for a given frequency and wave direction (kW) r c coupling radius (m) P isolated 2 × total power output of an isolated WEC array (kW) P farm total power output of a WEC farm (kW) T r resonance or natural period of an oscillating body (s) Z j complex amplitude of heave velocity of body j (m/s) ζ complex wave amplitude i (m) array effects = hydrodynamic effects of WECs in an array that produce a perturbation in the incident wave field perturbed wave = radiated + diffracted wave