Responses to Reviewer Comments on “Effect of tip spacing, thrust coefficient and turbine spacing in multi-rotor wind turbines and farms”

The authors employ large-eddy simulations (LES) to quantify the difference in wake effects between a four-rotor wind turbine and an equivalent single-rotor wind turbine, for both an isolated wind turbine and a row of five wind turbines. The present work is a continuation of a conference article, where more parameters are investigated. An engineering wake model is calibrated with the LES data and its performance for a row of five four-rotor wind turbine is investigated. The article is well written and provides interesting results. I have written a list of main and minor comments below:

approach and critical sense of scientific results.
Response: We have performed simulations at two further levels of refinement, labeled G4 and G5 in the revised manuscript. These simulations have been performed for 1-rotor as well as 4-rotor turbine with s/d = 0.05. The results show that the velocity deficits and added turbulent kinetic energy (TKE) values are better converged at these finer resolutions. It is further observed that the convergence of velocity deficit profiles is non-monotonic, a common observation in channel flow simulations (Meyers and Sagaut, 2007). The convergence of the resolved TKE is monotonic, as expected. Comparing the results of our working grid G2 with grid G5, we obtain differences in velocity deficit of 3.2% and 1.9% at x/D = 4 and 6 respectively. In view of the marginal difference in the velocity deficits and the enormous added cost (almost 16 times), we feel that the use of grid G2 for evaluating the effect of all problem parameters is justified. We would like to reiterate that our grid size is in keeping with several previous studies as outlined in our previous response. Also, as mentioned by the reviewer himself, it is unreasonable to expect one to perform a parametric sweep on a fully converged grid. We also agree with the reviewer about the need to highlight possible limitations of a study. We have included all of the above information in the revised manuscript (lines 173 -202). We hope this addresses the concerns of the reviewer and makes the article more useful to future readers.
-My previous comment was: "The authors state "It is seen that P2-5 is larger for all 4-rotor wind farms...". This is not correct. If you look at Figure 10(c) there is actually a cross-over for the 3rd turbine, where the single rotor produces more. Be careful, when you do the aggregate statistics, because it gets lost. Please rephrase." =¿ New comment: I understand that the data is aggregated, and the statement itself is not wrong. But it is a little "dangerous" to simply aggregate and conclude that the power production is larger for all multi-rotors wind farms(line 299-300) when the power difference is occassionally negative, i.e. sometimes it is better to have a single rotor. If you included more turbines the advantage migth disappear all together.

Response:
The text in lines 323 -325 has been reworded to reflect the reviewer's concern.
Finally, I also wish to point the authors attention to a newly published paper, which examines many of the same things and have comparable findings.

Response:
A reference to this paper is included in Section 1, lines 45 -49, in the revised manuscript.
We thank Dr. S. J. Andersen once again for his careful assessment and hope to have addressed all concerns.

Introduction
Wind energy is among the fastest growing sources of renewable energy worldwide. Understanding and mitigating the deleterious effects of interactions between wakes of multiple turbines is critical for efficient utilization of the wind resource. In large wind farms, the wake interactions can limit the power density, or the power extracted per unit land area. The turbulent wake interactions also determine fatigue loads on downstream turbines, which has a direct bearing on the levelized cost of energy. Previous work has shown that wake losses are closely tied to wind farm layout parameters such as inter-turbine spacing (Meyers and Meneveau, 2012;Yang et al., 2012), alignment between columns and the wind direction (Stevens et al., 2014a;Ghaisas and Archer, 2016), horizontal staggering between adjacent rows (Archer et al., 2013) and vertical staggering of similar or dissimilar turbines (Vasel-Be-Hagh and Archer, 2017;Xie et al., 2017;Zhang et al., 2019). 25 The idea of mounting multiple rotors per tower has been explored in recent years Branney, 2012, 2014;Chasapogiannis et al., 2014;Ghaisas et al., 2018;van der Laan et al., 2019;Bastankhah and Abkar, 2019). For example, Jamieson and Branney (2012) pointed out that the scaling laws for power and weight with the diameter of a turbine (the 'squarecube law') pose a challenge to upscaling the design of current single-rotor turbines to very large systems, but make multi-rotor turbines an attractive alternative. Structural considerations with designing a 20 MW multi-rotor system were investigated in 30 Jamieson and Branney (2014). Their results suggested that for a 45-rotor 20 MW system, the benefits due to reduced rotor and drive train costs would outweigh potential challenges associated with a more complicated tower structure. Chasapogiannis et al. (2014) studied the aerodynamics of a 7-rotor system, with the tips of the blades of adjacent rotors spaced 0.05 diameters apart. Interference due to adjacent rotors was found to lead to approximately 3 % increase in power, while about 2 % increase in the blade loading amplitude was observed.

