Dynamics of Large Scale Turbulence in Finite-Sized Wind Farm Canopy Using Proper Orthogonal Decomposition and a Novel Fourier-POD Framework

: Large scale coherent structures in the atmospheric boundary layer (ABL) are known to contribute to the power generation in wind farms. In order to understand the dynamics of large scale structures, we perform proper orthogonal decomposition (POD) analysis of a ﬁnite sized wind turbine array canopy in the current paper. The POD analysis sheds light on the dynamics of large scale coherent modes as well as on the scaling of the eigenspectra in the heterogeneous wind farm. We also propose adapting a novel Fourier-POD (FPOD) modal decomposition which performs POD analysis of spanwise Fourier-transformed velocity. The FPOD methodology helps us in decoupling the length scales in the spanwise and streamwise direction when studying the 3D energetic coherent modes. Additionally, the FPOD eigenspectra also provide deeper insights for understanding the scaling trends of the three-dimensional POD eigenspectra and its convergence, which is inherently tied to turbulent dynamics. Understanding the behaviour of large scale structures in wind farm ﬂows would not only help better assess reduced order models (ROM) for forecasting the ﬂow and power generation but would also play a vital role in improving the decision making abilities in wind farm optimization algorithms in future. Additionally, this study also provides guidance for better understanding of the POD analysis in the turbulence and wind farm community. visualization T.C.; and


Introduction
Wind farms in the atmospheric boundary layer (ABL) pose a complex dynamical system with turbulent phenomenology occurring at multiple length scales, mainly due to the interaction of atmospheric boundary layer turbulent eddies [1] and the turbulence generated by the wind turbine wakes [2][3][4]. These interactions are manifested not only in the small scale structures (of the order or smaller than the turbine rotor diameter) but also in the large scale structures (one or two orders of magnitude larger than the turbine rotor diameter) [5]. In the current paper, we are concerned with the understanding of these large turbulent structures which arguably can serve as mediators to the generation of large scale motions or very large scale motions [6] (LSMs or VLSMs) which has significant contribution to wind farm power [4,7]. Several researchers have shown before that large scale counter-rotating roll cell structures are formed in the atmospheric boundary layer [8,9]. From the perspective of understanding the phenomenology, Chauhan et al. [10] predicted the dynamics of such large structures using two-point correlations from field experiments in Utah's western desert. Hutchins et al. [11] predicted the universality of these structures in laboratory scale and an atmospheric surface layer. In particular, Ref. [11] also illustrated the similarity of neutrally-stratified atmospheric surface layer flows and canonical wall-bounded turbulence and the formation of large-scale counter rotating roll vortices in both flows. Further details of the roll-cells in planetary boundary layers can be found in the work of Etling et al. [8] and Yong et al. [9]. In the presence of wind farms, the ABL modulates such structures in and around the turbine rotors as illustrated in the field [12], laboratory-scale experiments [7,13] as well as LES simulations [2,5,14]. Proper orthogonal decomposition (POD) [15][16][17], has been popularly used by the community to perform detailed analysis of these large scale roll structures [15], because of its inherent nature to "rank" the coherent turbulent eddies based on their kinetic energy content (POD eigenvalues). In the current manuscript we are primarily interested in the POD analysis of a finite-sized wind farm where the flow is heterogeneous in the streamwise direction. Finite sized wind farms have multiple applications, e.g., fundamentally understanding the entrance region of very large wind farms [18] without being restrained by computational expense. These finite sized turbine arrays also find their place in utility scale farms e.g., in a distributed wind setup aimed for powering rural or remote sub-urban community [19]. One of the key problems, in the snapshot-based POD analysis [20,21] (more popular because it is computationally cheap), is that the snapshots need to be quite far spaced (more than three flow though times apart), such that the snapshots are linearly independent of each other for the generation of the correlation matrix in POD, and hence capture the large scale structures/modes (having large decorrelation times). Restrained by the computational expense, the number of snapshots in such a case involving finite sized wind farm canopy will be quite limited, which will detrimentally impact the scaling and the convergence of the POD eigenspectra. Previous work by turbulence/wind community have not addressed this issue explicitly [22]. This problem was circumvented by performing the analysis on a periodic wind farm [15]. In this way even though very few distant spaced (3 flow through times apart) snapshots are generated by solving the Navier-Stokes equation, the number can be artificially amplified by "shifting" the snapshot data across the streamwise and spanwise distance between the turbines. The shifting method only works if the wind farms are arranged in a completely matrix setting (n x × n y turbines) in a horizontally periodic framework. In this case the total number of snapshots generated for POD would be n x n y times the snapshots obtained from simulations. In the finite sized wind farm, as in the current manuscript, due to the heterogeneity of the farm layout and the geometry of the domain, such shifting is not possible, and at most we can increase the number of snapshots by two times by reflecting the data across the xz plane of symmetry passing through the middle row of turbines. Thus the data from finite-sized wind farm simulation serves as a perfect candidate for performing such POD analysis studies. As we will see later, this creates a limitation on the total number of snapshots we can generate that impacts negatively on the convergence of the POD eigenspectra. Additionally after performing detailed analysis of the eigenspectra and understanding the dynamics of the 3D eigenmodes, we adapt a Fourier-POD (FPOD)-based framework for better understanding of the eigenspectra (most importantly their convergence) as well as the modes. In this methodology, we perform POD analysis of laterally/spanwise Fourier transformed snapshots rather than the 3D snapshots themselves. As will be discussed in the subsequent section, we see that the FPOD modes can provide better insights to the scaling laws of the eigenspectra as well as elucidate on the turbulent structures which contribute to the problem of "slower convergence" of the spectra. Ideas for computing POD with a similar spirit exists in the literature, e.g., spectral POD (SPOD) by Towne and coworkers [23], where the temporal Fourier transform of the correlation function is used for POD eigendecomposition problem in order to understand the modal frequency content. Towne and coworkers also showed the connection between SPOD and dynamic mode decomposition (DMD) [24]. Recently, Zhang et al. [17] have utilized SPOD analysis specifically for the wind farm problem. Hamilton and coworkers [25,26] have introduced methodologies like dual POD as well as 2D-POD at each streamwise location in an effort to understand the streamwise dependencies of the coherent structures for wind farms. Glegg and Devenport [27] have also shown that for a problem involving turbulence-acoustics interference, POD can be performed in combination with linear stochastic estimator in the homogeneous direction (space or time for homogeneous/stationary flows) to describe flows and that such analysis requires less modes for a flow representation than the conventional 3D POD. The methods described above have utilized POD such that additional information of length scale dynamics and/or decoupling streamwise and spanwise length scales is possible. The concept of Fourier-POD itself is not uncommon in the turbulence community and can be found in many canonical flows having a homogeneous direction, the most recent being Rayleigh Benard convection [28,29]. Contrary to the previous methods, the FPOD methodology (operated on complex-valued two dimensional fields) and its novelty in the current work lies in its introduction to: (1) gaining physical insight to the convergence of the 3D eigenspectra and (2) showing a comparison of the reconstructed physical modes from both FPOD and 3D POD decomposition. Additionally, FPOD also helps us decouple length scales of the energetic eddies based on the spanwise and streamwise direction of the flow.
The aim of the current work is to understand the dynamics of the coherent structures in a finite-sized wind farm using POD analysis. In this context, we first evaluate a traditional 3D POD methodology in the form of a method of snapshots, and then propose an improvement to the traditional methodology, the so-called Fourier-POD method, which demonstrates a faster convergence of the eigenspectra owing to an imposition of analytically optimum modes (Fourier modes) in a periodic direction. The second objective of the manuscript is to propose a physical interpretation of the FPOD modes and compare them with the corresponding 3D-POD results. This work would help us better assess the potential of reduced order models (ROM) for a compact representation of the complex wind farm flow physics, and their ultimate utilization for a forecasting the flow and power generation in heterogeneous wind farms.
The paper is organized as follows. First, we briefly discuss the numerical setup with regards to the large eddy simulation (LES) framework used for the wind farm simulations. Second, we introduce the mathematical formulation of both 3D-POD, and a novel 2D Fourier-POD analysis framework. Subsequently, we present the results of 3D-POD and 2D FPOD methodology. Finally, we conclude by highlighting the merits of a 2D Fourier-POD methodology in understanding the physics of wind farm flows as compared to its 3D POD counterpart.

