Anisotropic double‐Gaussian analytical wake model for an isolated horizontal‐axis wind turbine

An anisotropic double Gaussian (DG) model for analytical wake modeling to predict the streamwise wake velocity behind an isolated non‐yawed horizontal‐axis wind turbine is proposed. The proposed model is based upon the conservation of mass and momentum inside a streamtube control volume. The wake growth rate parameters to distinguish the wake expansion rate between lateral and vertical directions were tuned based on numerical and measurement data of utility‐scale turbines. It was found that the proposed model can give feasible predictions within the full‐wake region under different inflow conditions. In addition, the other analytical models based on top‐hat shape and single Gaussian approaches were evaluated for comparison. The root‐mean‐square error statistical analysis was used to evaluate the performance of each examined model under different flow conditions. In general, the proposed model outperformed the other examined models in all wake region categories, particularly within the near‐wake region and the onset of the far‐wake region, which are beyond the scope of the conventional approach for analytical wake modeling. This advantage gives the potential for the proposed model to provide a better prediction for the wake flow estimation within tightly packed wind farms.

Analytical wake model provides fundamental insight into the physics where their formulations rely on the derivation of fluid flow conservation equations. The flow physics inside the wind turbine wake itself can be roughly divided into two main regions: (i) near-wake region, which is roughly taken as the area just behind the rotor, and (ii) far-wake region, which is the region beyond the near-wake. 39 Inside the near-wake region, the flow field characteristics are strongly influenced by the turbine geometry. Meanwhile, the flow characteristics of the far wake become more dependent on atmospheric and topographic conditions. 40 Since the wake is fully developed in this far-wake region and, in the hypothetical absence of ambient shear flow, the perturbation profiles of velocity deficit and turbulence intensity are assumed to be axisymmetric and have self-similar distributions in the wake cross-sections. These selfsimilar and axisymmetric properties then become the basis of analytical models to predict the wake velocity profiles using either top-hat shape [5][6][7] or Gaussian shape distribution. [8][9][10][11][12][13][14][15][16][17][18][19][20][21] One of the pioneering works in analytical wake modeling is the top-hat shape wake model proposed by Jensen in his technical report 5 and polished by Katić et al. 6 This model was derived by only considering the mass conservation equation. The wake velocity was described in a much-idealized way where the velocity inside the wake region was considered constant. Later work by Frandsen et al. 7 used the same top-hat shape assumption for the rectangular distribution of the wake velocity profiles. Their work also considered the conservation of mass and momentum to encompass velocity deficit for both small and large regular wind farms with rectangular shapes. Although those top-hat models were derived without violating the general principles of conservation equations, in fact, the actual wake shape distribution tends to follow a Gaussian shape function. 39,41,42 By considering this reality, some recently developed analytical models are based on the Gaussian function to estimate the streamwise flow field within the wake region. 10,15 Moreover, some Gaussian models consider the nonaxisymmetric shape of the wake area because the wake expands in an anisotropic manner under different atmospheric conditions. 12,13 However, those mentioned Gaussian wake models are based on a single Gaussian (SG) approach, which may be valid only for velocity deficit estimation within the far-wake region. While in the nearwake region, it was clarified both numerically and experimentally that the wake profiles resemble the DG distribution, which contains two local minima corresponding to maximum lift points along the blades. 43 Hence, it is worth considering the DG shape function to better approximate velocity distribution inside the near-wake region.
To date, the analytical investigations that focused on the wake development within the near-wake region are few. Therefore, further investigations will be needed to elucidate accurately the aerodynamics of tightly spaced wind farms, which might be out of the prediction capability from the standard wake models. Related to the near-wake study, Aitken et al. 44 conducted statistical analysis of wind turbine wake characteristics, including an analytical wake model for the near-wake region using the DG approach. However, the parameter for velocity deficit amplitude was not derived directly but was estimated using data fitting from the Doppler Lidar measurement instead. More recently, Keane et al. 14 proposed an analytical wake model based on the DG approach to predict the wake velocity for all downstream distances by applying conservation of mass and momentum. However, Schreiber et al. 21 found some issues in the original Keane model, primarily related to its inconsistency to satisfy the conservation of mass and momentum. Hence, some previously derived expressions were reformulated. In addition, a new wake expansion function was defined using a streamtube control volume, which is a simplified approach to analyze complex airflow through a horizontal-axis wind turbine (HAWT). Although recent work by Keane 45 improved his previously published DG-based wake model, however, none of those mentioned DG models explicitly consider the anisotropic behavior of wake expansion which is inevitable, particularly under different types of atmospheric stabilities. 12,13 Furthermore, the radial position of local minima of the DG wake profile needed to be tuned case by case, thus limiting its generality.
This study complements the previous study regarding the analytical wake modeling using the DG function 21 to model the streamwise wake velocity behind an isolated nonyawed HAWT. The anisotropic nature of the wake growth rate downstream of the wind turbine is considered, thus facilitating its usability under different atmospheric stabilities. Moreover, the additional empirical coefficient for utilityscale turbines is added to the original formula of far-wake onset 46 to calculate the streamtube outlet position instead of tuning. Meanwhile, the constant value of Gaussian minima r 0 for utility-scale turbines is suggested based on the rootmean-square error (RMSE) evaluation. Finally, the effectiveness of the proposed DG model within the full-wake region was validated against the CFD data sets and lidar measurements of utility-scale turbines. In addition, the prediction results from the top-hat shape and anisotropic SG analytical wake models were also included for comparison.
2 | METHODOLOGY 2.1 | Anisotropic DG wake model A simple approach is needed to consider the energy extraction process of the turbine, one of which is by employing the actuator disc theory. In this approach, the rotor is replaced by a circular disc where the airstream flows inside a streamtube control volume. Figure 1 illustrates this control volume scheme where the airflow in a streamwise direction comes from the inlet section A 0 with the undisturbed velocity U 0 and exits from the outlet section A w with the wake velocity U w where the air pressure has regained its undisturbed condition, p p = w 0 . The actuator disc with its permeable surface partially blocks the incoming flow for extracting the wind energy and causes the pressure drop ∆p over the disc area, proportional to axial thrust force T. The incoming velocity that decreases at the disc due to the thrust force T would be even smaller in the wake section A w . Since the mass flow rate ṁis conserved throughout the flow, the cross-sectional area increases at the disc and continues to expand at the wake area A w . By using the same approach by Frandsen et al. 7 where shear forces on the control volume and the acceleration, pressure, and gravity terms are neglected, in addition, if the pressure at the inlet and outlet of the streamtube is equal, the axial momentum equation on the control volume that assuring the conservation of mass and momentum can be formulated as follows: where ρ denotes the air density. The axial thrust force can also be expressed as the function of a dimensionless parameter for the rotor performance, trust coefficient C T , using the following correlation: is rotor cross-sectional area and D is rotor diameter. The velocity difference between U U − 0 w is the velocity deficit downstream of the turbine ∆U . Suppose the deficit is normalized against the undisturbed flow velocity and the anisotropic wake expansion of the wake profile in Gaussian shape is considered, an expression for normalized velocity deficit can be written as follows: is the standard deviation of the velocity deficit which represents the wake expansion function, and f r y z σ x ( ( , ), ( )) is the Gaussian shape function which represents the wake profile distribution.
In this study, the proposed analytical model for estimating wake development is intended to cover the wake prediction within the near-wake region. Therefore, a DG-based approach for the wake distribution is F I G U R E 1 Streamtube control volume of horizontal-axis wind turbine rotor modeled by an actuator disc SOESANTO ET AL.
| 2125 used. The DG shape function for the anisotropic wake expansion is formulated as follows: where r 0 denotes the radial position of Gaussian minima, which is assumed to be equal for both sides from the wake center at each lateral (y) and vertical (z) direction. In this study, its value was derived empirically from the CFD data and lidar measurements. Meanwhile, σ y and σ z are wake expansion function at the lateral and vertical directions, respectively. By referring to Schreiber et al. 21 for deriving isotropic DG wake model from momentum and mass conservation equations, the following equation for anisotropic DG wake model is obtained using the substitution of Equations (3) and (4) into Equation (1): After solving the integrations from the Equation (5), the following expression can be obtained and Further details regarding the integration steps of M and N can be found in Appendix A. Next, the amplitude function C σ x ( ( )) for the DG distribution, which maintains the momentum conservation can be defined by substituting Equation (2) into Equation (6).
By employing the quadratic formula, the realistic solution for the amplitude function is as follows: , then the solution is approximated using absolute value or modulus of a complex number as proposed by Keane 45 : Thus, for , the expression for normalized velocity deficit as formulated in Equation (3) becomes: C σ x f r y z σ x Δ = ( ( )) ( ( , ), ( )). c 0 (11) In this study, the wake expansion function in the lateral direction σ y and vertical direction σ z are defined explicitly instead of using the same expansion function σ for all directions. Thus, the wake expansion functions from the isotropic DG wake model 21 are reformulated as follows: Since the rate of wake expansion is considered anisotropic, thus the lateral growth rate of the wake in the lateral direction k y and vertical direction k z are set differently. Meanwhile, s out denotes the downstream position of the streamtube outlet where the pressure has regained its undisturbed flow (p p = w 0 ). This position represents the onset of the far-wake region which can be calculated based on the formula proposed by Bastankhah and Porté-Agel. 46 However, the parameter in this formula was derived empirically based on wind tunnel measurements, therefore its compatibility for the actual utilityscale turbine needs to be re-evaluated. This is mainly because the Reynolds number R e in wind tunnel measurements is significantly lower than the actual turbines. Hence, a coefficient c is added to the original formula as a scale factor to fit with the utility-scale turbines. Hence, in this study, the far-wake onset formula for an isolated non-yawed utility-scale HAWT is formulated as follows: where α = 0.58, β = 0.077. In this study, several cases from large eddy simulation (LES) simulations and lidar measurements for the wake flow behind utility-scale turbines (some information can be found in Table 1) were used to find the best fit values of the parameter r 0 and coefficient c using RMSE analysis. The result is shown in Figure 2. Figure 2 shows that r 0 = 0.26 and c = 1.2 have the best fit against the analyzed LES simulations and lidar measurements for wake flow behind isolated utilityscale wind turbines. Hence, these values are used in the present study for all evaluated cases and are also suggested for utility-scale wind turbines. By considering the mass conservation inside the streamtube control volume, the wake expansion at the outlet ε as derived in [21] is used. If we assume that there is no incoming flow into the control volume through its surface, the wake expansion at the outlet may be obtained according to the mass conservation equation within the wake region. Thus, by equating the total mass flow deficit at the position just behind the turbine and at the streamtube outlet m x ṡ ( = ) 2 o u t , the wake expansion at the outlet can be approximated using the following expression: where β, which is the ratio between the initial wake area immediately after the wake expansion A 1 (x = 0) and the rotor area A, is expressed as the function of thrust coefficient C T by The left side in Equation (15) is the mass flow deficit rate at downstream distance x = 0 obtained using the Frandsen model. 7 Meanwhile, the righthand side is derived from the mass flow deficit rate at the streamtube outlet from the DG model. Based on Equation (15), the quantity of ε would be complicated to be defined explicitly, thus the solution is obtained numerically.
After the parameters of amplitude and Gaussian functions are obtained using the previously described formulas, the wake velocity behind an isolated wind turbine can be estimated. Another essential factor that needs to be included in calculating wake velocity is the incoming flow condition. It should be noted that the conservation equations which underlie the proposed formulas are assumed to have uniform inflow distribution. Therefore, the wind shear function cannot be integrated directly into the formula. To overcome this limitation, a simple approach is employed to consider the nonuniformity of the incoming flow using superposition of incoming streamwise velocity and the predicted velocity deficit. By using this assumption, the normalized wake velocity is predicted as follows: where U 0 is the undisturbed incoming velocity at the hub height and U in denotes the incoming velocity at the relative height from the ground either in uniform, powerlaw, or logarithmic profile distribution.

