{'fraction_non_alphanumeric': 0.03199123650264917, 'fraction_numerical': 0.013100534305068074, 'mean_word_length': 4.470465941054177, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A three-dimensional dynamic mode decomposition analysis of wind farm flow aerodynamics

High-fidelity large-eddy simulations are suitable to obtain insight into the complex flow dynamics in extended wind farms. In order to better understand these flow dynamics, we use dynamic mode decomposition (DMD) to analyze and reconstruct the flow field in large-scale numerically simulated wind farms by large-eddy simulations (LES). Different wind farm layouts are considered, and we find that a combination of horizontal and vertical staggering leads to improved wind farm performance compared to traditional horizontal staggering. We analyze the wind farm flows using the amplitude selection (AP) and sparsity-promoting (SP method) DMD approach. We find that the AP method tends to select modes with a small length scale and a high frequency, while the SP method selects large coherent structures with low frequency. The latter are somewhat reminiscent of modes obtained using proper orthogonal decomposition (POD). We find that a relatively limited number of SP-DMD modes is sufficient to accurately reconstruct the flow field in the entire wind farm, whereas the AP-DMD method requires more modes to achieve an accurate reconstruction. Thus, the SP-DMD method has a smaller performance loss compared to the AP-DMD method in terms of the reconstruction of the flow field.


Introduction
Renewable and sustainable energy plays an increasingly important role in modern society's development [1,2,3]. Wind energy is one of the fastest-growing forms of renewable energy and will form a central part of the future energy system. It draws extensive attention from academia and industry, not only because of its fast technological advancement but also for business opportunities [4]. Most wind energy is produced in large wind farms in which wind turbines are clustered together. Unfortunately, this leads to wake effects, which can affect the performance of large offshore wind farms by 10% to 20% [5]. Therefore, it is essential to understand these wake interactions and to further develop wind farm control strategies that can be employed to mitigate wake effects in extended wind farms [6].
This work will use high-fidelity large-eddy simulations to study the flow dynamics in extended wind farms. Such simulation cannot be used directly to design wind farm controllers as they are too computationally intensive. However, these simulations capture the relevant flow dynamics accurately, which can be studied using popular flow decomposition methods [7,8] such as proper orthogonal decomposition (POD) and dynamic mode decomposition (DMD). While POD enables one to find the most energetic time-averaged flow structures [8] DMD allows one to capture the flow structures at different frequencies, which allows for a more detailed characterization of the relevant flow physics.
Verhulst & Meneveau used POD analysis of large eddy simulations (LES) to show that large-scale streamwise counter-rotating rolls contribute to the kinetic energy in entrainment of high-velocity wind into the wind turbine array [9]. It has also been reported that the energy of POD modes was associated with the characteristic wavelength of the flow between the turbines [10,11]. Generally, POD modes with high kinetic energy were associated with a large characteristic wavelength in wind farms [10]. Zhang & Stevens [11] used POD to study the flow structures within and above large wind farms. The POD modes of different flow patterns have been identified with physically-relevant frequencies. Hamilton et al. [12] constructed a reduced-order model for wind turbine wake based on POD modes and showed that their model can accurately predict the turbine wake dynamics in terms of large-scale structures.
Due to the enormous amount of data required for a detailed DMD analysis, previous DMD studies were limited to analyzing the flow around a single or limited number of turbines. Premaratne & Hu [13] used DMD to study the data from experimental particle image velocimetry measurements of the flow behind a wind turbine. Their analysis reveals that vortical structures evolve from axisymmetric to asymmetric shape before the vortices break up. A data-driven approach built on Hankel-based DMD was implemented by Ali & Cal [14] to study the wake characteristics behind a wind turbine array, and then a predictive model was built to forecast the flow evolution on small time scales. The model also captures the chaotic transition between the center of the rotor and the atmospheric streamflow caused by wake advection and vortex shedding. Iungo et al. [15] proposed a data-driven reduced-order model to study the wake behind a wind turbine using DMD. They showed that important DMD modes can be connected to, for example, the blade rotation frequency. Annoni et al. [16] and Cassamo & van Wingerden [17] used a so-called input-output Dynamic mode Decomposition (IO-DMD) method in which the DMD system is extended by considering input and output observables. These studies showed that it is possible to construct a model capable of reconstructing wind turbine wakes quite accurately using a limited number of modes, demonstrating the potential DMD approach for the design of wind farm controllers. In the present study, we will analyze the flow of high-resolution wind farm simulations, considering DMD analysis on a much larger scale than previously considered. The analysis will assess the potential of different DMD modes to reconstruct the flow fields in wind farms using a limited number of modes.
The wind farm layout has a strong influence on the wind farm performance and large-scale flow structures in the wind farms. The work by Lissaman [18] in 1979 symbolizes the start of wind farm layout optimization, and Mosetti et al. [19] pioneered the use of wind farm layout optimization using genetic algorithms. Also, various LES studies (see, e.g. Archer et al. [20], Wu & Porté-Agel [21], and Stevens et al. [22], and many others) have shown that horizontal staggering improves wind farm efficiency while significantly affecting the dominant flow structures. The use of vertical staggering to improve wind farm performance is less well explored. However, it is known from experiments (see e.g. Vested et al. [23], Chamorro et al. [24]), LES (see e.g. Vasel-Be-Hagh & Archer [25], Zhang et al. [26,27], Wu et al. [28]), and model calculations (see e.g. Wang et al. [29], Chamorro et al. [30], Stanley et al. [31,32]) that vertical staggering can improve wind farm performance. Besides, Abkar and Porté-Agel [33] and Xie and Archer [34] have shown that the wake expansion in the horizontal direction is faster than in the vertical direction. As the wake recovery is slower in the vertical direction, it is more effective to displace wind turbines in the vertical direction such that they are outside the wake of preceding turbines, although the vertical displacement of turbines is limited by the turbine height that can be achieved. Therefore, we will perform LES of different wind farm layouts to assess the ability of DMD to capture the flow structures in various cases. We consider four different wind farm layout, exploring the effect of horizontal and vertical staggering, which also reveals that even for large wind farms combining horizontal and vertical staggering is beneficial.
To summarize, previous DMD studies have focused on the flow structure around a single or a limited number of turbines. In contrast, in the present study, we focus on the ability of DMD to capture the essential dynamical modes in extended wind farms, thus providing an assessment of the potential of DMD to develop reduced-order models that are suitable to study wind farm dynamics and wind farm control strategies. In Sec. 2 we introduce the simulation method and DMD analysis techniques. The DMD analysis is provided for four distinctly different wind farm layouts, and in Sec. 3 we reveal that a combination of horizontal and vertical staggering is beneficial for wind farm power production. In Sec. 4 the DMD analysis is presented, focusing on the ability of amplitude mode selection (AP method) and sparsity-promoting (SP method) to capture the relevant flow modes and the ability to reconstruct the flow with a limited number of modes meaningfully. The main conclusions are presented in Sec. 5.

modeling methodology and analysis methods
Essential aspects of LES of wind farms are introduced in subsec. 2.1. The DMD application to characterize wind farm flows is discussed in subsec. 2.2.

LES for wind farm flows
We use LES to simulate the flow in large wind farms in a neutral atmospheric boundary layer (ABL). Turbines are modeled using the actuator disk model [35] with the assumption that they operate in regime II in which the turbine thrust coefficient C T is constant [36]. In this work, we assume that all turbines operate at C T = 0.75 [37]. The governing equations that are solved in the LES are: whereũ i represents the resolved velocities, f i denotes the forces from the turbines modeled by the actuator disk model,p * represents the filtered modified pressure written asp * =p/ρ + 1 3 τ kk − p ∞ /ρ. τ i j is the subgrid-scale stress term, whose deviatoric part (τ i j − 1 3 τ kk δ i j ) is modeled by the scale-dependent Lagrangian dynamic model [38,39]. The molecular viscosity term is neglected as we consider very high Reynolds number flows [38]. A pseudospectral discretization is utilized in the horizontal directions combined with a staggered second-order centered finite difference scheme in the vertical direction. On the top of the physical domain, a no stress and no flow-through condition is assumed. At the ground, the Monin-Obukhov similarity theory is employed to calculate the shear stress on the surface [40]. Time marching is performed using a second-order accurate Adams-Bashforth scheme. Thermal effects and Coriolis effects are not considered. The concurrent precursor method is used to generate realistic atmospheric inflow conditions [41,42]. The code [43] has been validated against the EPFL wind tunnel experiments [44,45], the Horns Rev wind farm measurements [42], and high-resolution wake measurements performed at Delft [46].
To assess the impact of the wind farm layout on the flow structures in the wind farm using DMD we consider four different layouts, namely: (1) Alignment in the horizontal and vertical direction;  staggering is achieved by elevating the height of turbines in odd rows and lowering the even rows. The horizontal staggering is obtained by shifting the position of turbines in even rows by half of the spanwise turbine spacing in the lateral direction.
All simulations are performed. on a 512 × 128 × 192 grid for a three-dimensional domain of 4πH × πH ×H in the streamwise, spanwise, and vertical direction, respectively. The domain height H = 1000 m,  [11,47] each wind farm has 108 wind turbines (N x = 18 rows in streamwise and N y = 6 columns in spanwise direction). Both the spanwise turbine spacing S y and the streamwise turbine spacing S x are constant. It is worthy of mention that for the horizontally staggered cases 'H' and 'HV', the effective distance in which the downstream turbines are mainly affected by the upstream wakes is two times the streamwise turbine spacing, e.g. the row 3 is now affected by the wakes of row 1 instead of row 2. H d refers to the difference in hub-height compared to the normal turbine height H t = 100 m. H d is 40 m, which means that the shorter turbines are 60 m and the taller turbine 140 m. Based on the data from Ref. [26] we find that that for the boundary layer under consideration the power output of the tall 140 m and the short turbine 60 m is about 10% higher than that for two turbines of average height 100 m. Table 1 summarizes the simulation parameters. The horizontal computational domain lengths L x and L y are normalized by H, and streamwise and spanwise turbine spacing, i.e. S x , and S y , are normalized by D t .

Sparsity-promoting dynamic mode decomposition
DMD is a powerful tool for flow analysis (see Ref. [48] among others). Like POD, a DMD analysis results in modes that can be used to describe the flow physics. Each mode evolves in time and provides compact dynamic information underlying the original physical system. These modes can be regarded as a reduced-order representation of the extracted structures from a physical domain sampled in a time sequence that contains the largest contribution to the evolution of the system. DMD enables one to cast complex flows onto subspaces that constitute a dynamic system whose degrees of freedom are reduced. Sparsity-promoting dynamic mode decomposition (SP-DMD) was proposed by Jovanovic et al. [47] to achieve a trade-off between the quality of the approximation and the number of required modes. In other words, it provides an objective method for one to determine quantitively at what order the decomposition can be truncated to achieve a given accuracy. We perform three-dimensional SP-DMD on the flow field created in different wind farm layouts in this work. The algorithm is briefly introduced below the following Ref. [47].
Firstly, two data matrices originating from the velocity field measured in temporal sequence are formulated X ∈ R 3N xyz ×N t , X ∈ R 3N xyz ×N t : where N t indicates the number of snapshots; N xyz accounts for the spatial degrees obtained by N x × N y × N z , which are the number of grid points in the three spatial directions. Therefore, each row vector of X, X has size N t , and each column of X, X has a size of 3 × N xyz , in which the factor 3 indicates that three velocity components are considered. Note that X is one step forward than X. Then, one can find a time-independent matrix operator which is capable of optimally and linearly relating the two matrices in time such that: where the A operator can be obtained by solving the optimization problem: In Eq. 5, F is the Frobenius norm and the † symbol indicates the pseudo-inverse, which is computed as where * denotes the complex conjugate transpose, and V , Σ, U can be obtained from singular value decomposition: It is necessary to clarify that directly solving Eq. 6 is computationally intractable since the spatial orders are large. Hence, one typically uses X X = V Σ 2 V * and conducts the eigenvalue decomposition to obtain V and Σ before obtaining U using U = XV Σ −1 .
Observing Eq. 4, reveals that A tends to have a large dimension. To make the computation more tractable, we consider only the first r leading orders by projecting the A matrix onto U. We obtain: where the F DMD matrix is the optimal reduced representation of the A matrix in the U space.
In the r-dimensional subspace, we should formally havex t = F DMDxt−1 , where the reduced-orderx t is related to the original state vector by x t = Ux t = UF DMDxt−1 . F DMD can be expressed in a diagonal coordinate with respect to its eigenvectors { y 1 ... y r } and the corresponding eigenvalues{ µ 1 ... µ r }: where Y and Z are unitary matrices with orthonormal columns, i.e. z * j y i = δ ji . φ i = Uy i indicates the DMD mode, and we can reconstruct the original state as: where α i = z * ix 0 can be considered as the amplitude of each DMD mode. Rewriting the above equation in a matrix form gives: where V and is the Vandermonde matrix of the eigenvalues of matrix F DMD . The determination of amplitude vector α is equivalent to solving the following optimization problem: In order to have a controllable trade-off between the number of modes and the quality of the approximation, the sparsity structure of amplitude vector is modified accordingly. The optimization problem can be rearranged as: where J(α) is obtained by plugging Eq. 7 and Φ = UY into Eq. 12, the second term calculates the L 1 norm of α i facilitating the convex optimization tools and γ is a user-defined positive regularization parameter. The above introduces the specific procedures of SP-DMD. We will consider two different mode selection criteria. This first is known as the AP method, which selects DMD modes based on the amplitude of α i in Eq. 10. Specifically, the DMD modes are rearranged based on their amplitude from the highest to the lowest amplitude. Order reduction is achieved by excluding the lower amplitude modes, while the modes with larger amplitudes are retained. The second mode selection criterion we use is the SP method. As introduced above, one can firstly raise the user-defined γ to penalize the non-zero elements in the amplitude vector by solving a convex optimization problem. In other words, larger γ results in more zero elements in the amplitude vector, which increases its sparsity. The global solution of the convex optimization problem can be solved by the alternating direction method of the multipliers algorithm. The remainder of the non-zero elements is adjusted by optimization to approximate the original space as closely as possible. As such, one can iterate the above process by inputting a set of γ ranging from low to high values and obtain a set of amplitude vectors ranging from a small number of zero elements to a large number of zero elements. However, it should be noted that the non-zero elements in different amplitude vectors might be different because after the sparsity is increased, the values of the remaining non-zero elements are optimized. For more details, see Ref. [47]. Fig.2 briefly summarizes the flow chart of the DMD analysis. Firstly, singular value decomposition (SVD) is implemented on raw data X, where U is the POD mode and Σ is the eigenvalue of POD mode, see Eq.7. Secondly, the mapping operator A is computed as shown in Eq.5 and Eq.6. It is further projected on POD modes to obtain a lower-order representation F DMD indicated in Eq.8. Next, eigenvalue decomposition is implemented on F DMD to obtain eigenvalue and corresponding DMD mode denoted as D µ and Φ respectively, see Eq.9 to Eq.11. Finally, sparse amplitude vector is calculated based on AP method or SP method.

Flow and performance analysis for different wind farm layouts
In subsec. 3.1 we will first discuss the average flow structures for the various wind farm layouts before discussing the wind farm power production of the various layouts in subsec. 3.2.

Flow structures in different wind farm layouts
In the 'A' case shown in Fig. 3 turbines are aligned in the horizontal and vertical direction. The figure shows the dominant streamwise velocity component of the flow, revealing the wake interaction in the streamwise direction. Clearly, the turbines experience strong wake effects in this layout as they are fully submerged in the wake of relatively close upstream turbines. u * = 0.45 m/s is the friction velocity and Take the singular value decomposition (SVD) of X: The matrix A from Eq.4 can be computed: Compute F DMD in Eq.8: Take the eigenvalue decomposition of F DMD to obtain Y and D µ , and DMD mode Φ=UY in Eq.10 Compute the sparse amplitude vector by AP method or SP method has been used to nondimensionalise the velocity in our simulations to get typical dimensional results that may be observed in a wind farm (in the sections below, dimensional velocities and frequencies will be shown and discussed. They can be nondimensionalised by this velocity scale and the relevant length scale). Fig. 4 shows the streamwise velocity profiles for the vertically staggered case 'V'. Note that, in this case, the turbines are horizontally aligned and vertically staggered, so two streamwise velocity profiles are sketched from different-height XY-planes (at z/H = 0.14 and z/H = 0.06). Turbines in odd rows are elevated to 140 m, and turbines in even rows are lowered to 60 m. Comparing Fig. 3(b) with Fig.  4(c) reveals that the flow in the vertical direction is more complex than for the 'H' case as the vertical staggering leads to a zigzag pattern in the average velocity field. Comparison between Fig. 3(a) and Fig.  4(a)(b) suggests that the wake effects in the 'V' case are smaller than for the 'H' case. This is confirmed in Fig. 4(c), which shows that the difference in hub-height between the short and tall turbines limits the impact of the wakes on the performance of downstream turbines. Fig. 5 shows the streamwise velocity in different planes for the horizontally staggered case 'H'. The visualization shows that the wakes from upstream turbines have significantly recovered before reaching  the next turbine. However, the velocity in the entire wind farm area is now affected by wakes, whereas the results for the 'A' case revealed the presence of high-velocity wind speed channels in between the turbines see also e.g. Wu & Porté-Agel [21] and Stevens et al. [22]. Fig. 6 presents the streamwise velocity profiles for the case 'HV'. Since, in this case, the turbines are horizontally and vertically staggered we present the average flow field at two different heights (panels a,b) and two different spanwise locations (panels c,d). The figure suggests that the wake effects are the smallest for this configuration as the velocity just in front of the turbines appears stronger than for the other wind farm layouts. This impression will be    confirmed in the next section, in which we analyze the power production of the various cases.

Effect of layout on power production
To study the effects of the wind farm layout, we calculate the time-average power output as P = − FU d , where denotes time averaging. F = − 1 2 C T ρU 2 d A accounts for the local force used in turbine model. Here, U d stands for the disk-average velocity, and A = πD 2 t /4 is the turbine area. Plugging the expression of F into that of P, we can write P = 1 2 C T ρ U 3 d A, see Refs. [27,26,49,22] for more details. In general, both horizontal and vertical staggering can be used to ensure that wake effects are limited. For the horizontally staggered layout, the downstream turbines are also located outside the wake of turbines in the preceding row. Fig. 7(a) presents the power production of each turbine row. Note that the power output is normalized by that of the first row of the corresponding case. We find that vertical staggering (case 'V' and case 'HV') exhibit a zigzag shape due to logarithmic ABL profile and the vertical staggering, while the curves of the case 'A' and case 'H' develop smoothly. Fig. 7(b) presents the cumulative power as a function of downstream turbine row, normalized with the reference aligned case 'A'. The figure shows that vertical staggering case 'V' performs significantly better than the reference case 'A' in the entrance region. However, the relative benefit of vertical staggering decreases further downstream as the 'V' case gradually asymptotes to the result of case 'A'. The results are in good agreement with that of Zhang et al. [27,26]. The reason is that in the entrance region, the taller turbines can benefit from the stronger winds at higher elevations. However, in the fully developed region, the vertical kinetic energy transfer is the main energy transfer mechanism, and this mechanism does not differ much in the fully developed region between 'A' and 'V' cases [50,22,51].
The figure reveals that the horizontally staggered layout 'H' is more effective in increasing the wind farm power production than the 'V' case. The reason is that in the horizontal direction there is more space than in the vertical direction to distribute the turbines such that wake effects are limited. Interesting, when horizontal and vertical are combined, i.e. the 'HV' case, it leads to a significantly better wind farm performance than using just horizontal or vertical staggering alone. The figure shows that also in the fully developed region the 'HV' case is beneficial over the 'H' case. Even though this is not our primary focus, it is interesting to point out that the 'HV' case layout outperforms all the other wind farm configurations. This signifies that combining the horizontal and vertical staggering (see Fig. 1) can further improve wind farm performance.

Flow characterization by DMD
In subsec. 4.1 we discuss the convergence criteria and the required sampling frequency for the DMD analysis. In subsec. 4.2 we analyze the flow structures obtained using the AP and SP method, and in subsec. 4.3 we discuss the accuracy of the flow reconstruction that can be obtained using both methods.

Convergence and frequency selection
This section presents three-dimensional DMD analyses of the flow structures for the four aforementioned wind farm layouts. It is worthy mentioning that in our DMD analysis (following Ref. [48]), the POD modes are first computed. As shown in Fig.2, the first step to implement SVD on raw data gives the POD modes and the corresponding eigenvalues, i.e. U and Σ respectively (see Eq.7). In order to obtain converged POD modes, it is necessary to determine the number of required flow snapshots by checking the convergence of eigenvalues corresponding to POD modes. To verify this we follow the approach by Zhang & Stevens [11] and Newman et al. [10]. We calculate the L2 norm of the normalized eigenvalue spectrum as follows. where λ i (1 ≤ i ≤ N t ) represents the eigenvalues of corresponding POD modes in U (not to be confused with the eigenvalue in Eq. 9). The relative error is then given by: where the reference value η ∞ is determined using N t = 5000 snapshots. The number of tested snapshots is 50, 100, 250, 500, 625, 1000, 1250, 2500, and 5000, respectively. The snapshots are selected from the available simulation length, which means that the sampling frequency is higher when more snapshots are used.
The result of the convergence analysis is presented in Fig. 8. Panel(a) demonstrates that E(N t ) gradually decreases when the number of snapshots is increased and reaches a relatively small value at 2500 snapshots. The eigenvalues in panel (b) can be interpreted as the kinetic energy in the POD modes [8]. Thus, similarly to our previous works [11], we can see that some POD modes (lower-ranking in i) are associated with higher kinetic energy than higher-ranking modes. Panel (c) shows the cumulative eigenvalues, which reveals that about 1000 POD modes are required to capture 60% of the kinetic energy.
In addition to determining the necessary amount of snapshots to obtain converged POD modes, we also need to ensure that we capture the dynamically important flow physics. Therefore, we compare the sampling frequency ( f s ) with the frequency corresponding to the average travel time of a fluid parcel between two turbines ( f t ) in the streamwise direction. Note that both of them have physical unit Hz. The following shows how to compute them. To compute f t , we first obtain the mean velocity of a fluid parcel traveling through the turbine spacing at the turbine height by taking the spatial average of the velocity in the regions encompassed by the green boxes as shown in Fig. 3 (the conclusion below is not sensitive to the actual size of the boxes). The overall mean velocity in these regions is calculated to be around 6.5m/s. Consequently, we can have f t = 6.5 m/s ÷ 524 m ≈ 0.012 Hz (524 m is the streamwise turbine distance in the aligned case). This frequency is halved for the cases in which the turbines are horizontally staggered. f s is selected based on the Nyquist criterion that the lower bound sampling frequency by which the process can be identified is at least twice as high as its inherent frequency. In principle, more snapshots in a fixed span of the flow data acquires more dynamic information with a higher frequency. However, this would be computationally expensive. Thus, as a compromise, we decided to use 2500 snapshots separated by 150 time steps for the DMD analysis below, which are sufficient to provide converged POD modes as mentioned above. This results in a sampling frequency f s = 1 ÷  Fig. 9. The 100 representative DMD modes expressed by circles with colors. Circles superposed by black plus: modes selected by the AP method; Circles superposed by the red cross: modes selected by the SP method. (a) the first 7 modes selected by the two methods respectively; (b) the first 15 modes selected by the two methods respectively. The color and the size of each mode indicate its amplitude. The mode numbers indicated in the zoom-in panels correspond to the leading order modes determined using SP method presented in Fig.10(c).
times higher than f t , which is high enough to capture the relevant flow physics in our wind farms. In the remainder of the manuscript, we will present the DMD results analysis, focusing on the difference between the AP and SP method. Fig. 9 displays the DMD modes selected by the AP and SP methods. As it is impossible to distinguish all 2500 DMD, we only present the 100 representative modes here. According to the color bar, both the size and the color of the circles indicate the amplitude of corresponding DMD modes; see Eq. 10. We have also used black plus and red cross signs to highlight the first DMD modes selected by the AP and SP methods, respectively. Panels (a) and (b) show the first 7 and 15 DMD modes selected by the two methods, respectively. One can immediately notice that the modes selected by the two methods are different.

Comparison of AP and SP method
Recall in subsec. 2.2 that the AP modes are selected based on the magnitude of the amplitude vector of corresponding element α i (see Eq.10) in the amplitude vector. The elements with higher magnitude in the amplitude vector will be chosen in priority. The AP method reflects a tactic of determining the modes that have strong influences on a system's response resulting from the typical initial condition in the snapshot collection [47,52], see the definition of α i in Eq. 10. Figure 9 shows that the modes selected by this criterion are prone to have a higher amplitude. The modes inside the unit circle are stable, and modes falling outside the unit circle are unstable. We note that some of the larger amplitude modes (meaning that their associated α i values are high) are closer to the origin. This means that these modes have a strong decay rate, which means that their strong influences are only manifest at the initial state of the physical process. Hence, we can learn from Fig. 9 that the AP method tends to select those fast-decaying modes which have a large amplitude. For the SP method the modes are selected based on the optimization constraint, see Eq.13. In other words, the SP method determines which modes should be discarded so that the remaining ones can give the best performance. Based this criterion, the SP method selects modes that reside close to the circle with low frequencies, keeping their influences in the whole physical evolution for a long time. This difference in selecting DMD modes between AP and SP method is more apparent as shown in Fig. 9(b) with 15 modes. The difference between the AP and SP is consistent with Jovanovic et al. [47]. In the remainder of this section, we compare the different flow structures captured by the leading AP and SP modes. Fig. 10(a)(b) presents the leading DMD modes selected by the AP method and Fig. 10(c)(d) those obtained using the SP method. We note that, as mentioned above, the magnitude of the AP modes determines their ordering. The ordering of SP modes is determined by gradually raising γ to see which modes remain. However, we find a set of 7 leading SP modes, which can not be further reduced even when high values of γ (see Eq.13) are used in the SP method. Therefore, here, and in the remainder of the paper, we decide to order these leading SP modes from low to high frequencies. The growth rate and frequencies of the SP medes presented in Fig.11 are given in Fig.9 the zoom-in panels. The differences between the obtained modes are pronounced. Modes selected by the AP method focus on small-scale structures that decay rapidly. In contrast, modes chosen by the SP method tend to present large-scale flow structures. In particular, mode 2 obtained by the SP method captures the mean flow pattern discussed in section 2.1, while the AP method completely ignores this mode. Besides, it is worth noting that the flow patterns captured by the SP method resemble the patterns grasped by POD modes of high energy, which are called streamwise roll structures [9,11]. Some of these leading POD modes are presented in the panels of Fig. 11 for reference. Just as the DMD modes, the POD modes are delimited by turbine rows in the spanwise direction and represent large-scale ABL structures. Besides, we note that POD mode 7 demonstrates the streamwise-varying rolls resulting from kinetic energy redistribution due to wake interaction. The POD result matches well with previous works [9,11,53]. From this point of view, we can understand the reason why DMD modes selected by the SP method are capable of reconstructing the original flow field in the wind farm more accurately than using the AP method. The reason is that the SP-DMD modes capture the long-term flow dynamics in the statistically steady state. That DMD modes resemble POD modes to some extent agrees with Tu et al. [54] who showed that the pattern of low-frequency DMD modes is similar to the leading POD modes. Fig. 12 presents the DMD modespectrum with the positive frequencies. Note that the frequency and growth/decay rate are calculated by Eq. 16 below.
where µ i can found in Eq. 9 and Im and Re indicate that the imaginary and real part of a complex number, respectively. The amplitude of each DMD mode is indicated by the circle size and the color, as shown in Fig. 12. In total, there are 1250 modes (N t = 2500, but we only show half of the spectrum) with the size and the color of the circle (see the color bar) both presenting the amplitude of the DMD mode. The black crosses represent the first 200 modes selected by the AP method, and the blue crosses the first 200 modes selected by the SP method. Again, like the results in Fig. 9, one can clearly see the two methods favor different modes. The AP method tends to select the fast decaying modes with a high amplitude. In contrast, the SP method tends to select low-frequency modes with a low amplitude. Plus, most of the 200 modes chosen by the SP method sit within a low-frequency range below approximately 0.025 Hz. As for the AP method, we can see a wide distribution of scattered frequencies from 0 to around 0.08 Hz. Consequently, it can be concluded that the modes selected by the AP method are connected with the wake dynamics of a high frequency as supported by flow patterns in mode 1, mode 5, and mode 6 in Fig. 10(a) represent the large-scale structures similar to streamwise-varying rolls. Also, the frequency of mode 4 in Fig. 10(a)(b) is close to f t so we can see that its flow pattern is depicted by the alternating positive and negative perturbation whose characteristic length scale is close to S x (streamwise turbine spacing) in the streamwise direction is displayed in XZ-plane of Fig. 10(b).
In this section, we compared the results obtained using the AP and SP methods for the aligned case 'A'. The results for other layouts are attached in Appendix A. It is worthy stating that a comparison of the horizontally and vertically aligned case (Fig.10) with the corresponding vertically staggered case (Fig.A.17) shows that the DMD modes are different inside the wind farm. However, the structure of the DMD modes above the wind farm turns out to be quite similar for the wind farms with and without vertical staggering. The reason for this is that the large-scale flow structures above the wind farm formed in the atmospheric boundary layer are not considerably affected by the vertically staggering in the wind farm.

Flow reconstruction using DMD
This section discusses the flow field reconstruction using the AP and SP methods based on Eq. 10. First, in order to gain a general understanding of how the DMD analysis performs, we reconstruct the fluctuating velocity component using a different number of the DMD mode by SP method, 1%, 10%, 20%, 30%, 50% and 90% of the available N t = 2500 DMD modes. As shown in Fig. 13, with the number of used modes decreasing, we can see a clear trend that when enough modes are used (e.g., 90%), the flow reconstruction is faithful in terms of both flow structure and amplitude. It is encouraging to see that even 10% of the modes can reconstruct the salient features of the large-scale flow structures in the fluctuation part of the velocity field (even though the amplitude of the fluctuation is smaller because the higher DMD modes are discarded). When using an even smaller number of DMD modes (1% modes) as shown in panel (f), one can only get the streamwise rolls residing between the turbine rows, similar to what is obtained using POD [9,11]. Fig. 14 shows the DMD reconstruction of the complete flow field (rather than the fluctuation part) for the AP and SP methods using 986 and 1985 modes out of N t = 2500. In the two cross-sections, we can see that the original snapshots manifest highly turbulent structures. The figure shows that the AP method apparently cannot reconstruct these flow structures with satisfactory accuracy, while the SP After gaining a visual impression of the performance of the AP and SP methods, we further quantify the performance of the methods by determining the performance loss, which is defined as [47] P% = X − ΦD α V and F X F × 100 (17) In Fig. 15, we show that increasing the number of DMD modes reduces the performance loss, i.e., the accuracy becomes higher for both the AP and SP method. When a limited number of DMD modes is retained, see the left end of the figure, the performance loss of the SP method is less than for the AP method. The performance loss of the AP method strongly depends on the number of retained modes. The reason is that the leading AP modes have a large decay rate, which does not capture the long-term system dynamics well. In contrast, the SP method prefers to select modes with a low decay rate, which captures the long-term dynamics better. In agreement with Fig. 14 we find that retaining, e.g., 40% or 80% of the DMD modes in the SP method does not affect the performance loss much. Thus, the SP method can be considered superior to the AP method when a limited number of DMD modes is retained. We can also observe that the results for the four wind farm layouts are similar, even for the 'HV' case in which complicated flow patterns are formed. This indicates that the excellent performance of the SP method in the flow reconstruction is robust and general. Finally, we attempt to use different numbers of DMD modes to compute the power production for the 'A' case by the SP method. First, the flow field is reconstructed using different numbers of DMD  modes (100%, 90%, 50% and 20%) with the mean mode included. Then, U d , the disk-average velocity, is obtained by conducting disk-averaging of the reconstructed flow field for different numbers of DMD modes, based on which the power production is calculated correspondingly. Fig. 16 shows the normalised Raw data 100% modes 90% modes 50% modes 20% modes Fig. 16. Power production normalized by the power of the first row in the original data, compared to that calculated using reconstructed flow field based on different number of DMD modes. Note that the mean mode is included in the reconstruction.
power production by the first row using the raw data, compared to the power production calculated using the DMD-reconstructed flow field. It demonstrates that a limited number of DMD modes allow for an accurate reproduction of the power production data once the mean mode is included in the reconstruction.

Conclusions
In this work, LES is used to study the turbulent flow structures in various wind farm layouts using DMD. The combination of horizontal and vertical staggering (the 'HV' case) has not been considered before, and our results indicate that it generates a higher power production than either vertical or horizontal staggering case. Subsequently, the turbulent flow structures in the wind farm are then analyzed by DMD. It is noted that unlike previous DMD studies, which focused on wake dynamics behind a single turbine or a small number of turbines [16,13,55], we concentrate our analysis on the complex flow dynamics in extended wind farms. We used two methods to select the DMD modes and found that the AP method tends to select modes with a high frequency. The frequencies of the leading modes can be as high as 0.08 Hz and these modes characterize the highly oscillating flow patterns corresponding to the small-scale vortex structures shedding from turbines. On the other hand, the SP method tends to select modes with dominant large-scale flow structures with low frequencies and decay rates. The flow patterns of the modes selected by the SP method are similar to the dominant flow structures identified by POD analysis, see also Ref. [9,11]. However, it is essential to emphasize that DMD, in contrast to POD, can resolve the temporal evolution of the structures.
The temporal information in DMD can be used to reconstruct the time evolution of the flows. To assess the potential of DMD for the development of reduced-order models that are suitable to study wind farm dynamics and study wind farm control strategies, we investigated how well a subset of DMD modes selected by either the AP or SP method can reconstruct the flow field in a wind farm. We find that the SP method provides higher accuracy in global physical space while the AP method offers more information on small-scale structures. It is promising that less than 10% of the SP modes allow for an adequate representation of the flow field in the wind farm. This highlights the great potential of DMD for flow analysis, prediction, and wind farm control. As the flows, we simulated here are highly turbulent, it is reasonable to observe that the nonlinearity makes it difficult to accurately capture the fluctuating velocity field using a small number of DMD modes. It also has to be mentioned that the DMD reconstruction works almost equally well for all the turbine layouts considered, showcasing the robustness of DMD in dealing with the complex turbulent flow structures created in wind farms.
To conclude, our work fills the gap of flow analysis of complex wake dynamics in large numerically simulated wind farms using DMD. There is an enormous potential and need for studies on DMD-based flow control for improving wind farm efficiency or more accurate flow field reconstruction.

Appendix A. DMD modes for H,V, and HV cases
This appendix presents the DMD modes for the wind farms with other turbine layouts (horizontal staggering, vertical staggering, and horizontal-vertical staggering). In the following analyses, we will briefly summarize the similarities and differences between different cases.
Appendix A.1. DMD modes for 'V' case  Fig. A.19(c)(d), which is related to a fluid parcel traveling through streamwise distance longer than the turbine spacing (longer length scales). The patterns of DMD modes chosen by the SP method also resemble those of POD modes of this case in Fig. A.20. What is different is the roll structures in POD modes are re-distributed and no longer delimited by turbine rows, for the turbines are horizontally staggered.

Appendix A.3. DMD modes for 'HV' case
Similarly, in the 'HV' case, the AP method tends to pick out high-amplitude modes while the SP method leads to selecting more persistently influential modes. The DMD modes obtained using the AP method show strongly oscillating flow patterns while the flow patterns determined by the SP method conform with the large-scale flow structures. The flow pattern whose characteristic length scale approximately equals 2 S x , is also discovered in Fig. A.21(a)(b) mode 6. The frequency of this mode is adjacent to that of mode 6 in the 'H' case, which is roughly 0.5 f t .

Appendix B. Flow field reconstruction
The flow field reconstruction results of the cases 'V', 'H', and 'HV' are presented below using the AP and SP method (roughly 40% and 80% of the total N t = 2500 modes). The results from these pictures are consistent with the conclusion in the main text: the AP method tends to capture the original flow field's small-scale structures while the SP method enables one to acquire the global flow field accurately. This shows that the SP-based DMD is capable of reconstructing the different original flow features for different layouts. The velocity deficit of the horizontally staggered layouts is more evident in the reconstructed snapshot as shown in Fig. B.24 and Fig. B.25. Conversely, the wake effects are more evident in Fig. B.23.