Numerical Setup
The database for the finite sized wind farm was obtained from spectral element simulations in LES carried out in our previous work [5,30,31]. In particular, we used the open-source exponentially accurate spectral element code Nek5000 [32,33] for setting up the LES simulation. Nek5000 solves the incompressible Navier-Stokes equation in a variational/weak formulation with tessellating the domain into 3D hexahedral elements. The variables, velocity, pressure etc., are expanded as higher order Legendre polynomials within each element and the gridpoints where the polynomials are defined correspond to the roots of the polynomials. For the velocity, Gauss-Lobatto-Legendre (GLL) points were used within each element while for pressure, Gauss-Lobatto (GL) points were used in the current formulation.
The domain consisted of a 3 × 3, finite-sized, wind turbine array (WT) layout in an inflow-outflow [30,34] setup. (Please see Figure 1 for the 3 dimensional schematic of the layout). The inter-turbine streamwise and spanwise distances were 7d and 3d respectively. The wind turbine rotors in the finite sized array, d were set to be 20% of the ABL thickness, H: d = 0.2H, and the hub-height being z hub = d = 0.2H. Such a setup of hub-heights is consistent with the wind-tunnel laboratory scale wind turbine tests carried out in the past [7,35]. The vertical height of the domain was H = 5d. The LES framework involved wall-modeled large eddy simulation with the subgrid scale closure designed as an algebraic wall-damped Mason and Thompson model [36][37][38] and the bottom wall boundary condition prescribed as a wall shear-stress (corresponding to the log-law of the wall). The validation of the wall-modeled LES framework for the precursor neutral ABL [5] (used as an inflow boundary condition to the WT array) compared against experimental data is illustrated in Figure 2. The streamwise velocity, streamwise variance and the integral length scale for the neutral ABL showed expected trends when compared against experimental data obtained from both wind tunnel [39,40] and field experiments [41], as well as high-fidelity numerical simulations [42]. Additionally, our LES model tested in a neutral ABL framework manifested correct logarithmic trends of the velocity profile at the log-layer and inverse scaling laws of the energy spectra [37] at large wavenumbers.   [42] (Re τ ∼ 2 × 10 3 ) and LWP-Leipzig wind profile, Lettau (1950) [41].
For more details related to the validation of the LES framework of neutral ABL and wind farm simulations in spectral elements, please refer to the authors' previous work [5,43]. The wind turbines are modeled as actuator-line (AL) models [39,[44][45][46] without having to resolve the blades. In the current AL simulations, Tjaereborg turbines have been used [44]. The lift and drag coefficients for the blades were taken from the NACA 44xx series aerofoil wind tunnel measurements performed on NACA airfoils at a chord Reynolds number, Re c = 6 × 10 6 . The equivalent time-varying thrust coefficient of the turbines is in the range of C T ∼ 0.7-0.8, which is similar to C T = 0.75 used in the actuator disk model of [2]. Each actuator line is discretized using 30 uniformly sized blade elements (i.e., per rotor radius) as in [44]. The AL model is of a higher fidelity than the actuator-disc model ( [2,15,47]) commonly used in numerical computations of wind farms, in its capability to capture the tip vortices being shed in the near wake ( [44,45]).
Moreover, the finite-sized layout is inherently heterogeneous due to the streamwise growth of the internal boundary layer, wake impingement effects etc. The top boundary condition of the finite-sized wind farm array is "symmetry" or no mass transport in nature which mimics an inversion layer as in a conventionally neutral boundary layer [48,49]. For the current setup, the artificial inversion layer is five times the hub-height (equivalently, five times the rotor diameter) of the turbines and can be thought of as low-lying [48]. From the perspective of canopy turbulence, it is apparent that the ABL above the wind turbine canopy is essentially the "canopy sub-layer" (vertical domain size 5d = 5z hub ). (The canopy sublayer is where the turbulent eddy effects of the canopy are still felt and is roughly five times the canopy height or even larger, see [50][51][52] for details of recent canopy/roughness sublayer fundamentals). Table 1 shows the specification of the domain size of the wind turbine array as well as the precursor neutral ABL which is concurrently simulated to generate inflow conditions [30,31]. Additionally, the grid counts and the average normalized grid sizes for the farm layout as well as neutral ABL (turbulent scales fed to the inflow) can be found in Table 2. This illustrates that the smallest resolved length scales is an order of magnitude smaller than the turbine rotor size, d. Note, here that the definition of lengthscales or wavelength resolved is based on the Nyquist limit of the coarsest grid size (twice the coarsest grid size). For more details of how the grid sizes were defined, see [5,46]. From Table 1 it is apparent that the precursor neutral ABL has a much uniform distribution of grids compared to the wind turbine (WT) domain. This is primarily because of the fact that the WT domain grids were build on top of the ABL domain, with grid refinements around the turbine rotors (∼30 gridpoints per actuator line blade) for capturing the wake turbulence accurately [30,34,46,53]. Table 1. Details of the computational grids for the precursor ABL and the wind turbine array domain. N e η is the number of spectral elements in the η direction. Eight Gauss-Lobatto-Legendre (GLL) nodes (Legendre polynomial order 7) have been used per element per Cartesian direction.