| Model validation
The effectiveness of the proposed model for the prediction of wake velocity profiles for both lateral and vertical directions was validated against CFD simulations and experimental measurements of utility-scale turbines. Table 1 shows the relevant cases used in this study to validate the proposed DG model.

| 2127
Several relevant cases for wake flow measurements behind HAWTs, as tabulated in Table 1, were used in this study for validation purposes. In this study, the utilityscale wind turbines in the typical MW-class turbine sizes range were chosen to represent the actual large turbines.
In Table 1, case 1 corresponds to high-fidelity LES of a utility-scale reference wind turbine from National Renewable Energy Laboratory (NREL), NREL 5 MW reference HAWT. Further details regarding the turbine can be found in Jonkman et al. 48 The LES simulation of full-scale analysis of the turbine was conducted using a noncommercial CFD solver FrontFlow/blue (FFB), 49 and was performed on Japan's Supercomputer Fugaku. The FFB solver uses the finite-element method (FEM) to solve the unsteady incompressible Navier-Stokes (NS) equations numerically under cartesian coordinates. A structured hexahedral mesh was used in the entire computational domain with a total number of 375,960,278 elements. Meanwhile, the blade geometry was directly modeled, thus the detailed vortical structures around the rotating blades were expected to be reasonably produced. By concentrating only on the details of the natural development of the wake flow field induced by the turbine, uniform inlet condition U 0 of 11.4 m/s without turbulence was set at the domain inlet. The mean of fully developed wake flow field data from the last six revolutions was used to validate the effectiveness of the proposed DG model, particularly within the near-wake region where the influence of turbine presence is dominant.
The second case in Table 1 was simulated using SOWFA (Simulator fOr Wind Farm Applications), an LES solver developed by NREL for wind farm simulations. 50 The indirect approach for blade modeling was used by employing the actuator line model (ALM). Since the ALM uses the aerodynamic information of each blade element, thus an accurate representation of the blade aerodynamics can be obtained. An isolated non-yawed utility-scale reference wind turbine INNWIND 10 MW 51 was used in the SOWFA simulation. The computations were carried out under CL-Windcon project, and the LES data are opened publicly as a public database upon request. 52 The averaged data of fully-developed wake flow from the last 600 s of the simulation time (sampling rate of 10 s) were processed as a numerical benchmark. In the simulation, the turbine was operated under the ABL with the reference velocity at the hub height U 0 = 7.87 m/s with the incoming streamwise turbulence intensity TI u of 5.1%. Pairs of streamwise wake velocity profiles in the lateral and vertical directions were sampled at several downstream distances (up to 22D) behind the turbine for validation purposes.
The cases 3-5 in Table 1 correspond to lidar measurements of streamwise velocity along the wake centerline at hub height behind the actual offshore wind turbine AREVA M5000. 53 The turbine is located within Alpha Ventus Offshore wind farm located in the North Sea 45 km north of the island of Borkum, Germany. Based on long-term measurements (2011-2018) from the meteorological mast FINO 1, which is located upstream of the AREVA turbine, the average incoming wind speed at the hub height of 90 m toward the turbine was about 9 m/s. 54 Thus, the lidar measurement (10-min averaging period) from that mean incoming velocity of 9 m/s was chosen. In addition, two resulting wake centerline profiles filtered from lidar data originating from the reference velocities of 7 and 11 m/s were also analyzed. In cases 3-5, the extracted data of streamwise turbulence F I G U R E 2 Case-averaged total RMSE against LES simulations and lidar measurements for pairs of parameters r 0 and coefficient c (Equation 14). LES, large eddy simulation; RMSE, root-mean-square error intensity from FINO1 followed the definition from IEC 61400-3 standard for offshore turbulence intensity, which is based on an approximation of the 90th percentile of the standard deviation of streamwise wind velocity at the hub height. 55 An open-source semi-automated tool WebPlotDigitizer 56 was used to approximate the cited data sets initially published in the image format, resulting in the 90th percentile of streamwise turbulence intensity TI u_90 ≈ 8.2, 7.3, 6.8% 55 and the turbine's Ct ≈ 0.82, 0.81, and 0.79 57 at the mean incoming velocities of 7, 9, and 11 m/s, respectively.