35
Analysis of the wake of a 4-rotor turbine was carried out in our previous work  using large eddy simulation (LES). It was shown that the multi-rotor turbine wakes recover faster compared to wakes of an equivalent singlerotor turbine. The turbulent kinetic energy added due to multi-rotor turbines was also lesser than that due to an equivalent single-rotor turbine. Wind farms comprising of five aligned turbines spaced four diameters apart were also considered in this study. The potential for reduced wake losses as well as reduced fatigue loads was clearly pointed out. 40 The results for the wake of an isolated turbine were confirmed recently in van der Laan et al. (2019) using a combination of field observations and numerical simulations. van der Laan et al. (2019) also studied the aerodynamics of individual and combined rotors. It was found that rotor interaction can lead to an increase of up to 2 % in the power generation, similar to that reported in Chasapogiannis et al. (2014). Isolated multi-rotor turbines were studied in detail in van der Laan et al. (2019), and potential benefits in multi-rotor wind farms were discussed. Bastankhah and Abkar (2019) also studied isolated multi-45 rotor wind turbine wakes and found similar wake recovery characteristics. Multi-rotor configurations other than the 4-rotor configuration studied in the present paper and elsewhere were considered. The effect of number and direction of rotation of the individual rotors on the rate of wake recovery was also studied, and was found to be negligible by Bastankhah and Abkar (2019).
Interactions between several multi-rotor wind turbines arranged in a 4×4 grid were studied using several Reynolds-Averaged 50 Navier Stokes (RANS) simulations and one large eddy simulation (LES) in van der Laan and Abkar (2019). The annual energy production of multi-rotor wind farms was found to be 0.3 − 1.7% larger compared to that of equivalent single-rotor wind farms. The benefit was confined to the first downstream turbine row and for cases where the wind direction was fully aligned with the turbine columns. This discrepancy with the results of Ghaisas et al. (2018) can be attributed to the large tip spacings considered in Ghaisas et al. (2018). In the present work, we study more realistic tip spacings, and observe consistent qualitative 55 and quantitative trends with the results of van der Laan and Abkar (2019).
In this paper, we extend our previous work  by considering a larger number and range of multi-rotor wind turbine and farm design parameters. A schematic of the multi-rotor turbine considered here is shown in Fig. 1(b). Four rotors with identical diameters, d, are mounted on a tower with height H T (Fig. 1b). The tips of the rotors are separated by s h and s v in the horizontal and vertical, respectively. As a result, the rotors are centered at H T ± (s v + d)/2, and the mean hub-60 height is H T . The multi-rotor configuration (henceforth referred to as 4-rotor turbine) is compared to a conventional turbine with a single rotor (referred to as 1-rotor turbine) with diameter D = 2d per tower with height H T (Fig. 1a). The total frontal rotor area is πD 2 /4 in each case.
The primary aim of this paper is to quantify the benefits associated with the wakes of multi-rotor turbines for a wide range of tip spacings, thrust coefficients and inter-turbine spacings using LES. A second aim is to develop an analytical modeling This paper is organized as follows. The LES methodology, details of the simulations and the analytical framework are described in Sect. 2. Results of isolated 4-rotor turbines are described in Sect. 3, while results of wind farms comprised of 4-rotor turbines are described in Sect. 4. In each case, LES results are presented followed by predictions of the analytical modeling framework. Sect. 5 presents a brief summary and the conclusions. The LES-filtered incompressible Navier-Stokes equations are solved on a structured uniform Cartesian mesh using Fouriercollocation in x and y directions, sixth-order staggered compact finite-differences in the z direction and a total variation diminishing (TVD) fourth-order Runge-Kutta time-stepping scheme. Non-periodicity is imposed in the x direction using a fringe region technique (Nordström et al., 1999). Partial dealiasing is achieved by applying the 2/3 rule in x, y and the use of skew-symmetric form for the convective terms in the z direction. The governing equations and numerical discretization 85 details may be found in Ghate and Lele (2017) (Appendix A). The effect of sub-filter scales is modeled using the Anisotropic Minimum Dissipation (AMD) model (Rozema et al., 2015). Wind turbine forces are modeled as momentum sinks using the actuator drag-disk model (Calaf et al., 2010). The turbine forces in the LES are defined in terms of the disk-averaged velocity and a 'local thrust coefficient', C T . The local thrust coefficient (assuming validity of the inviscid actuator-disk theory) is related to the nominal thrust coefficient, C T , through the relation C T = 16C T / (C T + 4) 2 , or equivalently, through the relations where a is the axial induction factor. Algebraic wall models based on the Monin-Obukhov similarity theory are used to specify the shear stresses at the bottom wall. Viscous stresses in the rest of the domain are smaller than the sub-filter scale stresses by around 8-10 orders of magnitude and, hence, are neglected in these simulations. The code has been validated over several previously published studies Ghaisas et al., 2017;Ghate et al., 2018).