Case
Geometry Table 2. Maximum, minimum and average wavelengths captured, for the ABL and the wind turbine array domain. d is the turbine rotor diameter.

Neutral ABL WT Array
Direction λ η,res max λ η,res min λ η,res λ η,res max λ η,res min λ η,res The number of snapshots obtained from the simulations was 3285, which were spaced 1/5T e apart (T e = 15πd/U ∞ is the flow-through time). Since, the domain/layout is symmetric about y = 2.5πd, we created 2 × 3285 = 6570 snapshots by reflecting and shifting of the snapshot data about xz plane at y = 2.5πd (similar to the shifting method by Verhulst et al. [15]).
In the subsequent section we discuss the mathematical details of the 3D POD and the Fourier-POD methodology that we use in our analysis in the manuscript.

3D POD-Method of Snapshots
The POD analysis was carried out using the method of snapshots [20] developed in the spectral element code Nek5000. The 3D velocity vector field was represented as is the space of square-integral functions in the 3D real space, for which the L 2 , or energy norm, can be defined. The velocity field can be decomposed into a set of orthonormal basis functions where the turbulent velocity fluctuation field u (x, t) = u(x, t) − u(x, t) T , and u(x, t) T is the time average of the velocity field, (ϕ i , ϕ j ) = δ ij ∀i, j. The POD problem can be cast as a constrained variational problem (See [54] for details), with the maximization of the objective function , T is a temporal averaging procedure. A necessary condition of the extrema dictates that the functional derivative vanish for all variations using temporal correlation instead of a spatial correlation. Mathematically, the POD method of snapshots [20] of the velocity field dataset arises when solving for the projection of the dataset P r : V → V r of fixed rank r, minimizing the error N t −1 ∑ j=0 ||u j − P r u j || 2 in the least-squares sense with the constraint ||ϕ|| = 1 (||.|| is the norm corresponding to the inner product (, ) ∈ V). The temporal snapshots of the velocity field u(x, t j ) ∀j = 1, . . . , N t have been written as u j in the error expression for brevity. The projection P r can be written as The correlation matrix in indicial notation can be obtained from the inner product of the snapshots and is given as This method ensures that the eigenvalue problem arising is independent of the size of V which is equal to Ω(R 3k ) (k: number of coordinates in u(x, t) at discrete grid points). In the eigenvalue problem in Equation (4) below, Λ is the eigenvalue corresponding to the turbulent kinetic energy and v is the eigenvector. [ The POD eigenmode can be constructed from the eigenvalues and eigenvectors as for some coefficients b k . Using Equation (2) for u (x, t j ), Equation (5) can be expanded as Since, (ϕ k , ϕ l ) = δ kl , it is straightforward to see from Equation (6) that b k (t j )a l (t j ) = δ kl /N t ∀j. Mathematical details analysing the relation of b k , a k can be found in [15]. It is important to note, that Parseval's identity can be applied in POD (orthonormal basis functions) which from the inner is consistent with the orthonormality of ϕ.

Fourier-POD Methodology
For the Fourier-POD methods, analogous to the projections methods of 3D POD, we can write the projection P f ,r as In a similar spirit, the correlation matrix can be written aŝ whereû is the Fourier transform in the lateral/spanwise y direction of the velocity fluctuation vector Since the snapshots in FPOD methods are complex valued rather than being real-valued as in the 3D POD case, the modes are expected to be complex in nature as well. Note, however that the correlation coefficient (3D POD as well as 2D Fourier-POD) are based on the inner products of snapshots. If χ is a matrix containing columns of snapshot χ i = [χ 1 , χ 2 , χ 3 , . . . , χ n ], then the correlation matrix can be written as χχ T while for complex snapshots (e.g., 2D FPOD), the inner product can be written as χχ H , where [] T and [] H represent transpose and Hermitian (complex conjugate transpose) respectively. Subsequently, it is apparent that while the eigenvectors and hence the modes can be complex for FPOD, the eigenvalues (or the diagonal eigenvalue matrix) are given by diag(λ) = ΣΣ H (since svd can be complex for real systems ), where Σ is the matrix containing singular values at diagonal entries from the SVD of the snapshot matrix χ. Consequently, the eigenvalues of 3D POD as well as complex 2D Fourier-POD are always real-valued and represent kinetic energy content the mode. The FPOD analysis has been carried out in an open-source python code MODRED [56] written as a high-level class interface using spectral methods.

3D POD
In this section we present the results obtained from the POD analysis of a finite-sized wind turbine array. Unless otherwise mentioned, the results (also shown in the plot labels) were normalized with the free-stream velocity scale, U ∞ and rotor diameter d as deemed necessary. The snapshots for the layout were each T e /5 (T e is the flow through time) snapshots apart, which were much frequent than the snapshots "3T e apart" as reported in Ref. [15]. This essentially means that the current POD analysis were carried out in the framework of "smaller time scales". The snapshots ∼ 3T e apart in the previous literature [15] in the context of asymptotic wind farms ensure that the temporal autocorrelation of the velocities completely decay to zero. In the current manuscript we resorted to snapshots separated by time T e /5, when the velocity correlations were roughly ∼0.2. However, we also used snapshots which were 2T e /5, 3T e /5 and 4T e /5 times apart to test the convergence of the eigenspectra. The time extent of 600T e was used for the analysis of the LES database. Note, that during convergence study, we chose to keep the time extent of the database for POD analysis fixed and only varied the time spacing between the two snapshots. Ref. [15] performed POD analysis on periodic wind farm layout, and hence after generating a set of snapshots, 3T e apart, artificially increased the number of snapshots by an order of magnitude by exploiting the stationarity of the flow and the large scale statistical symmetry of the flow around every turbines in the streamwise and spanwise direction (method of "shifting"). Figure 3 shows the normalized eigenspectra of the finite-sized farm layout compared with the spectra of a periodic 8 × 6 array with domain size 20πd × 10πd × 10d (2πH × πH × H) [5], 3D POD modes were calculated using one shifting due to a spanwise symmetry to increase the number of snapshots. (a) (b) Figure 3. (a) Eigenvalue spectrum normalized by the first eigenvalue (highest energy) for different wind turbine layouts. Lines: Blue-periodic wind farm (asymptotically infinite) with 8 × 6 array, 20πd × 10πd × 10d [5], Magenta-current finite sized wind farm (3 × 3 array, 3πH × πH × H). Black circle-periodic wind farm (4 × 6 array, 10πd × 10πd × 10d), LES simulation by Verhulst et al. [15]. Red triangle-periodic wind farm (16 × 12 array, 40πd × 20πd × 10d), LES simulation by Zhang et al. [17]. Blue square-developing wind farm (12 × 12 array, 40πd × 20πd × 10d) [17]. H is the boundary layer thickness. (b) Spectral convergence of the eignevalue spectrum λ m /λ 1 for the current finite sized wind farm. Snapshot separation, black-4T e /5, red-3T e /5, blue-2T e /5 and magenta-T e /5.
The eigenspectra of the 3D POD modes was validated for the periodic and developing wind farm simulations with the results from the previous literature in Figure 3a. In particular, we considered the periodic wind farm layout from Verhulst et al. [15], and periodic as well as developing large wind farms from Zhang et al. [17]. The geometry for all three reference simulations is summarized in Table 3. The developing wind farm geometry of Zhang et al. [17] was the same as the geometry of their periodic wind farm simulations cited in Table 3, albeit with a lesser number of turbines (12 versus 16) in a streamwise direction. Note that the inter-turbine distance was the same for all three references [5,15,17] but different from the current finite-sized wind farm study. Additionally, the turbine diameter d was set as d = 0.1H, and the hub-height z hub = d = 0.1H, for all the three studies [5,15,17].
We note that our periodic wind farm case matches better with the periodic case of Zhang et al. [17], underscoring the fact that both these studies performed the POD eigendecomposition with similar number of snapshots. We also illustrate the convergence behaviour of the eigenspectra for the finite sized farm in Figure 3b. Note, instead of comparing the eigenvalues for different POD problems λ m , we compared the normalized eigenvalues λ m /λ 1 (m is the number of mode), which ensured better convergence of the scaling trends.The current convergence plots were shown for a POD computed for a same time extent of 600T e , but the snapshots were spaced for different fractions of flow-through times apart. However, similar trends in convergence could be observed when performed with the POD for fixed spacing of the snapshots (T e /5) and progressively adding more snapshots with each case (results not added for brevity). The eigenspectra and their convergence trends clearly indicated the following features, (i) convergence is observed by increasing N t , the number of snapshots, (ii) with increasing number of snapshots, the scaling law of the modes for m > 10, clearly changed from m −0.5 [16] to m −0.8 which is extremely close to m −0.9 as observed by the periodic wind farms in [15]. For the same number of snapshots, Zhang et al. [17] also illustrated the m −0.8 scaling for the m ∼ 10 1 − 10 2 for periodic farms and m > 10 for developing farms. Additionally we can also comment that obtaining the convergence of m −0.9 scaling law iwas primarily dependent on the number of snapshots rather than the spacing of the snapshots. Incidentally, a similar scaling law m −1.2 was also observed in experimental studies, particularly for 2D POD eigenspectra performed at different streamwise locations of wind turbine array by Hamilton and coworkers [26]. We hypothesize that with increasingly more number of snapshots, the scaling law for modes m > 10, should approach m −0.9 . However, restrained by the bottleneck of the computational expense we cannot add further snapshots from simulation or take advantage of the symmetric shifting to generate more "artificial" snapshots. Consequently, in order to gain more insights on the scaling law of the eigenspectra and their convergence, we propose a different methodology based on Fourier-POD analysis, which essentially deals with the POD of the snapshots of the lateral/spanwise Fourier transformed velocities for different wavenumbers instead of the 3D snapshots themselves. We provide a dedicated section (Section 3.2) for introducing the mathematical formulation and the results and insights of the FPOD analysis, but for now we continue our discussion related to the results of the 3D POD analysis.
The first 3D-POD mode of the periodic wind farm case (8 × 6 array [5]) is shown in Figure 4 for validation purposes. The mode illustrates 2 pairs of counter-rotating roll-cells (global sweeps and ejections) as the most energetic features in the periodic farm. Similar roll-cell features were also observed in the periodic farm POD computations in [15,17]. Verhulst et al. [15] noticed that the gap between the roll-cells depends on the global dynamics of wall-bounded turbulence and is independent of the domain size. We observe the most energetic mode in [15] had 1 pair of counter-rotating rolls for a domain size of 10πd × 10πd × 10d, and 3 pair of rolls (or features reminiscent of rolls) in [17] for a domain size of 40πd × 20πd × 10d. Note, the differences in the energy contribution of the modes having 1 and 2 pairs of roll-cells were 1%, as reported by [15]. Table 3 describes the roll shape and the features observed in the POD computation of the periodic wind farms in the literature as well as the current authors [5], manifesting good similarity in their topology. Additionally, we also observed that the roll cells approximately covered around two turbine columns similar to the observation by [17], and the value of the streamwise modal velocity was on par with [15]. Analyzing the number of roll cells per domain width and their size obtained for different simulations [5,15,17], we can conclude it was unlikely that the domain width played a role in determination of these quantities, similar to the conclusion of [15]. Again, our simulations [5] were closer to that of [17] rather than that of [15], due to a closer similarity in the number of snapshots used in both [5] and [17] studies. The fact that the energy difference between the mode-1 and mode-2 (in terms of the number of counter-rotating vortex pairs) was less than 1% further testifies to a relatively equal importance of these two most dominant modes in the wind farm energetics. Table 3. Roll features in POD computation of the periodic wind farms, d = 0.1H. Roll-size denotes spanwise extent from center to center of roll vortices (zero crossing of streamwise modes). Roll density is number of roll pairs/L y . Roll-height-height of the center of roll-vortices above ground. Here we discuss the results of the finite-sized wind array. Figure 5 shows the frontal (yz plane) picture of the POD modes for the first eight eigenvalues. In particular, the figures illustrate the colour contours of the streamwise modes overlayed on the top of in-plane (vertical, spanwise) modes represented as vectors. The modes clearly manifested circulatory features, reminiscent of the counter-rotating roll-cell structures and the downdrafts and updrafts of these circulations coincide with the positive (higher energy) and negative (lower energy) streamwise modal structures. Interestingly, the finite-sized layout also showed the most energetic mode containing two pairs of roll-cells similar to the periodic computations by the authors [5]. Additionally, the roll-cells in the finite-sized turbine array also manifested the "mode-pairing" feature [15] where two modes with the same number of roll-cell paired with similar energy content can be observed. The number of roll-cells observed in the layout was related to the modes and hence the eigenvalues/kinetic energy of the flow domain. The height of the roll cells for m = 0, 1 mode was ∼ 2.5d = 0.5H, which was similar to ∼ 4d = 0.4H observed in periodic wind farms. Additionally, we note that the height of the roll-cells was not constant and varied with the mode rank and the number of roll-cells in the layout. Figure 6 gives a 3D perspective to the circulatory roll-cell features we discussed above, in particular, mode m = 0. The isosurface of the streamwise velocity modes indeed showed the long counter rotating roll-cell feature spanning the whole domain. This is corroborative of the fact that these roll-cell features are a phenomenology of the atmospheric boundary layer turbulence [11,15]. The streamtube picture depicted in Figure 6b further illustrates the three-dimensional nature of the large-scale structure manifested by the most dominant POD mode. The three sets of the circulatory features depicted by the streamtubes are formed in and around the 3 columns of the turbines. In order to get a better understanding of how the wind farm modulates these large scale structures we took a closer look at the streamwise variation of the modes in Figure 7 illustrated as the contour plots of the POD mode velocity magnitude. As we will see later that the streamwise plots also assisted us to make a one-to-one comparison with the Fourier-POD modes discussed in the subsequent sections. We observed the streamwise (almost) homogeneous streaks which were essentially footprints of the roll-structures for lower modes. Interestingly, in those footprints we could observe "wake like features" which were clearly manifestation from the wind turbine array. For modes m < 3, we observed that the wake like features (velocity-deficits) extended for scales ∼ 7d (inter-turbine distances) embedded in the roll-cell footprints. Particularly, for mode m = 3, the "wake" footprints could be conspicuously observed with high-velocity regions near the turbine wake rotors and is a clear manifestation of the modulation of the flow structure modes by wake-turbulence. This is reminiscent of the generation of turbulent kinetic energy at the wakes due to vertical entrainment at the inner layer. The mode m = 3 was a clear example of the mode that was entirely due to wind-farm-ABL interactions and would not form in a pure ABL case. Apart from the rolls, at higher modes (m 7), we also observed large structures starting from the wall and inclined at an angle of 15-20 degrees [57,58]. We believe that these inclined structures were footprints of clusters of hairpin vortices which together formed the framework of "attached eddies" [59,60].  In the next section, we discuss the results from the Fourier-POD methodology for complex Fourier-transformed velocity snapshots.

Fourier-POD
Before moving into the results related to Fourier-POD analysis, we present the spanwise Fourier energy spectra of the velocity snapshots ( Figure 9). The spanwise spectra (averaged in the streamwise, vertical direction as well as time) can be given as where x, y are the streamwise and spanwise direction respectively. Note the presence of the λ 11/3 y or k −11/3 y law near the tail of the spectrum. This appears to be a filtering of the k −5/3 y spectrum, as G(k y )k −5/3 y , with G(k y ) = k −2 y spectrum. Such filtering is possibly an outcome of vertical averaging of the spectra, streamwise averaging of the heterogeneity owing to wind turbine array and is mainly an observational documentation. Note for performing the Fourier transform a total of 512 uniform grid points are chosen in the spanwise/lateral direction which ensures a total of lateral resolved length scale (Nyquist limit) of 5πd/256 = 0.06d, which is lower than what the LES simulation with Navier-Stokes solver could resolve (1). Such a high number of grid points for the FFT ensures that the aliasing error could be avoided during the reconstruction of the 3D POD mode from the 2D complex Fourier-POD modes. Figure 9 not only indicates energy cascade from the larger to the smaller length scales in the spanwise spectra, but also illustrates that the larger length scales were associated with higher degrees of uncertainties (larger decorrelation time scales) than their smaller scale counterpart. Figure 10 manifesting the Fourier-transformed velocity magnitude serves as a sanity test of the Fourier-transform itself, with the snapshot at k y d = 0 representing the temporal snapshot of spanwise averaged velocity magnitude, while the same for k y d = 2, manifests wake-imprints from finer scale fluctuations. From Equations (2) and (7), we can define the complex velocity magnitude mode from the POD analysis as ||u|| = |û| 2 + |v| 2 + |ŵ| 2 and |û i | = [R(û i ) 2 + I(û i ) 2 ] 1/2 , ∀i = 1, 2, 3... Additionally, the phase of the complex POD mode can also be defined as ψ u i = arctan(I (û i )/R(û i )) Figure 11a illustrates the eigenspectra of the 2D complex field snapshots (spanwise Fourier transformed velocity) documented at different wave numbers. Expectedly, we observed that the eigenspectra had a slower decay with modes as we moved to higher wavenumbers or smaller spanwise length scales (laterally thin structures). This manifested by the spectra getting flatter for higher wavenumbers. This indicated that for turbulent structures which were laterally thin (higher k y ), the energy content did not vary significantly based on the POD ranks and approached towards a uniform distribution of eigenspectra at k y → ∞ (isotropization of small scales). Figure 11b shows the convergence of the Fourier-POD modes for a representative wavenumber of k y d = 2. While similar results can be obtained for other wavenumbers, it is striking to note that the FPOD eigenspectra converged much better than the 3D POD eigenspectra at all modes. Since the number of snapshots and the spacing between the snapshots are the same for the 3D POD and 2D Fourier-POD analysis, the only way this was possible was when most of the uncertainty manifested at the larger length scales are contained in the Fourier energy spectra as is evident in Figure 9. It is worth noting that the Fourier-modes are analytical descriptors of the POD modes with periodic boundary conditions [55]. Consequently, a FPOD type of decomposition for a wind farm with spanwise periodic boundary conditions, by construction, picked up the correct spanwise modes without running into convergence issues. The 3D POD modes were supposed to converge to the FPOD modes asymptotically with the increase in the number of snapshots.
In Figure 12a-c, we illustrate the behaviour of the FPOD eigenspectra for lower and higher wavenumbers. What is striking, is that we could capture a m −1 scaling law for modes spanning a decade from m = 10-100 for wavenumbers k y d < 2.0. For the number of snapshots considered, we observed a scaling law of m −0.8 in the 3D POD eigenspectra. Interestingly, we also observed this scaling of FPOD eigenspectra at 25% of the total number of snapshots (∼ 1600 snapshots) used in eigen-decomposition, where we were only able to capture m −0.5 scaling for the 3D POD eigenspectra [16]. This was one of our most crucial observations in the whole analysis. We observed that the scaling exponent gradually decayed as we consider the eigenspectra at higher wavenumbers indicating a flatter spectral tail as was discussed above (more isotropization of structures). The results have two implications; (i) The m −1 (m −0.9 [15].) scaling is not an "artifact of the lack of convergence" and is observed in laterally wide turbulent structures for FPOD modes m > 10. Interestingly, similar scaling laws (m −1.2 ) were also noted by Hamilton and coworkers [26] for 2D POD of wind turbine arrays at different streamwise locations. (ii) For FPOD eigenspectra, as the turbulent eddies become thinner, the scaling exponents become larger and tend closer to zero (for higher wavenumbers). Note, in 3D POD we were not able to capture scaling laws of m −1 or m −0.9 , owing to the lack of enough number of snapshots to cover several decorrelation times of the larger scales of interest (higher uncertainties in the larger coherent scales). In the FPOD analysis, we were not only able to decouple the length scales in the spanwise (Fourier) and streamwise(POD) direction, but were also able to decouple the uncertainties-with the Fourier energy spectra carrying most of the uncertainty. Consequently we were able to capture the scaling laws of the FPOD modes (less uncertainty due to better convergence) accurately. This analysis shows that while FPOD expectedly cannot essentially improve the POD results, they have the potential to provide crucial insights to the lack of convergence in the 3D POD eigenspectra and hence their scaling laws. In Figures 13-16, we illustrate the Fourier-POD eigenmodes for wavenumbers k y d = 0, 0.4, 0.8, 1.2. Interestingly, for all modes m = 0, ∀k y , we saw homogeneous/quasi-homogeneous streaks in the streamwise direction which were evidently the footprints of the "roll-cells" modulated by wake turbulence. Apart from the "roll cells", modes m > 0, ∀k y ≥ 0.8 also displayed large scale inclined structures (k y d = 0.8, m = 2-5, k y = 1.2, m = 2-5) which are essentially manifestations of "attached eddies" [6,59]. The FPOD's manifest that the attached eddies form at some threshold spanwise wavenumbers which could not be identified in 3D POD modes. These structures also display similarity to the modes computed from the 3D POD eigen decomposition. Additional to the "roll-cells" and "attached eddy" foot-prints, we observe another type of mode, which are reminiscent of wave like features (k y d = 0.4, m 4, k y d = 0.8-1.2, m = 6). The feature is most conspicuous for k y = 0.4, m = 4. Such modes might be manifestation of wave modulation of large turbulent structures at the particular wavenumber k y = 4, but further studies are needed to speculate the hypothesis. Interestingly, we also observe that as k y d increases (laterally thin structures) the streamwise size of the inclined structures/attached eddies remain approximately the same for a fixed mode number and the streamwise length scale progressively decreases with increasing mode size. This manifests disintegration/cascading of larger eddies to their smaller counterpart. Figure 17 illustrate the phase of the complex streamwise FPOD modes for mode m ≤ 4. The streamwise FPOD modes are dominant compared to their lateral and vertical counterparts. It is interesting to observe that the mode phase ψ changes sign at the edge/boundaries of the large scale eddies and hence is expectedly a great method to visualize the edges of the turbulent structures.

Reconstruction of 3D Modes from Fourier-POD
In this final section we deal with the reconstruction of the three-dimensional modes obtained by performing inverse Fourier transform of the complex Fourier-POD 2D modes. The three dimensional reconstruction of the modes, ϕ j (x, y, z) can be obtained as Note mathematically, ϕ = ϕ (3D POD) as will be established later from the structure of the correlation matrix.
The correlation matrix for FPOD can be written as Note for the inner-product in the physical space ( Ω ()() * dV), the complex conjugate (denoted by * ) of the real velocity snapshot is the velocity snapshot itself. For brevity of analysis, we use the symbol k instead of k y for spanwise wavenumbers in subsequent derivations. Since , we have a relation between C mn andĈ mn as follows.
RealizingĈ mn (k) =ˆ C mn (k, k = k), with the assumption that the temporal correlation of the snapshots are the strongest when the wavenumbers k = k . Equation (14) can be simplified as Note, the integration does not have a closed form. However, the function in the integrand e i∆kL y − 1)/i∆k is bounded between 0 (k max −→ ∞ and finite complex constant L y [1 + i]. Thus, we can comment thatĈ mn (k) is bounded if C mn is bounded as well. Additionally, since the (x, z) spatial structure of C mn andĈ mn are the same (from Equation (13)), the eigenvalue scaling laws for the FPOD modes are similar to its 3D counterpart indicating that the scaling laws are a manifestation of the streamwise dominance of the large scale eddies. We present the reconstruction of the 3D POD modes by inverse Fourier transform (IFFT) of the FPOD modes as illustrated in Equation (11). Figure 18 represents the reconstructed mode by IFFT of FPOD modes, m = 0, m = 1 for k y d = 0.8 (k y = 2). Figure 19 illustrates similar such reconstruction, but for k y = 1.6 (k y = 4). Finally, Figure 20 illustrates the reconstructed modes by performing IFFT of m = 1, for a set of wavenumbers k y d = 0-1.2 (k y = 0-4) and k y d = 0-1.6 (k y = 0-5). From a more quantitative perspective, the reconstructed 3D POD modes from the FPOD that are illustrated in this paper, can be given as Here, W(k y ) is a weight function of wavenumbers which act as a filter to the Fourier transform. For Figures 18 and 19, W(k y ) = 1 for k y = ±2 and k y = ±4, respectively, while W(k y ) = 0 for all other wavenumbers. For  Figure 20 represents summation of the first several most dominant spanwise modes for a particular m. The reconstruction methodology is motivated by the fact that the first two dominant modes (highest contributors of kinetic energy) have two-pairs of counter-rotating roll-cells or four vortical structures corresponding to an "apparent wavenumber" of 2 [15]. The partial reconstructions (via IFFT), shown in Figure 20, further illustrate the importance of the 3-4 pairs (k y = 3, 4) of counter-rotating roll-cells in the dominant mode shapes contributing to the reconstruction, despite having a dominant Fourier mode k y = 2 identified as the most energetic POD mode. Note, also that the partial reconstructions give rise to variable size and shape of the modes as opposed to reconstructions by IFFT involving a single wave number or modes obtained from the 3D POD decomposition. This can be probably attributed to the exponential premultiplier term e i∆kL y − 1)/i∆k referred in Equation (15) that alter the y variation of the FPOD correlation matrix as opposed to the 3D-POD correlation matrix. While examining the m = 0 and m = 1 mode composition for k y d = 0.8 and k y d = 1.6 that correspond to the most dominant 3D POD modes (Figures 18 and 19), an interesting fact can be observed. It can be seen that the modes m = 0 illustrate the global mechanisms of momentum transfer via ejections and sweeps (updrafts and downdrafts) across the entire boundary layer depth, i.e., the global interactions between the inner and outer layer, as 3D POD modes capture as well. However, m = 1 modes illustrate a "bi-layer" structure, corresponding to a momentum transfer between the wind turbine wake region and the inner/outer layer respectively. These modes, which characterize the important energetic mechanisms in wind-farm/ABL interactions, are not picked up by the 3D POD decomposition, while they are by FPOD. The analysis involving the 3D-POD, 2D Fourier-POD and the reconstruction of the 3D modes from FPOD reveals that, while one-to-one mapping of the 3D POD modes and the reconstructed 3D modes by the mode ranks, m, is difficult at this stage, it shows the importance of counter-rotating roll-cell structures involving ejections and sweeps in wind farms and the atmospheric boundary layers in general.