| Model calibration
In this study, the parameters of wake growth rate in both lateral and vertical directions (k y and k z ) in the proposed-DG model were calibrated against the benchmark data sets. The calibration minimized the residuals between the measured and the predicted values by optimizing the adjusted parameters using the MATLAB function fminsearchcon. 61 In addition, the Jensen model, which is popular as an industry-standard wake model, was also evaluated for comparison. Meanwhile, the anisotropic SG model, which considers the anisotropic wake expansion as proposed by Xie et al. 12 was also examined for the same purpose. Here, the formula of initial wake expansion ε = 0.2 β in the original anisotropic SG model led to complex solutions of wake velocity within the near-wake region. Therefore, the initial wake expansion formula was redefined to ε = 0.25 β , which is resulted from analytical derivation, 10 thus resulting in real solutions as expected. The wake expansion parameters of those analytical models were also calibrated to get their optimal performances for the given cases. The tuning results are shown in Table 2, and the details of flow conditions and turbine characteristics of each case are informed in Table 1.
All those analytical models were compared using the benchmark data sets from cases 2-5, which cover the full-wake regions. It should be noted that the calibrated values in Table 2 are not the universal values of the related parameters, thus it needs to be adjusted case by case. This limitation will become a subject of future work to derive the universal expression of the expansion parameters under different flow field conditions.