Cases Simulated
Half-channel (HC) simulations are carried out using the concurrent precursor-simulation methodology (Stevens et al., 2014b) on domains of length L x , L y , L z in the three coordinate directions. A schematic of the simulation domain is shown in Fig. 1(c).
All simulations use (L y , L z ) = (π/2, 1) H, while L x = πH or 1.25πH, depending on the case. Here H is the height of the half-channel. The flow in the 'precursor' simulation is driven by a constant imposed pressure gradient, −u 2 * /H, where u * is the 100 friction velocity at the bottom wall. The HC configuration is used as a model for the neutrally-stratified atmospheric boundary layer (ABL) with the Coriolis forces neglected (Stevens et al., 2014a;Calaf et al., 2010), and we use the terms HC and ABL interchangeably. The surface roughness height at the bottom wall is z 0 = 10 −4 H. This corresponds to rough land, and has been used in previous wind turbine studies (Calaf et al., 2010  in one simulation. All isolated turbines, and the most upstream turbine in the five-turbine cases, are located at x = 0, where the domain inlet is at x = −4D. The turbine towers are located at y = Ly/2 in the spanwise direction and the tower height is H T = 0.1H for all turbines. The domain size in the x-direction is increased to 1.25π to accommodate larger axial spacings for the cases with S X = 5D or 6D. Field measurements and simulations reported in van der Laan et al. (2019) show that the bottom pair of rotors has a slightly 125 larger thrust coefficient than the top pair of rotors. However, for simplicity, the same thrust coefficient is used for all rotors in one simulation. The methodology of keeping thrust forces identical across all rotors of the multi-rotor turbine was adopted by van der Laan et al. (2019) as well in the part of their study that focused on comparing wakes of multi-rotor and single-rotor turbines. The effect of variable operating conditions for the top and bottom pairs of rotors can be studied systematically in the future. Finally, the appropriateness of considering a single-rotor turbine with the same total frontal area, thrust coefficient and 130 mean hub height as that of the multi-rotor turbine is evaluated in Appendix B.

Analytical Model
An analytical modeling framework based on the model by Bastankhah and Porté-Agel (2014) is evaluated for the multi-rotor configuration in this paper. The model assumes that the velocity deficit in the wake decays in the streamwise (x) direction, and follows a Gaussian profile in the radial directions. The deficit due to turbine rotor i located at (x i , y i , z i ) at a downstream point 135 (x, y, z) is given as for x > x i . The length scale d 0 equals D for 1-rotor and d for 4-rotor cases. The argument of the square-root in eq. (2) is set to zero whenever it is less than zero, which happens very close to the turbines.

140
The combined effect of multiple turbine rotors has been modeled in the past using several empirical techniques. Primary among these are addition of velocity deficits (implying linear addition of momentum deficit), square-root of sum of squares (implying addition of kinetic energy deficit; also termed as quadratic merging), and considering the largest deficit to be dominant. In this study, a hybrid between the first two approaches is found to give best results. Appendix A presents brief comments justifying the hybrid approach. The hybrid approach involves linear merging of wakes originating at the same x location, with 145 quadratic merging of wakes originating at different x locations. This can be written as N xt is the number of unique axial locations where a turbine is located. N r (x t ) is the number of rotors at the location x t . In this paper, N xt is 1 for the isolated turbine cases, and 5 for the wind farm cases. Furthermore, since we only consider either an isolated turbine or a wind farm with one column of turbines, N r is 1 for the 1-rotor cases and 4 for the 4-rotor cases. Finally, the mean velocity at each point in the domain is calculated according tō u(x, y, z) =ū up (z) − ∆ū tot (x, y, z). (4) The upstream velocity is assumed to follow the logarithmic profile,ū up (z) = (u * /κ)ln(z/z 0 ), with κ = 0.4.
This modeling framework involves two empirical parameters, k * and σ 0 . Comments regarding selecting these parameters are provided in the appropriate sections below.

Grid Convergence and Baseline Cases
Precursor ABL simulation results are shown first in Fig. 2. These results are averaged over time and the horizontal directions.
As expected, the mean streamwise velocity profiles follow the logarithmic law of the wall, particularly in the lower 20% of the domain. The total shear stress profiles also follow the expected line with slope equal to -1. This indicates that the vertical 160 transport of momentum by the ABL is correctly represented by the numerical method and AMD subgrid-scale model, and that the ABL simulations are statistically stationary. Figure 2 also shows that the spatial resolution employed is adequate for these ABL simulations, since the results are almost independent of the grid size.
Results of an isolated 1-rotor turbine and an isolated 4-rotor turbine with s/d = 0.05 are shown in Fig for the 1-rotor configuration. The 4-rotor configuration has two mid-span planes, Y cen = L Y /2 ± (1 + s h )d/2. Results at only one of these, at L Y /2 − (1 + s h )d/2, are shown, since both planes are statistically identical.
Figure 3(a) shows that the velocity deficit profiles for the 1-rotor turbine have a single peak close to z/H = 0.1. Two distinct peaks, close to z/H = 0.1 ± (1 + s v )d/2, are seen for the 4-rotor turbine wake in Fig. 3 Simulations with varying grid sizes (the IT*-s cases) show that the differences between the results reduce as the grid is refined. In general, the sensitivity to grid resolution is larger for the 4-rotor case as compared to for the 1-rotor case. This is expected because the 4-rotor configuration involves smaller length scales, associated with the smaller diameter of the individual 175 rotors, and the tip-spacing. The differences between the velocity deficits obtained using grids G3, G4 and G5 are not easily discernible on the scale of Figures 3(a,b). Differences between the results of grid G2 and those of finer grids are easily apparent only at x/D = 2 for the multi-rotor configuration. The double-peaked shape of the velocity deficit at this location is not fully resolved using grid G2, and is better resolved using grids G3 and finer. The velocity deficit values, averaged over the rotor disk regions, for different grid sizes are used to assess grid convergence. Taking the results of grid G5 as reference, the errors in 180 velocity deficits obtained using grid G2 are 3.2% and 1.9% at x/D = 4 and x/D = 6, respectively.
The added TKE profiles in Figures 3(c,d) show greater sensitivity to grid size than the mean velocity deficits. The resolved portion of the TKE is expected to increase with increasing grid resolution. It should be noted that the resolved TKE cannot be supplemented with a subgrid contribution in an LES using an eddy-viscosity model, where only the deviatoric part of the stress is modeled. Except for a small region close to z/H = 0.15 at x/D = 6, over most of the domain, the resolved portion of the 185 added TKE is also found to increase with increasing resolution. The turbulence intensity averaged over the rotor area is found to change by around 15% at x/D = 4 and 6 between grids G2 and G5. Between grids G2 and G3, the disk-averaged turbulence intensity values vary by 6.5% at x/D = 4 and by 3.5% at x/D = 6.
A change of 3.2% in the disk-averaged velocity deficit on doubling the grid resolution (from G2 to G5) implies a change of approximately 9.9% in the averaged power. The results pertaining to estimates of power, in particular the comparisons 190 between LES and analytical model predictions, presented in this manuscript should be interpreted keeping this limitation in mind. The computational costs per simulation were approximately 4400 CPU-hours and 70000 CPU-hours on grids G2 and G5, respectively. Even with near-perfect scaling, as was obtained with very careful attention to parallel implementation in our code, in view of the large parameter space to be evaluated, it was decided to conduct all further simulations on grid G2. For the wind farm cases with domain size increased to 1.25π in the x direction, the number of points in the x direction is increased 195 to 320 to retain the same resolution. This grid is labeled as G2A in Table 1. The grids G2/G2A imply that the smaller rotor disk (diameter d) is resolved by 4 × 8 points and the composite wake of the multi-rotor turbine (diameter D) is resolved by 8 × 16 points in the y − z plane. The details in the region between the rotor tips are obviously missed. However, as shown in the next subsection, the overall effect of varying tip spacing is captured, because the actuator disk model appropriately adjusts the distribution of forces across the discretization points. It should be noted that the level of resolution of the composite wake 200 is in keeping with the recommendation in Wu and Porté-Agel (2011), and is comparable to the grid resolution used by several previous studies (Calaf et al., 2010;Stevens et al., 2014a).

Effect of Tip Spacing
Isolated 4-rotor turbines with varying tip spacings, s h = s v = s, are studied in this subsection (IT2-s cases). Contours of the mean streamwise velocity deficit and the TKE (Fig. 4) in the mid-span planes show that one large wake immediately 205 downstream of a 1-rotor turbine is replaced by four smaller wakes immediately downstream of the four rotors of the 4-rotor turbines. Comparing Figures 4(a,c,e), it is clear that the wake of a 4-rotor turbine at any downstream location (e.g. at x/D = 4), is weaker in magnitude than that of the 1-rotor turbine. This is also seen in the profiles shown in Fig. 5. In other words, the wake of a 4-rotor turbine is seen to recover faster than the wake of a 1-rotor turbine with the same thrust coefficient and rotor area. Figure 5 also shows that greater the tip spacing of the 4-rotor turbine, faster is the wake recovery. This is also indicated 210 by the shortening of the contour lines corresponding to ∆u/u * = 1 and 2.5 in Fig. 4 with increasing tip spacing.
An intuitive explanation for the increasing rate of wake recovery with increasing tip spacing is as follows. The characteristic length scale of the wake of the 1-rotor turbine is diameter D, while that for the individual wakes of the 4-rotor turbines is the smaller diameter d. Furthermore, the spacing between the tips of the 4-rotor turbine allows for greater entrainment of low-momentum fluid into the 4-rotor turbine wakes. As a result, the rate of wake recovery is larger for the 4-rotor turbine as 215 compared to the 1-rotor turbine, and increases with increasing s.
The wakes of the individual rotors of a 4-rotor turbine expand with downstream distance, and eventually merge to form a single wake. The axial distance where individual wakes of the four rotors may be considered to have merged increases with increasing s. This is seen clearly in Fig. 5, where two peaks in the velocity deficit profiles are not seen at x/D = 4 for the s/d = 0.1 turbine, while two peaks are clearly visible at x/D = 6 for the s/d = 0.5 turbine.

220
The contour plot of TKE shown in Fig. 4(b) is strikingly similar to those reported previously (e.g. Fig. 18 in Abkar and Porté-Agel (2015)) for an isolated 1-rotor turbine. The TKE contours in Fig. 4(b) are similar in shape to those in Fig. 4(d) beyond approximately x/D = 4, but are quite dissimilar to the contours in Fig. 4(f). This is further evidence for the observation that the wake-merging distance increases with increasing s. The rotors of the 4-rotor turbine behave independently up to increasingly larger downstream distances with increasing s.

225
A succinct representation of the effect of tip spacing on the wake of an isolated 4-rotor turbine with respect to that of an isolated 1-rotor turbine is shown in Fig. 6, where rotor-disk-averages of four quantities is plotted as a function of the axial distance . The rotor-disk averages are calculated at each axial (x/D) location and over different regions in the y − z plane depending on the turbine configuration. The averages are computed over one disk of diameter D, centered at (L Y /2, 0.1H) for the 1-rotor turbine, and over four disks of diameters d each, centered at (L Y /2 ± (1 + s h )d/2, 0.1H ± (1 + s v )d/2), for the 230 4-rotor turbines. The disk averaged TI is actually the ratio of the square-root of the disk-averaged TKE and disk-averaged mean streamwise velocity, being slightly different from the disk-average of the point-wise turbulence intensity. The disk-averaged added turbulence intensity is defined as,  The disk-averaged added turbulence intensity can be compared to that reported in Figures 18(b,d,f)  x/D=4 x/D=6 1-Rot s=0.0d s=0.1d s=0.2d s=0.5d

Effect of Thrust Coefficient
The IT2-C T cases, along with two cases from the IT2-s set of simulations, are compared to study the effect of thrust coefficient.

250
Only one 4-rotor configuration, with tip spacing s/d = 0.1, is considered here. Figure 7 shows that the trends observed for C T = 4/3 hold for the other two thrust coefficients studied as well. The disk-averaged velocity deficits are smaller for the 4-rotor turbine than for the corresponding 1-rotor turbine. The added TKE (not shown) and T I disk are also smaller for the 4-rotor turbine than for the 1-rotor turbine for all the thrust coefficients studied.

255
The analytical modeling framework predicts the mean velocity deficits of the 1-rotor and 4-rotor turbines accurately. Empirical parameters values k * = 0.025 and σ 0 /d 0 = 0.28 were found to lead to accurate predictions for all the cases investigated. Here, In particular, Fig. 5 shows that the radial profiles of the velocity deficit at several downstream locations, and for turbines with 260 different tip spacings, are predicted quite accurately. Slight under-predictions or over-predictions are observed very close to the turbine, but the overall predictions are accurate, particularly beyond x/D = 2. Disk-averaged velocity deficit profiles are also predicted accurately, but are not shown on Fig. 6(a) to avoid clutter. Figures 7(a-b) show that the Gaussian analytical model is reasonably accurate at predicting the disk-averaged velocity deficit for all thrust coefficients beyond the very-near-wake region, i.e. approximately beyond x/D = 2.

4 Multi-Turbine Simulation Results
Wind farms comprised of a line of five turbines aligned with each other and with the mean wind direction are studied here.
These cases are labeled WF* in Table 1.

Effect of Tip Spacing
The effect of tip spacing on the contours of velocity deficit and TKE are seen in Fig. 8. The axial spacing between different 270 turbines in the wind farm is kept fixed at 4D and the thrust coefficient is 4/3 for all rotors of all turbines. It is clear that the velocity deficits are significantly different between the 1-rotor and 4-rotor wind farms, as well as between 4-rotor wind farms with different tip spacings. The single wake behind the turbines in the 1-rotor wind farm are replaced by four smaller wakes behind the turbines in the 4-rotor wind farms. The wakes move further apart in the radial directions as the tip spacing increases.
Similar to the TKE distribution behind an isolated 1-rotor turbine, the TKE values are largest around the top-tip height of the 275 turbines.
The effect of tip spacing on 4-rotor wind farms is quantified in Fig. 9. Focusing on Fig. 9(a-b), the profiles of the velocity deficits averaged over the rotor disk and T I disk have local maxima close to the turbine locations, i.e. at x/D = 0, 4, 8, 12 and 16. The velocity deficit profile for the 1-rotor wind farm has a maximum close to turbine 2 (located at x/D = 4), as seen in have approximately equal magnitudes. The T I disk profiles in Fig. 9(b) show similar behavior for the 1-rotor wind farm.
The velocity deficits of the 4-rotor turbines are seen in Fig. 9(a) to be smaller than those of the 1-rotor turbine for the first two turbines (x/D = 0, 4). In this region, x/D < 8, the deficits decrease with increasing tip spacing, which is consistent with the observations for isolated turbines (Fig. 6(a)). The deficits accumulate and the disk-averaged profiles for all 4-rotor wind farms are almost equal to that for the 1-rotor wind farm for turbine rows 3 onward (for x/D > 8). The turbulent intensity profiles are 285 smaller for the 4-rotor wind farms than for the 1-rotor wind farm, and decrease with increasing s/d. This sensitivity to the tip spacing persists downstream of all turbines, unlike the velocity deficits, which are sensitive only downstream of the first two turbines.
The relative powers of the turbines are shown in Fig. 9(c). The power of the first (or front) turbine is used for normalization in each wind farm. Thus, the relative power for turbine i is calculated as P i /P 1 = u 3 i /u 3 1 , where the overhead bar represents time-290 averaging and subscript i denotes the location of the turbine within the wind farm. The relative power of turbine 2 (x/D = 4) in the 1-rotor wind farm is minimum, and the relative power profile shows a slight recovery for turbines 3-5. This is consistent with the maximum for the velocity deficit at turbine 2, seen in Fig. 9(a). The relative powers of turbines in the 4-rotor wind farms are sensitive to the tip spacing as well as the turbine location. For s/d = 0.1, only turbine 2 has larger relative power than turbine 2 of the 1-rotor wind farm, while for s/d = 0.5, turbines 2-4 have larger relative powers than the corresponding 295 turbines of the 1-rotor wind farm. All these trends are consistent with the velocity deficit profiles seen in Fig. 9  Interaction between the effects of tip spacing and axial spacing are also seen on comparing Figures 9(c) and (f). For instance, the relative powers of turbines 2 and 3 of the wind farm with s/d = 0.5 are appreciably larger than the corresponding turbines of the 1-rotor wind farm, when the axial spacing is 4D. However, relative power of only turbine 2 of the wind farm with tip spacing s/d = 0.5 is appreciably larger than that of the corresponding 1-rotor wind turbine, when the axial spacing is increased 310 to 6D. Thus, tip spacing has a greater effect on the relative power in a closely spaced wind farm. Figure 10 shows that the trends observed for C T = 4/3 hold for other values of thrust coefficient as well. The velocity deficit and turbulence intensity are larger for cases with larger thrust coefficient. For each value of C T , the velocity deficit of the 4-rotor wind farm is generally smaller than that of the 1-rotor wind farm downstream of the first two turbines (for approximately x/D < 8) and are almost equal beyon this. Since the tip spacing of the 4-rotor wind farm is s/d = 0.1, only turbine 2 shows 315 a larger relative power in the 4-rotor wind farm compared to the 1-rotor wind farm, consistent with the observation made in Figure 9. For C T = 2, the velocity deficit profiles cross over, and the 4-rotor profile is larger than the 1-rotor profile, in a small region upstream of turbine 3. As a result, the relative power of turbine 3 is smaller in the 4-rotor wind farm compared to the 1-rotor wind farm. However, this crossover in power is smaller in magnitude than the values for turbine 2, such that the collective relative power of the downstream turbines is larger for the 4-rotor wind farm than for the 1-rotor wind farm.

320
The effect of all governing parameters (s, S X , C T ) on the wake losses in multi-rotor wind farms is presented succinctly in Fig. 11. Figure 11(a) shows the average power of turbines 2 through 5 (P 2−5 = (1/4) 5 i=2 P i ), normalized by the power of the front turbine in each wind farm. Aggregation of relative powers across all downstream rows, as done here, can hide negative power differences (associated with the crossovers referred to above) that might occur at individual turbine rows. Despite this, the aggregated relative power is a useful measure of the overall wake losses associated with a particular wind farm. It is seen 325 that P 2−5 /P 1 is larger for all 4-rotor wind farms than the corresponding 1-rotor wind farm with the same thrust coefficient and axial spacing. The benefit increases with increasing tip spacing.
Each data point in Fig. 11(a) is normalized by the power of the front turbine in the respective wind farm. The front turbine power is expected to be similar to that of an isolated turbine, and hence, is expected to be dependent on the thrust coefficient, but not on the axial spacing. This is seen to be the case in Fig. 11(b), where the power of the front turbine extracted from the 330 different wind farm cases are shown. For comparison across cases with different thrust coefficients, all powers are normalized by the power of the front turbine in the 1-rotor wind farm with the same thrust coefficient. The front turbine powers are independent of the axial spacing, and lines corresponding to S X = 5D and 6D lie on top of the line corresponding to S X = 4D. Figure 11(b) also shows that the front turbine power in 4-rotor wind farms is weakly dependent on the tip spacing. As the tip spacing varies over s/d = 0.1 to 0.5, the front turbine power varies by 3.5%, 2.7% and 3.2%, with the thrust coefficients fixed 335 at 1, 4/3 and 2, respectively. We note that this variation cannot be explained by the variation in power potential due to different tip spacings (see Appendix B), and is likely caused by the effects of turbulent mixing in the wake (Nishino and Wilden, 2012), which are different for different tip spacings.
To account for the differences in the front turbine power, the average power of turbines 2 through 5 is replotted in Fig. 11(c), with only the 1-rotor front turbine powers used for normalization. The same qualitative conclusions can be drawn from 340 Fig. 11(c), as were drawn from Fig. 11(a), although the magnitudes of the benefit are larger. Finally, the differences between the relative powers of the 4-rotor and 1-rotor configurations are plotted in Fig. 11(d). This plot is directly derived from Fig. 11(c) by subtracting the data points corresponding to the 1-rotor wind farm from the 4-rotor wind farm data, i.e.
. This quantity measures the extent by which wake losses in a 4-rotor wind farm are smaller than wake losses in a 1-rotor wind farm with the same inter-turbine spacing and with all rotors operating with the same thrust 345 coefficient. The benefit of 4-rotor wind farms increases with increasing tip spacing and with decreasing thrust coefficient. The effect of axial spacing on the benefit is slightly ambiguous. For a fixed thrust coefficient and tip spacing, the benefits are largest for S X = 4D, and are almost equal for S X = 5D and 6D. Appendix C shows that the conclusions drawn above are not affected by the fact that the first turbine powers are significantly different between the 1-rotor and 4-rotor wind farms.

Analytical Model
Predictions of the analytical modeling framework for wind farms comprised of a line of five turbines are examined in this section. The parameter k * , which controls the growth rate of the wake, is extracted from all the 1-rotor wind farm LES. First, the wake widths in the y and z directions are calculated using the definition outlined in Bastankhah and Porté-Agel (2016).
where (Y cen , Z cen ) = (L Y /2, 0.1H) are the mid-span and mid-vertical planes of the 1-rotor wind turbine wakes, and ∆ū max (x) is the maximum of the velocity deficit at location x. The wake width is then calculated as the geometric mean of the wake widths in the two transverse directions, σ = √ σ y σ z . Wake widths extracted from three 1-rotor LES with fixed S X = 4D and varying thrust coefficient are shown in Fig. 12(a).
Turbines are located at x/D = 0, 4, 8, 12 and 16 in this plot. Moving downstream from one turbine location, the wake widths generally increase, until the effect of the next downstream turbine is felt. The wake width profiles show dips close to the turbine locations, followed by regions of growth. Regions where the wake widths grow approximately linearly are identified with black solid lines in Fig. 12(a). These black solid lines are linear fits to the data, and the extents of the linear fitting region are identified visually. The slopes of these lines yield the wake growth rate parameter, k * .

365
The wake growth rate parameter values for all turbines in the 1-rotor wind farm simulations are compiled in Fig. 12(b).
The k * values are plotted against the streamwise turbulence intensity, I x , at each turbine rotor disk. As observed in previous studies, the wake growth rate increases with increasing turbulence intensity. The solid blue line fits the data with a correlation coefficient of 0.8. In subsequent model runs for 1-rotor and 4-rotor wind farms, this linear regression model is used to determine k * , with I x extracted from the LES results.

370
Model predictions are compared to LES results for two cases in Fig. 13. The sensitivity of the model predictions to the second tunable parameter, the initial wake width σ 0 , is seen in this figure. Figure 13(a) shows that the disk-averaged velocity deficit is over-predicted by the analytical model with σ 0 /D = very close to the turbines, while it is under-predicted (to a lesser degree) with σ 0 /D = 0.32. Farther away from the turbines, approximately between 1D to 3D downstream of each turbine, using σ 0 /D = 0.28 yields good agreement with the LES results, while using σ 0 /D = 0.32 continues to yield under-predictions.

375
The power predictions shown in Fig. 13(b) also show sensitivity to the value of σ 0 . The relative power of turbine 2 is captured accurately with σ 0 /D = 0.28, while the relative powers of further downstream turbines are under-predicted by around 10%.
With σ 0 /D = 0.32, the relative power of turbine 2 is over-predicted, while that of further downstream turbines is in better agreement with the LES results. Similar conclusions can be drawn from the results of the 4-rotor turbine with s/d = 0.1, shown in Figures 13(c) and (d). In summary, σ 0 /D = 0.28 leads to better prediction of the mean velocity deficit in the wake 380 region (1D−3D downstream), while σ 0 /D = 0.32 leads to better prediction at the turbine locations, as evidenced by the better predictions of the power. Thus, the combination of model parameters which leads to accurate predictions in the wake does not necessarily lead to accurate predictions of power, for which, the values at and very close to the turbines need to be predicted accurately.
The influence of using spatially constant values for the wake growth rate parameter on the model predictions is shown in It is seen that using a spatially non-varying k * leads to a gradual decrease in the relative power with turbine number. The LES results show the characteristic feature of recovery of the relative power after turbine 2 in the 1-rotor wind farm and after turbine 3 in the 4-rotor wind farm. This feature is not captured for any combination of σ 0 /D and non-varying k * . Comparing 390 Figures 14(a,b) and Figures 13(b,d) respectively, it is clear that the power degradation recovery is better captured using k * that varies spatially depending on the local turbulence intensity. Similar observations were reported previously for 1-rotor wind farms (Niayifar and Porté-Agel, 2016), and are seen here to hold for several 4-rotor wind farms as well. It is possible for some cases, particularly the s/d = 0.5 wind farms, where the relative power continues to gradually decrease until the fifth turbine (see Figures 15 and 16), to be better predicted using a spatially constant k * value. However, no single combination of spatiallyconstant k * and σ 0 /D values was found that resulted in good predictions for all cases. In view of the cases investigated here, we prefer the use of a spatially-varying k * dependent on the local turbulence intensity, consistent with previous studies for 1-rotor wind farms (Niayifar and Porté-Agel, 2016).
Relative power predictions for all the wind farm cases are compared to LES results in Figures 15 and 16. The average error in predicting the relative powers of turbines 2 through 5 are shown in each case. The k * values are obtained as outlined above, 400 while σ 0 /d 0 = 0.28 is used for all cases, where d 0 equals D for the 1-rotor cases and equals d for the 4-rotor cases. The absolute errors in relative power averaged over turbines 2 through 5 ((1/4) 5 i=2 | (P i /P 1 ) LES − (P i /P 1 ) model |) are shown in Figures 15 and 16. It should be noted that this level of accuracy is similar to that observed in previous studies (Stevens et al., 2015(Stevens et al., , 2016 of wind farms that are finite in axial as well as spanwise directions, and where the wind is directed along only one direction, or averaged over a very narrow (less than 2 • ) sector.

405
The errors are seen to be smallest for the 1-rotor cases. For 1-rotor wind farms, typically, the power of the second turbine is smallest, and there is a slight recovery for turbines 3, 4 and 5. This behavior is reproduced well by the analytical model.
In the 4-rotor cases, the relative power saturates farther into the wind farm, typically at the third row for s/d = 0.1 and 0.25.
For s/d = 0.5, the power continues to decrease until the fifth row for most cases. The model predictions, on the other hand, typically saturate by the second row. Thus, the errors are largest for the second row, although the relative power level of turbines 410 in the fourth and fifth rows is typically well captured.
In conclusion, the analytical modeling framework is capable of reproducing LES results of 1-rotor and 4-rotor wind farms with reasonable accuracy, comparable to previous results for 1-rotor turbines (Stevens et al., 2015). Improved prediction of the region very close to the turbine is needed to further improve the accuracy of the model at predicting the power degradation and wake losses in wind farms.

Discussion and Summary
This paper is devoted to studying the turbulent wake of a multi-rotor wind turbine configuration, and to comparing it with a conventional single-rotor wind turbine wake. The potential benefits offered by this configuration, with four rotors (with diameters d = D/2) mounted on a single tower, over the conventional single-rotor turbine (with diameter D) are studied in detail. Large eddy simulation is used as the primary tool for this work, Applicability of an analytical modeling framework 420 based on the assumption of Gaussian radial profiles of velocity deficits to the multi-rotor configuration is also examined.
The LES results outlined in Sect. 3 show that an isolated 4-rotor turbine wake recovers faster compared to an isolated 1-rotor turbine wake. The isolated 4-rotor turbine wake also shows smaller TKE levels in the rotor disk region. A simple physical reason for this faster wake recovery and lower TKE levels is that the greater perimeter-to-area ratio of the multi-rotor turbine allows for greater entrainment of low momentum fluid into the wake. The behavior of the wake is sensitive to the tip spacing 425 (s/d), with faster wake recovery seen for for larger s/d. This is consistent with the simple physical reasoning presented above, since if s/d is very large, each rotor of the multi-rotor turbine behaves independently of other rotors, and the wake of each in accurately reproducing these trends are partly due to the fact that the Gaussian wake model is valid only beyond a certain distance downstream of a turbine, and is not valid immediately upstream and immediately downstream of a turbine. Thus, this study points to the need for better analytical modeling of the region very close (upstream as well as downstream) to the turbine.
The actuator drag-disk model provides a crude representation of the processes occurring very near the turbine disks. While this crude representation is sufficient for the purposes of capturing the interactions between the turbines and the atmospheric 455 boundary layer, future studies should focus on using the actuator disk/line models with rotation of the blades included. Potential benefits associated with co-rotation and counter-rotation of the rotors in the multi-rotor configuration can be studied. Recent work by Andersen and Ramos-Garcia (2019) suggests that interaction between tip vortices of the individual rotors of the multirotor turbine aids in breakdown and recovery of the wake. These beneficial interactions might be missing from multi-rotor turbines with very large tip spacings, thus slowing down the rate of wake recovery. This issue can also be studied in the future. should also be studied in the future. Finally, developing better analytical models for both, 1-rotor and multi-rotor, configurations continues to be a persistent challenge in wind energy research, and will be pursued in future work.

Appendix A: Hybrid Linear-Quadratic Wake Superposition Methodology
A brief justification for following the hybrid linear-quadratic methodology of wake merging is provided in this appendix. (∆ū j (x, y, z)) p   1/p , with p = 1 and 2 corresponding to linear and quadratic merging, respectively. It is clear that linear merging gives better agree-470 ment with LES results compared to quadratic merging. Thus, for wakes originating at the same x location (i.e. 'adjacent' wakes), linear merging is preferred. Figure A1(b) compares LES results and model predictions for the s/d−0.1, C T = 4/3 and S X = 4D wind farm. Here, linear superposition of adjacent wakes is assumed, and superposition of these combined wakes originating at different x locations is examined.The choices evaluated here are with, once again, p = 1 and 2 corresponding to linear and quadratic merging. For this case, N xt = 5 and N r = 4 for all x t . Figure A1(b) shows that linear merging (p = 1) leads to a continuous increase of the velocity deficits, which is unphysical.
Quadratic merging leads to velocity deficits that saturate a few turbines into the wind farm, and is in better qualitative and quantitative agreement with the LES results. Thus, quadratic merging is preferred for wakes originating at different x locations.

480
Thus, a hybrid linear-quadratic merging strategy is seen to give best results. It should be noted that this is an empirical choice, and a physics-based/first-principles approach for wake superposition is a topic of active research.

Appendix B: Potential Power of Multi-Rotor Wind Turbines
Finding an appropriate single-rotor turbine, which can be considered as a reference against which a multi-rotor turbine can be compared, is not straightforward. This is because the lower and upper pair of rotors in the 4-rotor configuration are subjected 485 to different wind speeds and turbulence levels as compared to each other and to the single rotor in the 1-rotor configuration.
In this work, we consider a single-rotor turbine with the same total frontal area, same thrust coefficient and same mean hub height as a multi-rotor turbine to be a reference. To test the appropriateness of this assumption, the potential power, computed as P pot = πD 2 /8 C P U 3 0,disk , is shown in Table B1. Here, U 0,disk is obtained by averaging the logarithmic inflow profile (shown in Fig. 2a) over the rotor disks. The potential power normalized by that of the 1-rotor turbine, P pot /P 1−Rot pot , is also 490 shown in Table B1. A representative value of C P = 0.5625 is used, but this precise number does not matter when we compare the normalized potential powers. The normalized potential powers are seen to be almost equal to 1 for all the tip spacings, and slightly reduce as the tip spacing increases. This indicates that the net effect of shear and the chosen dimensions of the turbines is such that the effect of the reduced wind speed seen by the lower two rotors dominates the effect of th larger wind speed seen by the upper two rotors. This effect is not very strong, being only 2.4% for s/d = 0.5. For s/d = 1, the effect is larger, at 5.5%.

495
The same conclusion is reached if we use the hub height velocities instead of the disk-averaged velocities in computing P pot .
For the present study, the chosen 1-rotor configuration may be considered to be appropriate as a reference, since its potential power varies by less than 2.4% for the majority of the multi-rotor configurations.  pairs. C T = 1.14, 1.61 and 2.47 for the runs labeled 1R − C T , corresponding to C T = 1, 4/3 and 2 respectively.

Appendix C: C T -Matched 1-Rotor Wind Farms
Single-rotor and multi-rotor turbines with the same rotor area, same mean hub height and same thrust coefficient have been 500 considered to be equivalent in the main body of this paper. This equivalence was based on the 'local' thrust coefficient, C T .
Assuming validity of the inviscid actuator-disk theory, imposing a local thrust coefficient implies imposing an induction factor, a, and a thrust coefficient, C T . These quantities are related by The classical actuator-disk theory, however, is not valid for the turbine disks subjected to the sheared, turbulent boundary layer 505 inflow in this study. Consequently, given a value of C T , the implied values for a and C T are different from those predicted by eq. (C1). Furthermore, since the single rotor in a 1-rotor turbine and the four individual rotors in a 4-rotor turbine are subjected to different values of shear and turbulence intensity, the implied values of a and C T are different for the 1-rotor and 4-rotor turbines. As seen in Fig. 11(b), the power of the front turbine in 1-rotor and 4-rotor wind farms is different although identical C T values are used for all rotors.

510
In this appendix, three additional 1-rotor wind farm simulations are reported with S X = 4D and with C T adjusted such that the resulting C T is closer to those of the corresponding 4-rotor turbines. Through a trial-and-error approach, C T = 1.14, 1.61 and 2.47 were found to yield C T values that are within 1.5% of those of the 4-rotor wind farms with C T = 1, 4/3 and 2, respectively. These simulations are denoted as 'C T -matched' runs, and are labeled as 1R − C T in Figures C1 and C2 here. Figure C1 is a reproduction of Fig. 9(a-c) appended with the additional 1-rotor wind farm simulation with C T = 1.61. The 515 disk-averaged velocity deficit and turbulence intensity profiles are larger than for the 1-rotor wind farm, particularly at x/D = 4 (turbine 2). The resulting power degradation (Fig. C1c) is more severe at turbine 2, and almost identical to the 1-rotor wind farm for further downstream turbines. Figure C2 is a reproduction of Fig. 11 appended with results from all three 'C T -matched' runs. Focusing on the black line with squares in Fig. C2(b), it is seen that the power of the front turbine in the additional 1-rotor wind farm simulation (labeled

530
In summary, this appendix ensures that the qualitative conclusions regarding the benefits of the 4-rotor wind farms remain unchanged, regardless of whether '1-Rot' (C T -matched) or '1R-C T ' (C T -matched) 1-rotor wind farms are used for reference.
Code and data availability. The LES code used for these simulations is available on GitHub at https://github.com/FPAL-Stanford-University/PadeOps.
Data can be made available upon request from the corresponding author.