Conclusions
In the current manuscript we have analysed the dynamics of large turbulent structures in a heterogeneous finite-sized wind canopy using three-dimensional proper orthogonal decomposition. Large counter-rotating roll-cell structures as well as inclined wall-attached structures have been identified in this analysis. The current analysis further reveals that substantial number of snapshots are required to obtain the convergence of the scaling trends of the POD eigenspectra, or in particular, the m −0.9 law. In a heterogeneous wind-farm, where artificial snapshots cannot be created exploiting the domain homogeneity (shifting in periodic wind farms) [15], the eigenspectra does not converge well beyond mode m > 10. The lack of convergence is attributed to the uncertainty (higher decorrelation times) in the large scale structures which are still present even in high-order POD modes m > 10 (cf., e.g., Figure 3). Consequently, the scaling trend is slightly deviated to m −0.8 . This led us to adapt a novel Fourier-POD methodology (FPOD), to gain further insights on the convergence of eigenspectra as well as the dynamics of large scale modes. FPOD essentially performs the POD eigendecomposition of the laterally Fourier-tranformed two dimensional complex velocity snapshots at each wavenumber as opposed to the three dimensional physical velocity for the 3D POD. The Fourier-POD analysis helps us gain valuable insights on the convergence of the eigenspectra by decoupling the length scales in the spanwise and streamwise direction. In particular it shows that the laterally wider structures are responsible for the m −0.9 /m −1 scaling laws, while the spanwise thinner structures manifest m −β where β < 0.9. Additionally, we show excellent convergence of the Fourier-POD eigenspectra, indicating that the uncertainty of the larger turbulent scales are mostly contained in the Fourier energy spectra, rather than the FPOD modes. Finally, we look into the reconstruction of the 3D modal structures by performing inverse FFT operations on the 2D FPOD modes. From the mathematical analysis of the functional form of the correlation matrix, we provide deep insights about the similarities in the eigenspectra scaling and the modal shapes of the FPOD and 3D-POD. Eventually, our study reconstructs 3D POD modes from the 2D FPOD modes, which further provides guidance towards the understanding of the modal structure in wind farms. From the fundamental perspective, we have seen that the roll-cells are phenomenologically rudimentary structures contributing to the global sweeps and ejections. While [15] have predicted the contribution of such structures towards the kinetic energy entertainment in infinite wind farms, our studies corroborate such structures to be rather a fundamental property of rough ABL flows [17] and are modulated by the wind turbines at rotor scales. Finally, even though our study was performed in the context of a finite-sized wind turbine array, it introduces a novel framework of Fourier-POD modal analysis, which can be useful for analysis of turbulent flows in other flow domains and configurations, as long as they possess a periodic direction. Future work could include development of reduced-order dynamic models based on either 3D-POD or FPOD decomposition and adapting them to a prediction of flow and power generation in wind farms. Additionally, machine learning techniques [63] could potentially be utilized to generalize the ROMs to varying wind conditions. C mn m, n index of correlation matrix in 3D -POD C mn m, n index of correlation matrix in Fourier -POD Σ Summation ||.|| L 2 norm x min , x max lower, upper limit of the streamwise size of the domain y min , y max lower, upper limit of the spanwise size of the domain z min , z max lower, upper limit of the vertical size of the domain k y Spanwise wavenumber T e flow through time L x Domain size in x direction L z Domain size in z direction L y Domain size in y direction U ∞ Freestream velocity at the top of the ABL U hub Velocity at the hub-height f Temporal frequency E u i Energy spectra of the fluctuating velocity u i T avg Extent of averaging time span, = 600T e R 3 3D real space L 2 (R 3 ) Space of square-integrable functions defined in 3D real space Ω Topological manifold in which the N-S equations are solved T time averaging operator k max Maximum wavenumber used in the Fourier-reconstruction (, ) Inner product W(k y ) Wavenumber dependent weight functions acting as filters to fourier transforms ϕ j 3D j th POD mode vector from velocity snapshot ϕ j 2D j th Fourier-POD mode vector from Fourier-transformed velocity snapshot ϕ j Reconstructed 3D POD mode vector from 2D j th Fourier-POD mode () u,v,w Superscripts used in POD modes denoting streamwise, spanwise and wall-normal direction Reynolds number based on free-stream velocity and chord length C T Thrust coefficient in wind turbine rotor Re τ Reynolds number based on friction velocity