| Model performance evaluation
To evaluate the effectiveness of the proposed model and the other analytical models quantitatively, the RMSE, the standard deviation of residuals or errors, was employed. The general formula is as follows: where u r exp is the measured velocity, u est is the estimated velocity, and the overbar (¯) means the average of residuals from evaluated data points at the specific direction. In the present study, the directional RMSE, which is based on the evaluated directions at each downstream distance x/D, is calculated as follows: where dir covers lateral, vertical, and streamwise directions. Meanwhile, the average RMSE, which represents the directional RMSE along the evaluated downstream distances within the evaluated region, is expressed as follows: Average RMSE Directional RMSE = .
x D / Finally, the total RMSE, which represents the average RMSE for overall evaluated incoming velocities of the respective cases, is defined as follows: The high-fidelity LES simulation was conducted to provide the details of the wake flow behavior, particularly streamwise velocity distribution within the nearwake region. To focus only on the turbine-induced wake, therefore, mean reference velocity U 0 at the inlet was conditioned uniform without turbulence. Since the NREL 5 MW HAWT model was simulated according to its full-scale size, thus the Reynolds number was on the order of 10 7 , which is typical for MW-class large wind turbines. Pairs of wake profiles in the lateral and vertical directions resulting from the simulation and -their predictions by the proposed DG model are shown in Figures 3 and 4, respectively. By referring to the CFD results as depicted in Figures 3 and 4, the DG shapes was clearly observed within the near-wake region. The appearance of these DG profile shapes has been explained by a very small lift generated around the blade root as the result of its suboptimal shape and the lift is further reduced at the blade tip due to tip vortices. 47,59 This causes a less significant velocity deficit around these regions. Consequently, the wake velocity is still high around the blade root region and almost regained to its undisturbed velocity, particularly around the tip region. Meanwhile, the optimum lift is generated about the mid of the blades, causing a significant velocity deficit around those positions. Thus, minimum lift around the root and tip of the blade and maximum lift about the mid-blade position create two local minima that form the DG velocity profile distribution within the near-wake region.
At the downstream distance just behind the turbine x/D = 1.0, the streamwise wake velocity around the hub center position (y/D = 0) was slightly lower than its surrounding due to the nacelle blockage effect. However, this near-wake appearance was not clearly observed for the vertical direction at the same downstream distance. As the downstream distance increased, the effect of the blockage was slowly faded, thus resulting in a more smooth transition of the wake shape between the hub center and its surroundings. This occurrence created a clear double-Gaussian (DG) shape within the near-wake region. The same behavior was still observed until the domain's outlet. Moreover, it was found that the wake recovery was so slow since the DG shape distribution was still salient at the domain's outlet. This occurrence is reasonable since the onset of the far-wake region is inversely proportional to the incoming streamwise turbulence intensity. Therefore, the transition from the DG to the SG shape profile becomes farther from the turbine, particularly when the incoming turbulence intensity is zero. It also makes the wake growth rates are slow.
From those CFD results, it was confirmed that the proposed DG model showed its effectiveness in predicting streamwise wake profile within the nearwake region. In this case, the best fit of wake growth rates in lateral and vertical directions were about 0.0025 and 0.0024. This means that the wake spread anisotropically.

| Streamwise wake velocity of INNWIND 10 MW HAWT within the full-wake region
Pairs of streamwise wake velocities in the lateral and vertical directions were used to evaluate the effectiveness of the proposed DG model within the full-wake region behind an isolated HAWT. Therefore, the LES simulation data of a utility-scale INNWIND 10 MW HAWT were used. The turbine was operated under ABL with incoming streamwise turbulence intensity ≈5%. In addition, the Jensen and anisotropic SG models were examined for comparison. Those models are known for their robustness for wake velocity estimation within the far-wake region. The predictions of mean wake profiles of INNWIND 10 MW HAWT in the lateral and vertical directions are shown in Figures 5 and 6, respectively.
Within the near-wake region, the proposed DG shape could better reproduce the wake shape distribution in both lateral and vertical directions among the other examined models. However, at the downstream distance up to about x/D = 2, the DG model underpredicted the wake deficit value in both directions. This underprediction was mainly due to the acceleration term in the momentum equation, where the model relied upon, being neglected for the sake of simplicity. Therefore, this underprediction was reasonable since within the distance just behind the turbine, acceleration is more likely to happen because of the significant effect of internal turbulence generated by the turbine. As the downstream distance increased, the internal turbulence effect on the wake flow field gradually decreased, thus allowing the DG model to better estimate the velocity deficit, as it can be seen starting around x/D = 2.5.
The transition from the DG into SG shape was found to occur at about x/D = 4, which represented the onset of the far-wake region. The predictions from the proposed DG and the anisotropic SG models were almost similar at the wake centerline around the hub height. Nevertheless, the Jensen model overpredicted their values up to about x/D = 12. This was mainly due to the isotropic approach being used for the wake expansion in this model, which caused an overprediction of the wake area, thus resulting in an overestimation of the streamwise wake velocity.
Hence, the anisotropic behavior of the wake expansion becomes essential to be considered for a better estimation of the wake velocity, particularly within the far-wake region.
Basically, the proposed DG model and the anisotropic SG model rely upon the same simplification of the momentum equation and by ensuring the F I G U R E 3 Streamwise wake velocity profiles of NREL 5 MW HAWT in the lateral direction within the near-wake region. CFD, computational fluid dynamics; DG, double Gaussian; HAWT, horizontal-axis wind turbine conservation of mass and momentum. This assumption makes both models almost have the exact estimation for the wake profiles within the far-wake region, which can be observed from Figures 5 and 6 at the downstream distance of x/D ≥ 6. Meanwhile, for the downstream distance x/D < 6, the SG model overpredicted the wake shape distributions, azimuthally in the area between the blade mid-span and tip. On the contrary, the model underpredicted their values when approaching the wake centerline, which can be observed up to x/D = 2.5. Afterward, the accuracy of velocity deficit predictions along the centerline improved. These results re-emphasize the major drawback of the anisotropic SG model regarding its lack of accuracy within the near-wake region, where the wake flow is highly affected by the detailed features of the turbine, such as turbine characteristics and geometry. This drawback is reasonable since the SG-based models are based on a pragmatic approach where the wind turbine micrositing is commonly placed inside the far-wake region, thus the nearwake region is not of interest. To address this shortcoming, the proposed DG model considers the influence of the rotor presence within the near-wake region resulting in a better estimation of the formed DG shape and also taking into account the anisotropic behavior of wake expansion, thus increasing its usability under different atmospheric stabilities.

| Wake velocity centerline of utility-scale AREVA 5 MW HAWT within the full-wake region
Inside the wake region, the wake centerline was located approximately downstream of the rotor center at the hub height level. Along this centerline, the maximum velocity deficit occurred, particularly within the near-wake region. The processed lidar data of the mean velocity profile (10 min averaged) along the wake centerline behind utility-scale AREVA M5000 HAWT were used to investigate this behavior further. In this study, the lidar data were utilized as an experimental benchmark to evaluate the effectiveness of the proposed and the other examined models for the wake centerline prediction. Since the analyzed wake data were only a function of the downstream distances x/D, one-dimensional fitting was sufficient for each analytical model to estimate the wake centerline profiles at the hub height level. The measurement results and their estimations from each analytical model are shown in Figure 7.  Figure 7, the wake velocity centerline profiles had the maximum velocity deficit of about x/D ≈ 2 downstream of the turbine. However, their magnitudes increased slightly when the inflow velocity was higher. Notable results were also found regarding the wake expansion under different incoming streamwise turbulence intensities. For the analyzed lidar cases, as the inflow velocity increased, the incoming turbulence and the wake growth rate became weaker. This behavior can be seen from the centerline profiles, where the wake   recovered faster at the higher turbulence intensity. In this case, the fastest recovery happened at the incoming velocity of 7 m/s. The dependence of turbulence intensity on the wake recovery rate can be represented by the onset of the far-wake region as formulated in Equation (14). Since the onset of the far-wake is conversely proportional to the incoming turbulence, therefore for higher turbulence intensity, the wake recovers faster. This behavior has been confirmed by the analyzed utilityscale wake measurements, as shown in Figure 7 where the onset of the far-wake region became closer to the turbine as the turbulence rose. This occurrence is marked by black-dashed vertical lines in the figure.
In general, the wake growth rates from all examined analytical models show the same tendency, where their tuned values increased with the higher turbulence intensity. This means that the turbulence intensity is worth to be considered for analytical wake modeling, particularly to get a better understanding of the wake recovery rate. From Figure 7, it was observed that within the near-wake region, particularly from the downstream distance just behind the turbine to the distance where the maximum velocity deficit occurred (x/D ≈ 2), the DG model provided better estimations of the wake centerline profiles compared to the other analytical models. The onset of the far-wake region estimated by Equation (14) in the analyzed cases is at 3.0 < x/D < 3.6, which is still in the range for the far-wake onset ranging from 2-4D downstream of the turbine. 38 Around the onset of the farwake region, the DG model outperformed the other models for the wake centerline predictions. However, its predictions slightly deviated observably starting from x/D ≈ 6, where the model underestimated the centerline profiles. Afterward, the accuracy improved from x/D ≈ 11 for U 0 = 7 m/s and x/D ≈ 13 for U 0 of 9 and 11 m/s. Meanwhile, around the onset of the far-wake region, F I G U R E 7 The wake centerline profiles downstream of the AREVA 5 MW HAWT at different incoming streamwise velocities. DG, double Gaussian; HAWT, horizontal-axis wind turbine; SG, single Gaussian both Jensen and Anisotropic SG models significantly overestimated the wake centerlines from the overall reference velocities. Subsequently, the accuracy of these models was found improved started from x/D ≈ 6 for U 0 = 7 m/s and x/D ≈ 7 for U 0 of 9 and 11 m/s. Yet, the same trend of deviation for wake centerline predictions within the far-wake was also observed in the Jensen and Anisotropic SG models, specifically from x/D ≈ 6 for U 0 of 7 and 11 m/s and x/D ≈ 7 for U 0 of 9 m/s up to x/D ≈ 12.5 for the overall reference velocities.
Qualitatively, the evaluation of the analyzed wake models based on the available lidar data revealed that the DG-based model could give a better estimation of the wake centerline profiles within the near-wake to the farwake region about x/D ≈ 6. Surpassing that downstream position, all the evaluated models offered a feasible prediction with relatively trivial errors against the measurement data. These errors might also be caused by the presence of downwind turbines or overlapping wakes at the far-field, thus leading to a deviation from the typical wake flow field behind an isolated HAWT. The possibility of these trivial errors made the prediction accuracies from the examined analytical models within the far-wake region were more acceptable.

| Statistical evaluation
This evaluation aims to quantify the residuals between the proposed model and the other examined models against the benchmark wake data sets from the utilityscale turbines using the RMSE statistical model. The previously evaluated cases from LES simulations and lidar measurements were used to further analyze the performance of the proposed DG and the other analytical models quantitatively. For the CFD results, only the wake data sets from SOWFA simulation for wake simulation behind an isolated INWIND 10 MW HAWT were evaluated further since it covers a full-wake region. Perfect wake model performance would give an RMSE value of 0.
In this study, the RMSE analyses from the examined wake models were categorized into four main regions: near-wake, far-wake, full-wake, and practical regions. Here, the near-wake region refers to the downstream distance between the turbine and the onset of the farwake region, which is approximated using Equation (14). Meanwhile, the far-wake region is started from the onset of the far-wake region and beyond. Next, the full-wake region comprises the overall wake region ranging from the near-wake to the far-wake region. The last region is the practical region, which is based on the actual interrow spacing between wind turbines in operational wind farms. This region ranges from x/D = 2.4, used as the interrow spacing between the installed turbines in 20 MW Middelgrunden offshore wind farm, 58 up to x/D = 11.0, which is around the maximum interrow spacing between the turbines in specific wind direction (10.4 D at θ wind = 312) within Horns Rev offshore wind farm. 60 It is expected that the set range for the practical region (2.4 ≤ x/D ≤ 11) could cover all of the actual wind turbine micrositings inside the operational wind farms nowadays thus the effectiveness of the examined wake models can be evaluated comprehensively.
The first case to be evaluated is from lidar measurements for wake centerline velocity profiles behind utilityscale offshore wind turbine AREVA 5 MW HAWT. The performance of each model was examined into four main categories within the wake region (up to x/D = 15) resulting from different incoming velocities. In the present case, the onsets of the far-wake region according to Equation (14) are ≈3, 3.3, and 3.56 D for the incoming velocities of 7, 9, and 11 m/s, respectively. The evaluation results in terms of average RMSE for the near-wake, farwake, and full-wake regions at the incoming velocities of 7, 9, and 11 m/s are shown in Figures 8-10, respectively.
From the RMSE results, the proposed DG model gave the most accurate prediction for the wake centerline profiles among the other examined models within the near-wake region. This is mainly due to the model's ability to reproduce the velocity degradation of the wake centerline profiles. Meanwhile, It is interesting to note that the top-hat shape from the Jensen model outperformed the centerline prediction by the anisotropic SG model at each of the incoming velocities in the nearwake region. In contrast, the anisotropic SG model showed its robustness within the far-wake region by outperforming the prediction from the Jensen model at each incoming velocity and only at U 0 of 11 m/s for the proposed DG model. On average, the RMSE results within the far-wake region showed comparable performances among the evaluated models, except for U 0 of 9 m/s, where the proposed DG was much more precise than the others. For the full-wake analyses, the proposed DG provided the best estimation of the wake centerlines, then followed by the anisotropic SG model and the Jensen model afterward. Finally, the total RMSE of each examined model is plotted in Figure 11.
Total RMSE in Figure 11 was calculated based on the mean deviation for the wake centerline predictions under overall evaluated incoming velocities. The results revealed that the proposed DG had the best predictions within the near-wake, far-wake, and full-wake regions. The most notable performance from the proposed DG was its prediction accuracy within the near-wake region, which was significantly better than the other examined SOESANTO ET AL.
| 2137 F I G U R E 8 Average RMSE against lidar measurement (case 3) within the near-wake, far-wake, and full-wake regions at the incoming reference velocity of 7 m/s. DG, double Gaussian; RMSE, root-mean-square error; SG, single Gaussian F I G U R E 9 Average RMSE against lidar measurement (case 4) within the near-wake, far-wake, and full-wake regions at the incoming reference velocity of 9 m/s. DG, double Gaussian; RMSE, root-mean-square error; SG, single Gaussian F I G U R E 10 Average RMSE against lidar measurement (case 5) within the near-wake, far-wake, and full-wake regions at the incoming reference velocity of 11 m/s. DG, double Gaussian; RMSE, root-mean-square error; SG, single Gaussian F I G U R E 11 Total RMSE against lidar measurements (cases 3-5) within the nearwake, far-wake, and full-wake regions. DG, double Gaussian; RMSE, root-mean-square error; SG, single Gaussian models. This evaluation result also reaffirms the limitation of the Jensen and anisotropic SG model, which are unsuitable for the wake centerline estimation within the near-wake region. It is also interesting to note that the proposed DG model had the best accuracy among the other models, even in the far-wake region. This was plausible due to the DG model's capability to roughly model the wake centerline transition between the nearwake and the far-wake regions. As a result, its accuracy around the onset of the far-wake became better than the other models, thus enhancing its overall performance within the far-wake region and likewise in the full-wake region. It was observed that within the full-wake region, the Jensen model was least accurate among the other models. However, its prediction was not so different from a more complex anisotropic SG model and still relatively good for a simple top-hat shape approach being used.
The last category is the practical region, which covers the minimum to maximum interrow spacing between the actual wind turbines in operational wind farms. The average RMSE evaluation for each evaluated turbine at each incoming velocity within this region is shown in Figure 12. Meanwhile, the total RMSE resulting from each model at the overall evaluated incoming velocities is shown in Figure 13. The proposed DG model had better performance than the other examined models for the given incoming velocities and followed by the anisotropic SG model and the Jensen model, respectively. Averaged from the practical region, the proposed model gave much better performance over the other models as expected. This was mainly due to its robust performance, particularly around the onset of the far-wake region, which significantly strengthens its overall performance. Otherwise, the inability of both Jensen and Anisotropic DG models to reproduce the centerline profiles around the onset of the far wake region made those models work effectively after x/D ≈ 6 based on the analyzed lidar cases. The trend was the same, where both models significantly overestimated the wake centerline profiles at each evaluated incoming velocity around the onset of the far-wake region. Since the available wind energy is the cube of wind velocity, thus this significant overestimation may bring a significant error in predicting the energy output, particularly in tightly packed wind farms.
It should be noted that the lidar data were obtained from the wake centerline measurements in the streamwise direction to represent the wake velocity recovery. Therefore, for the sake of completeness, further wake evaluations that consist of RMSE analyses of the wake expansion in both lateral and vertical directions along a full-wake region are still needed to comprehensively examine the evaluated model performances. Therefore, pairs of directional streamwise wake expansions (in the lateral and vertical directions) resulting from CFD simulation behind an isolated utility-scale INWIND 10 MW HAWT were used to achieve this purpose. The wind turbine was simulated under ABL with incoming turbulence. The onset of the far-wake region was located at x/D = 4.41, according to Equation (14). The wakes were numerically measured along the downstream distances ranging from 1 ≤ x/D ≤ 22 on cross-sectional planes, which centralized from the hub center toward the blade span of both sides (+ and −) of each of the lateral and vertical directions. The evaluation results of directional RMSE for the wake profile predictions are shown in Figure 14. Figure 14A,B show the directional deviations of the examined wake models against the CFD data for the F I G U R E 12 Average RMSE against lidar measurements (cases 3-5) within the practical region at different incoming streamwise velocities. DG, double Gaussian; RMSE, root-mean-square error; SG, single Gaussian F I G U R E 13 Total RMSE against lidar measurements (cases 3-5) within the practical downstream distances. DG, double Gaussian; RMSE, root-mean-square error; SG, single Gaussian wake profile predictions in the lateral and vertical directions, respectively. Meanwhile, Figure 14C displays the average deviations in both directions along the evaluated downstream distances. Starting with the simplest wake model, the top-hat shape approach from the Jensen model had the least accurate prediction of the wake profiles among the other models. It was observed that the residual gaps between this top-hat model and the Gaussian-based models were significant until x/D = 18, where the largest gap area was found in the lateral direction. This occurred mainly due to its top-hat simplification for the wake distribution, which in reality tends to form a Gaussian shape. Therefore, if the wake velocity has yet to recover significantly, the overall disparity of the wake distribution toward the blade span direction will remain prominent. However, its prediction for the streamwise wake centerline may improve faster after the onset of the far-wake region as in the previous lidar cases.
On average, the proposed DG performance within the near-wake region was better than the other models. Its prediction almost coincided with the anisotropic SG model at further downstream distances but with lower residuals, observably starting from x/D = 6. These statistical results also pointed out the effectiveness of the proposed model within the far-wake region. It is also interesting to note that the anisotropic SG model showed its robustness from x/D = 6.0, which is consistent with the analyzed lidar cases. This observation emphasizes the SG model compatibility for wake prediction that started from the downstream distance of 6D behind the turbine. Furthermore, its accuracy, especially the Jensen model, was not so good around the onset of the farwake region either in the present CFD case or the previous lidar cases. This drawback could become another advantage of the proposed DG model among the other examined models for its effectiveness around the onset of the far-wake region.
The average RMSE evaluation of the wake model performance based on the wake regional categories is shown in Figure 15. Within the near-wake region, the proposed DG had the best accuracy for wake profile estimations in the lateral and vertical directions among the other models, followed by Jensen and anisotropic SG models, respectively. Meanwhile, the anisotropic SG model outperformed the Jensen model performance in the far-wake region. There was less difference in the overall deviations between the proposed DG and anisotropic SG models within the far-wake region, mainly due to the similarity of the Gaussian approach being used, particularly when the DG shape profile has transformed into the fully SG shape profile. Nonetheless, the proposed DG model still provided the best accuracy among the other models in this region. Next, the RMSE results for wake profiles within the full-wake region showed that the proposed DG had the best accuracy along the evaluated wake region behind the turbine and was followed by the anisotropic SG and the Jensen models, respectively. In addition to the DG model's ability to model a single wake shape within the far-wake region, the model capability to predict the wake shape transition, which tends to form a flatlike shape from around the blade mid-span to the wake center, becomes its another advantage that might be difficult to reproduce by the other examined models.
Finally, the average RMSE within the practical region is presented in Figure 16. The average deviation along the evaluated region showed that the proposed DG model had the best fit against the benchmark CFD data sets. It was observed that the inaccuracy of the Jensen model prediction within the practical region was inevitable due to the top-hat shape assumption being used, which is not consistent with the actual wake profile shape. Meanwhile, the performance of the SG model in this region was much better than the Jensen model. However, the SG model's significant deviations from the beginning of the practical region to around the onset of the far-wake negatively affected its whole performance in the practical region. Thus, the proposed DG model still outperformed the SG model performance. As mentioned previously, the proposed DG model has an advantage regarding its higher accuracy among the other examined models, particularly within the near-wake region to around the onset of the far-wake region, which is beyond the scope of the conventional approach for analytical wake modeling. Hence, the proposed model could be employed to fill this gap by offering a better estimation for the wake prediction, without exception around the onset of the farwake region.

| CONCLUSION
A DG model for analytical wake modeling is proposed to complement the previously developed DG model for predicting the streamwise velocity distribution behind an isolated non-yawed HAWT. The anisotropic nature of the wake growth rate downstream of the wind turbine is considered, thus facilitating the model's usability under different atmospheric stabilities. Moreover, constant values of the radial position of Gaussian minima r 0 and scale factor c in Equation (14), which were determined empirically, are suggested for the DG-based wake modeling behind utility-scale HAWTs. Finally, the effectiveness of the proposed model was validated against F I G U R E 15 Average RMSE against LES simulation of INWIND 10 MW (case 2) within the near-wake, far-wake, and full-wake regions. DG, double Gaussian; LES, large eddy simulation; RMSE, rootmean-square error; SG, single Gaussian F I G U R E 16 Average RMSE against LES simulation of INWIND 10 MW (case 2) within the practical downstream distances. DG, double Gaussian; LES, large eddy simulation; RMSE, root-mean-square error; SG, single Gaussian the LES simulations and lidar measurements of the wake flow field behind utility-scale wind turbines.
The validation results showed that the proposed model could well reproduce the DG shape of the wake profiles within the near-wake region. Meanwhile, a SG-like wake shape distribution within the far-wake region could also be well estimated by the proposed model, thus confirming its compatibility within the full-wake region. Furthermore, a comprehensive evaluation was performed to examine the model performance under different inflow conditions and wake regions using RMSE statistical analysis. In general, the proposed DG model provided a feasible prediction within all wake region categories and outperformed the other examined wake models. It was found that the model has an advantage regarding its higher accuracy within the nearwake region to around the onset of the far-wake region, which is beyond the scope of the conventional approach for analytical wake modeling. Thus, the proposed model may have the potential to give a better prediction for the wake flow field within tightly packed wind farms with relatively low computational cost.
In the proposed model, the wake growth rate parameters are still needed to be tuned case by case. This limitation will be considered in future work to determine the universal expression for the parameters based on the different flow field conditions. In addition, some wake flow fields resulting from high-fidelity LES simulations will be employed to achieve this goal.

ACKNOWLEDGMENTS
This study was supported by Ministry of Education, Culture, Sports, Science and Technology (MEXT) Japan as "Program for Promoting Research on the Supercomputer Fugaku" (Digital Twins of Real World's Clean Energy Systems with Integrated Utilization of Super-simulation and AI) (Grant/ Award Number: hp210175).

CONFLICTS OF INTEREST
The authors declare no conflicts of interest.

DATA AVAILABILITY STATEMENT
The data presented in this article are available from the corresponding author upon request.