text:
\section{Introduction}\label{sec:Intro}

A deeper understanding of the physical processes governing the operation of modern wind farms is instrumental to enable reduction of the levelized cost of energy (LCOE) and increase of wind energy penetration in the global power system \cite{SanzRodrigo2017}. Specifically, many mechanisms governing the aerodynamics of wind turbines immersed in the atmospheric turbulent boundary layer remain still elusive nowadays \cite{Veers2019}. By extracting kinetic energy from the incoming flow, wind turbines generate wakes, which can persist for several kilometers downstream \cite{Christiansen2006,Platis2018,Lundquist2019}, and cause power losses and enhanced fatigue loads for downstream wind turbines \cite{Porte-agel2019}. Wind industry stakeholders identify turbine wakes as one of the major causes for energy losses for onshore power plants in North America \cite{Bailey2014}. Quantitative studies of wake effects based on the analysis of supervisory control and data acquisition (SCADA) data showed energy losses between 2-4\% for onshore sites \cite{El-Asha2017} and 10\% to 20\% of the annual energy production (AEP) for offshore wind farms \cite{Barthelmie2010}. Nonetheless, a recent analysis of data covering ten years of operations for the Lillgrund offshore wind farm revealed a reduction in power capture as high as 28\% of the nominal capacity \cite{Sebastiani2020}.

The prediction of wake impact on wind farm performance is complicated by the influence of the atmospheric stability on wake recovery \cite{Magnusson1994,Iungo2014,CarbajoFuertes2018,Zhan2020}, which results in significantly reduced power losses for high turbulence intensity, convective atmospheric conditions compared to low turbulence intensity, stable atmospheric conditions \cite{BarthelmieWE2010,Hansen2012,El-Asha2017}. An additional hurdle for predicting wake-induced power losses is represented by the coalescence of multiple wakes, which creates a spatially heterogeneous and highly turbulent flow within the wind farm boundary layer. Furthermore, the pressure field created by the thrust of the rotors induces blockage \cite{MeyerForsting2017,Wu2017,Ebenhoch2017,Bleeg2018,Branlard2020,Schneemann2021} and speedups \cite{Nishino2012,McTavish2014,MeyerForsting2016}, which can significantly impact wind turbine performance, yet generally not modeled in engineering approaches due to the excessive computational burden that a coupled solution of the continuity and momentum equation entails when including the pressure field.

During the past four decades, several wake and wind farm models have been proposed, spanning a breadth of approaches and levels of sophistication. Wake models can be broadly classified into two categories: analytical models and computational fluid mechanics (CFD) methods \cite{Porte-agel2019}. The simplest analytical models assume a top-hat velocity profile and calculate wake width and velocity deficit based on mass \cite{Jensen1983} or momentum \cite{Frandsen2006} conservation. Other analytical wake models assume self-similarity of the wake velocity profiles, in analogy to the shear-flow theory \cite{Schlichting1979,Abramovich1963}, to predict maximum velocity deficit and wake width as a function of the downstream distance from the turbine rotor \cite{Lissaman1979,Voutsinas1990,Larsen1988}. The self-similarity assumption is also at the basis of the Gaussian wake model \cite{Bastankhah2014} and the subsequent evolutions, such as the Ishihara model \cite{Ishihara2018}, the elliptic-Gaussian model \cite{Xie2015}, the Gaussian plume dispersion model \cite{Cheng2018}, the double-Gaussian model \cite{Schreiber2020}, and the super-Gaussian models \cite{Shapiro2019,Blondel2020}. It is noteworthy that the self-similarity hypothesis is strictly valid only if the wake velocity deficit is significantly smaller than the freestream velocity, which is typically true only in the very far-wake \citep{Pope2000}.

The extremely low computational cost and the simple implementation contributed to the popularity of engineering wake models, which are still included in several state-of-the-art tools for wind farm design and modeling \cite{Crasto2012,Mortensen2016,Gebraad2016}. The main drawbacks of most analytical wake models are the need for tuning parameters regulating the wake recovery and the shape of the velocity profiles \citep{Zhan2020WES}. Furthermore, engineering wake models typically require the use of superposition principles in the case of multiple wakes, which introduce a further inaccuracy \cite{Gunn2016,Zong2020}. Nevertheless, if properly calibrated, these models exhibit a mean absolute error in the prediction of wake power losses smaller than 15\% and, unexpectedly, comparable to that of more complex CFD models \cite{Sebastiani2020,Barthelmie2006,Barthelmie2009,Moriarty2014,Archer2018}.

CFD models, especially Reynolds-Averaged Navier-Stokes (RANS) and Large Eddy Simulations (LES) have been increasingly employed to simulate wind farm flows by solving the governing equations for turbulent flows \cite{Santoni2018,Santoni2020}. Specifically, LES have provided insightful information on the turbine wake flow physics (see reviews by \cite{Mehta2014,Breton2017}); however, their high computational cost ($\sim 10^3-10^4$ CPU hours \cite{Porte-agel2019}, depending on the Reynolds number\cite{Sanderse2011}) represent a serious hindrance for their utilization in computationally intensive tasks, such as online control, real-time diagnostic applications, and layout optimization. Furthermore, the sub-grid model needs to be carefully selected to avoid inconsistent results near the ground \cite{Bou-Zeid2005}, for stratified flows \cite{Wan2011} and in the proximity of the turbines \cite{Porte-agel2019}.

Although very sensitive to the type of turbulent closure adopted \cite{Rethore2009,Breton2017}, RANS models have computational costs generally three orders of magnitude smaller than LES \cite{VanDerLaan2015}. The time required for a RANS simulation is intrinsically dependent on the simplifying assumption applied to the Navier-Stokes equations and the solver implemented. Linearization of the advection term is a commonly adopted approximation for early RANS models \cite{Taylor1980,Sforza1981,Crespo1985} and received renewed interest more recently thanks to the achieved low computational costs. Modern linearized models include FUGA \cite{Ott2011}, ORFEUS \cite{Ebenhoch2017}, the curled wake model \cite{Martinez-Tossas2019}. Another popular class of models is the parabolic RANS, where, by leveraging to the boundary-layer approximation, computational costs are reduced by solving the RANS equations from the inlet moving in the downstream direction by neglecting the elliptic nature of subsonic flows. The classic models of Liu \cite{Liu1983}, Ainslie \cite{Ainslie1988}, and UPMWAKE \cite{Crespo1985} belong to this category. More recently, a parabolic axisymmetric solver for turbine arrays was optimally calibrated based on high-fidelity simulations \cite{Iungo2015,Iungo2017,Santhanagopalan2018}. The commercial tool WakeBlaster \cite{Bradstock2020} is also a parabolic 3D wind farm model that can solve the flow of a medium-size wind farm in a few seconds. Furthermore, the curled wake model in his fully linearized version \cite{Martinez-Tossas2019} and the newest release \cite{Martinez-Tossas2020} use a parabolic solution to reduce the computational time. The main drawback of parabolic models is the breakdown of the boundary layer approximation in the near wake \cite{Pope2000} and the inability to simulate the effects of pressure (namely blockage and speedups) on the wind field \cite{Bleeg2018}.

A more detailed flow characterization can be achieved through elliptic RANS models, which solve the full set of RANS equations at the cost of a significantly higher computational requirement. Crespo and Hend\'andez\cite{Crespo1991} formulated an elliptic implementation of UPMWAKE, whose results showed only slight differences with respect to the original parabolic version. Masson et al.\cite{Masson1997} formulated a steady laminar axisymmetric RANS solver based on the actuator disk theory, subsequently enhanced through the inclusion of a turbulence model\cite{Leclerc1999} and a 3D formulation \cite{Ammara2002}.
In an attempt to reduce the dimensionality and, thus, computational burden, several authors resolved the wind farm flow at hub height in a purely 2D fashion \cite{Soleimanzadeh2014,Annoni2015,King2016,Adcock2018}, which however can result in an excessive curvature of the streamlines due to the unphysical vertical confinement \cite{King2017}.

A crucial aspect of a RANS model is the modeling of the turbulent Reynolds stresses, aimed at capturing the complex role of atmospheric turbulence, blade- and wake-generated turbulence, and wake dynamics, such as wake meandering \cite{Larsen2008,Keck2012}. The vast majority of the RANS models adopt the linear turbulent eddy viscosity hypothesis \cite{Boussinesq1897}, with just a few examples of Reynolds Stress transport \cite{Gomez-Elvira2005,Cabezon2011} and non-linear eddy viscosity closures \cite{VanderLaan2013}. It is worth mentioning that the eddy viscosity approach has been disproved based on theoretical arguments and experimental evidence for flows undergoing rapid distortion \cite{Pope2000}, and particularly for the flow near wind turbines \cite{Rethore2009}. Nevertheless, the eddy viscosity can provide a useful characterization of mean quantities for simple shear flows \cite{Pope2000}.

Algebraic closures for turbine wakes specify the value of the turbulent eddy viscosity or mixing length a priori \cite{King2017}, based on the Monin-Obukhov similarity theory \citep{Liu1983,Ammara2002,Ott2011,Ebenhoch2017,Martinez-Tossas2019}, wake characteristics \cite{Sforza1981} or both \cite{Ainslie1988,Lange2003,Bradstock2020}. A fairly complex algebraic closure is implemented in the dynamic wake meandering and includes effects of atmospheric stability \cite{Keck2014a}, shear, and added turbulence \cite{Keck2015}. Some authors pursued a more data-driven approach, by calibrating the turbulence model against LES \cite{Iungo2017} or experimental data collected with light detection and ranging (LiDAR) systems \cite{Iungo2018,Adcock2018}. The standard two-equation models (namely $\kappa-\epsilon$ and $\kappa-\omega$), despite their well-documented accuracy in the solution of other classic fluid mechanics problems, are known to be over-diffusive in the near wake \cite{Cabezon2011,Prospathopoulos2011}, where the rapid flow changes, high velocity, and pressure gradients undermine the fundamental assumptions of the eddy viscosity model \cite{Rethore2009}. This has spurred significant research work in the development of advanced turbulence closure for turbine wake flow \cite{ElKasmi2008,VanDerLaan2015,VanDerLaan2015b,Hennen2017}.

The foregoing appraisal of the capabilities and limitations of the state-of-the-art turbine wake models inspired the formulation of a novel mid-fidelity model for wind farm flows: the Pseudo-2D RANS (P2D-RANS hereafter). This tool is conceived to provide a comprehensive description of the time-averaged wind farm aerodynamics, including multiple wake interactions and the effects of atmospheric stability and pressure field while keeping the computational costs low enough to simulate a medium-size farm in a few CPU seconds. With these aims, three main novelties are introduced: $i)$ the 3D flow equations are reduced to a 2D mathematical problem by adopting the shallow-water model with tailor-made corrections for near-wake vertical fluxes and dispersive terms; $ii)$ the partially-parabolized Navier Stokes (PPNS) equations are solved by retaining the full elliptic pressure field yet ensuring numerical efficiency of the iterative marching scheme of the numerical solver; $iii)$ the turbulence closure and actuator disc model are optimally calibrated based on LiDAR data clustered based on incoming wind speed and turbulence intensity \citep{Zhan2020,Zhan2020WES}. We will show in the following that the P2D-RANS is also capable of accurately simulating highly challenging wind conditions where the pressure field induced by the turbine rotors can modify the incoming wind leading to local speedups, and, in turn, increased power capture, which has been confirmed experimentally. Accuracy of the P2D-RANS is assessed for the case study of an onshore wind farm located in North Texas including 25 2.3 MW turbines. The current implementation of the code can simulate the steady-state operation of the investigated wind farm in around 80 s on an i5 single-core laptop computer.

The remainder of this manuscript is organized as follows: the first two sections describe the conceptual (Sec. \ref{sec:Conceptual}) and numerical (Sec. \ref{sec:Computational}) modeling of the depth-averaged, partially-parabolized Navier-Stokes equations. Sec.  \ref{sec:Calibration} provides a description of the experimental site and the dataset used for the tuning of the LiDAR-driven actuator disk and turbulence closure, as well as a detailed outline of the calibration procedure. The novel approach is then verified (Sec. \ref{sec:Verification}), validated against real power data, and compared to other models of wind turbine wakes (Sec.  \ref{sec:Validation}). Finally, concluding remarks are discussed in Sec. \ref{sec:Conclusions}.

\section{Conceptual model}\label{sec:Conceptual}
Under steady conditions, the flow around wind turbines can be described through the RANS equations for incompressible turbulent flows:
\begin{equation}\label{eq:RANS}
\begin{cases}
\frac{\partial \overline{u}_j}{\partial x_j}=0\\
\overline{u}_j\frac{\partial \overline{u}_i}{\partial x_j}=-\frac{\partial \overline{p}}{\partial x_i}+\nu  \frac{\partial^2 \overline{u}_i}{\partial x_j \partial x_j}- \frac{\partial \overline{u'_i u'_j}}{\partial x_j} + \overline{f}_i,
\end{cases}
\end{equation}
where $\overline{\mathbf{u}}$, $\overline{p}$, and $\overline{\mathbf{f}}$ are the Reynolds-averaged velocity vector, pressure, and body force vector, respectively ($\overline{f}_1$ representing the turbine rotor thrust), and $\nu$ is the molecular air viscosity. Within the present notation, repeated indices imply summation and bold quantities indicate vectorial quantities. The $x_1$-axis (or $x$) is aligned with the turbine axis and points downstream, $x_3$ (or $z$) is the vertical direction, and $x_2$ (or $y$) is the traversal direction defined according to a right-handed reference system (Fig. \ref{fig:Turbine_sketch}). The velocity components in the $(x,~y,~z)$ directions are $(u,~v,~w)$, respectively. Throughout the paper, all the physical quantities have been made non-dimensional through the turbine diameter, $D$ (length), incoming wind speed, $U_\infty$ (velocities), undisturbed dynamic pressure $\rho U_\infty^2$, (pressure), and the purely dimensional terms $U_\infty D$ (viscosity) and $\rho U_\infty^2/D$ (forces per unit volume).
\begin{figure}[b!]
\centerline{\includegraphics[width=0.65\textwidth]{Turbine_sketch.pdf}}
\caption{Schematic of the reference system of the P2D-RANS and the shallow-water decomposition of velocity.\label{fig:Turbine_sketch}}
\end{figure}

The Boussinesq hypothesis allows to express the deviatoric part of the Reynolds stress tensor to be proportional to the strain rate tensor through the turbulent eddy-viscosity, $\nu_T$, as follows\cite{Boussinesq1897}:
\begin{equation}\label{eq:EV}
    \overline{u'_i u'_j}= \frac{1}{3}\delta_{ij} \overline{u'_i u'_i}
    -\nu_T \left (\frac{\partial \overline{u}_j}{\partial x_i}+\frac{\partial \overline{u}_i}{\partial x_j} \right),
\end{equation}
where $\delta_{ij}$ is the Kronecker delta. By incorporating the first term on the RHS of Eq. (\ref{eq:EV}) into a modified pressure and neglecting the molecular viscosity, we obtain:
\begin{equation}\label{eq:RANS_EV}
\begin{cases}
\frac{\partial \overline{u}_j}{\partial x_j}=0\\
\overline{u}_j\frac{\partial \overline{u}_i}{\partial x_j}=-\frac{\partial \overline{p}}{\partial x_i}+\nu_T  \frac{\partial^2 \overline{u}_i}{\partial x_j \partial x_j} + \overline{f}_i.
\end{cases}
\end{equation}

For the P2D-RANS, the dimensionality of the problem is reduced from 3D to 2D by adopting a depth-averaged (or shallow water) approximation  of the RANS equations, which requires the horizontal length and velocity scales of the mean flow to be significantly greater than their respective vertical quantities \citep{Cea2007}. For wind farm flows, this assumption may not be satisfied in presence of significant vertical shear (which is at the moment not considered and can eventually be added a posteriori through a superposition model \cite{Leclerc1999,Larsen2008}) or complex terrain, whose modeling is beyond the scope of this work.

We define a generic Reynolds-averaged scalar, $\overline{\xi}$, as the sum of a vertically-averaged value, $\hat{\xi}$, and a $z$-varying (or dispersive) component, $\tilde{\xi}$, as sketched in Fig. \ref{fig:Turbine_sketch}. The depth-averaged parameter is defined as:
\begin{equation}\label{eq:DepthAverage}
    \hat{\xi}(x,y)=
\frac{1}{z_2-z_1}\int_{z_1}^{z_2} \overline{\xi}(x,y,z) dz
\end{equation}
where $z_1$ and $z_2$ are the vertical coordinates of the bottom and top boundaries, respectively, of the integration domain. In this work, these boundaries are selected as $z_2=-z_1= D/2$ to include the entire rotor layer, which is convenient to study wind turbines with uniform hub height and rotor diameter on flat terrain.

Application of the operator defined in Eq. (\ref{eq:DepthAverage}) to both sides of Eqs. (\ref{eq:RANS_EV}) in conservative form, yields:
\begin{equation}\label{eq:RANS_SW}
\begin{cases}
\frac{\partial \hat{u}}{\partial x}+\frac{\partial \hat{v}}{\partial y}+\frac{\overline{w}_2-\overline{w}_1}{z_2-z_1} =0
\\
\frac{\partial \hat{u}} {\partial x}\hat{u}+\frac{\partial \hat{u}}{\partial y}\hat{v}+\frac{\overline{u}_2\overline{w}_2-\overline{u}_1\overline{w}_1-\hat{u}(\overline{w}_2-\overline{w}_1)}{z_2-z_1}
=-\frac{\partial \hat{p}}{\partial x}+
\nu_T \left( \frac{\partial^2{\hat{u}}}{\partial x^2}
+\frac{\partial^2{\hat{u}}}{\partial y^2}\right)\\
+\hat{f_x} + \frac{\nu_T}{z_2-z_1} \left(\frac{\partial{\overline{u}}}{\partial z}\bigg|_2-\frac{\partial{\overline{u}}}{\partial z}\bigg|_1 \right)
-\frac{\partial\widehat{\tilde{u}\tilde{u}}}{\partial x}-\frac{\partial \widehat{\tilde{u}\tilde{v}}}{ \partial y}
\\
\frac{\partial \hat{v}} {\partial x}\hat{u}+\frac{\partial \hat{v}}{\partial y}\hat{v}+\frac{\overline{v}_2\overline{w}_2-\overline{v}_1\overline{w}_1-\hat{v}(\overline{w}_2-\overline{w}_1)}{z_2-z_1}
=-\frac{\partial \hat{p}}{\partial y}+
\nu_T \left( \frac{\partial^2{\hat{v}}}{\partial x^2}+ \frac{\partial^2{\hat{v}}}{\partial y^2}\right)\\
+
\hat{f_y} +
\frac{\nu_T}{z_2-z_1}  \left(\frac{\partial{\overline{v}}}{\partial z}\bigg|_2-\frac{\partial{\overline{v}}}{\partial z}\bigg|_1 \right)
-\frac{\partial\widehat{\tilde{u}\tilde{v}}}{\partial x}-\frac{\partial \widehat{\tilde{v}\tilde{v}}}{ \partial y},\\
\end{cases}
\end{equation}
where the subscript 1 and 2 refer to 3D quantities evaluated at $z_1$ and $z_2$, respectively. Eqs. (\ref{eq:RANS_SW}) include depth-averaged ($\hat{.}$), 3D ($\overline{.}$), and dispersive ($\tilde{.}$) terms. While the first are explicitly solved by the P2D-RANS model, the last two contributions need to be modeled as a function of the depth-averaged field. By neglecting the wake swirling motion, which is only present in the near wake \cite{Iungo2013} and has negligible effects on the streamwise velocity and power \cite{Meyers2010}, and assuming an axisymmetric wake velocity field, which can be recovered for cases with no yaw misalignment of the turbine rotor by subtracting the incoming vertical profile of the incoming streamwise velocity \cite{Medici2006,Chamorro2009}, it is possible to calculate these fluxes following the approach proposed by Boersma et al. \cite{Boersma2018}.
Specifically, this allows defining linear operators that map the depth-averaged field $(\hat{u},\hat{v})$ onto the corresponding axisymmetric streamwise velocity field $(\overline{u}_x,\overline{u}_r)$, as follows (see Appendix I for details):
\begin{equation}
\label{eq:3D_corr}
\begin{cases}
\overline{u}_x=M_u^{-1}(\hat{u})\\
\overline{u}_r=M_v^{-1}(\hat{v}).\\
\end{cases}
\end{equation}
Once the axisymmetric velocity components corresponding to a given provisional depth-averaged field are known, the vertical fluxes and the dispersive stresses in Eqs. (\ref{eq:RANS_SW}) can be numerically evaluated and fed as source terms into the model in an iterative way, as it will be detailed in Sec. \ref{sec:Computational}. In the case of side-wise multiple wake interaction, which causes a breakdown of the assumption of axisymmetric flow, an appropriate superposition of the single-wake vertical and dispersive fluxes is performed, as it will be described in Sec.  \ref{sec:Verification}.

\section{Computational model}\label{sec:Computational}
The depth-averaged RANS equations (\ref{eq:RANS_EV}) are solved in terms of velocity components using a marching scheme approach in the streamwise direction, $x$, while the pressure field is calculated iteratively due to its inherently elliptic nature. Problems where the streamwise diffusion is negligible but the pressure gradient is not, are referred to as Partially-Parabolized Navier-Stokes (PPNS) \cite{Pletcher2012}, and provide advantages in terms of computational cost and data usage compared to fully elliptic solvers \cite{Pratap1975}. An essential feature of PPNS algorithms is the implementation of a pressure correction necessary to calculate the pressure field and enforce continuity to guarantee a divergence-free velocity field.

The PPNS scheme adopted for the P2D-RANS is inspired by previous works for ducted flows \cite{Briley1974,Pratap1975,Moore1979,Dodge1977}, with the main difference that no velocity correction is enforced throughout the sweeps. The global pressure field is indeed corrected at the end of each parabolic sweep \cite{Chilukuri1980}. The provisional velocity field does not satisfy mass conservation until the pressure correction achieves convergence (see Fig. \ref{fig:P2D_sketch}).
\begin{figure*}
\centerline{\includegraphics[width=1\textwidth, trim=0 0 0 0,clip]{P2D_sketch.pdf}}
 \caption{Workflow of the P2D-RANS.}\label{fig:P2D_sketch}
\end{figure*}

\begin{figure*}
\centerline{\includegraphics[width=0.35\textwidth, trim=0 0 0 0,clip]{P2D_grid.pdf}}
 \caption{Numerical grid of the P2D-RANS.}\label{fig:P2D_grid}
\end{figure*}

The numerical scheme adopted for the momentum equations is identical to that documented in \citep{Chilukuri1980}, except for the inclusion of the streamwise turbulent diffusion terms
in Eqs. (\ref{eq:RANS_SW}), which are retained to enhance accuracy and numerical stability. The scheme is fully implicit and the non-linear terms are treated via the lagging coefficient technique \cite{Pletcher2012}. The uniform grid is staggered so that each velocity node is placed between two pressure nodes (Fig. \ref{fig:P2D_grid}). The momentum equations in their discretized form are:
\begin{equation}\label{eq:uv_discretized}
    \begin{cases}
         \hat{u}_{i,j} \frac{\hat{u}_{i+1,j}-\hat{u}_{i,j}}{\Delta x}+\hat{v}_{i+1,j} \frac{\hat{u}_{i+1,j}-\hat{u}_{i+1,j-1}}{2\Delta y}+\hat{v}_{i+1,j+1} \frac{\hat{u}_{i+1,j+1}-\hat{u}_{i+1,j}}{2\Delta y}=\\
         =\frac{\hat{p}_{i+1,j}-\hat{p}_{i+2,j}}{\Delta x}+\nu_T \frac{\hat{u}_{i+2,j}-2\hat{u}_{i+1,j}+\hat{u}_{i,j}}{\Delta x^2}+\nu_T \frac{\hat{u}_{i+1,j+1}-2\hat{u}_{i+1,j}+\hat{u}_{i+1,j-1}}{\Delta y^2}+ \hat{f}_{x,i+1,j}\\
         \frac{\hat{u}_{i,j}+\hat{u}_{i,j-1}}{2} \frac{\hat{v}_{i+1,j}-\hat{v}_{i,j}}{\Delta x}+\hat{v}_{i,j} \frac{\hat{v}_{i+1,j+1}-\hat{v}_{i+1,j-1}}{2\Delta y}=\\
         =\frac{\hat{p}_{i+1,j-1}-\hat{p}_{i+1,j}}{\Delta y}+\nu_T\frac{\hat{v}_{i+2,j}-2\hat{v}_{i+1,j}+\hat{v}_{i,j}}{\Delta x^2}+\nu_T \frac{\hat{v}_{i+1,j+1}-2\hat{v}_{i+1,j}+\hat{v}_{i+1,j-1}}{\Delta y^2}+ \hat{f}_{y,i+1,j},
    \end{cases}
\end{equation}
where the indices $i$ and $j$ refer to the $i^\text{th}$ streamwise and $j^\text{th}$ spanwise locations, respectively. It is noteworthy that the 3D correction terms are tacitly included in the volumetric forces $\hat{f}$. The elliptic terms $\hat{u}_{i+2,j},\hat{v}_{i+2,j},\hat{p}_{i+2,j}$ are evaluated from the previous sweep. As indicated in Fig. \ref{fig:P2D_sketch}, Eqs. (\ref{eq:uv_discretized}) are solved parabolically marching from the inlet to the outlet and twice for each sweep, first without streamwise diffusion and then with the streamwise diffusion activated.

The pressure correction is derived from first principles in Appendix B. The final form of the Poisson equation for the calculation of the corrective pressure field, $p^c$, that is applied at the end of each sweep to correct the provisional velocity is:
 \begin{equation}\label{eq:Press_correction1}
 \begin{array}{l}
    \frac{\partial^2 \hat{p}^c}{\partial x^2}+\frac{\partial^2 \hat{p}^c}{\partial y^2}
     =\hat{u}\frac{\partial \nabla}{\partial x}+\hat{v}\frac{\partial \nabla}{\partial y}-\nu_T\left( \frac{\partial^2 \nabla}{\partial x^2}+\frac{\partial^2 \nabla}{\partial y^2}\right)
 \end{array}
 \end{equation}
 or, more concisely:
   \begin{equation}\label{eq:Press_correction2}
 \begin{array}{l}
    \nabla^2 \hat{p}^c
     =\frac{D\nabla}{Dt}-\nu_T\nabla^2\nabla,
 \end{array}
 \end{equation}
where $\nabla$ is the divergence of the provisional velocity field, including the correction for vertical mass fluxes. From a numerical standpoint, the pressure correction is discretized as:
 \begin{equation}\label{eq:p_discretized}
    \frac{\hat{p}^c_{i+1,j}-2\hat{p}^c_{i,j}+\hat{p}^c_{i-1,j}}{\Delta x^2}+\frac{\hat{p}^c_{i,j+1}-2\hat{p}^c_{i,j}+\hat{p}^c_{i,j-1}}{\Delta y^2}=s_{i,j}
 \end{equation}
 where:
   \begin{equation}
   \begin{array}{l}
         s_{i,j}=\hat{u}_{i-1,j} \frac{\nabla_{i,j}-\nabla_{i-1,j}}{2\Delta x}+\hat{u}_{i,j} \frac{\nabla_{i+1,j}-\nabla_{i,j}}{2\Delta x}+\hat{v}_{i,j} \frac{\nabla_{i,j}-\nabla_{i,j-1}}{2\Delta y}
         \\+\hat{v}_{i,j+1} \frac{\nabla_{i,j+1}-\nabla_{i,j}}{2\Delta y}
         -\nu_T\left(\frac{\nabla_{i+1,j}-2\nabla_{i,j}+\nabla_{i-1,j}}{\Delta x^2}+\frac{\nabla_{i,j+1}-2\nabla_{i,j}+\nabla_{i,j-1}}{\Delta y^2}\right)
     \end{array}
 \end{equation}
 and:
  \begin{equation}
    \nabla_{i,j}=\frac{\hat{u}_{i,j}-\hat{u}_{i-1,j}}{\Delta x}+\frac{\hat{v}_{i,j+1}-\hat{v}_{i,j}}{\Delta y} + \frac{\overline{w}_{2,i,j}-\overline{w}_{1,i,j}}{z_2-z_1},
 \end{equation}
where the last term represents the vertical mass flux obtained as explained in Sec.  \ref{sec:Conceptual}.

The Poisson equation for the pressure is solved through the $backslash$ Matlab solver for sparse linear systems \cite{HaglundElGaidi2018}. To speed up the convergence, the pressure correction is applied alternatively to all the cells and in blocks of 2$\times$2 cells \cite{Murthy1990}. Convergence is considered to be achieved when the mass residual, defined as the integral of the absolute value of the mass sources and sinks normalized by the incoming mass flow rate is below 1\% \citep{Chilukuri1980}.

Dirichlet boundary conditions on velocity ($\hat{u}=1, \hat{v}=0$) and pressure ($\hat{p}=0$) are prescribed at the inlet and outlet, respectively. Non-homogeneous Neumann conditions are applied to the side boundaries by imposing the spanwise gradients of the velocity components to be equal to their respective values at the nearest $y$ locations, which minimizes the back-propagation of boundary effects in the inner domain  \cite{Iungo2018}. Finally, homogeneous Neumann boundary condition is applied to the remaining boundaries for pressure, which is consistent with the Navier-Stokes equations \cite{Pletcher2012}. The spatial resolution of the grid is selected as $\Delta x=0.125$ and $\Delta y=0.05$, while the domain extends 15$D$ away from the turbines in all directions. These parameters are selected after a sensitivity analysis to ensure a maximum relative error on the streamwise component smaller than 1$\%$ in every point of the domain for both single-turbine and full wind farm simulations.

\section{Calibration}\label{sec:Calibration}
\subsection{Experimental Dataset}
To enhance accuracy in predictions of the intra-wind-farm velocity field and wind turbine power capture, the P2D-RANS has been developed by encompassing a fully data-driven actuator disk model and turbulence closure, both calibrated through wind LiDAR measurements of individual wakes generated by utility-scale wind turbines. The onshore wind farm under examination is located in North Texas and includes 25 identical 2.3 MW turbines with a rotor diameter of $D=108$ m and a hub height of $H=80$ m (Fig. \ref{fig:Panhandle_map_v2}). A WindCube 200S pulsed Doppler wind LiDAR equipped with a rotating scanning head was deployed from August 2015 to March 2017 (see Fig. \ref{fig:Panhandle_map_v2}). This section provides a concise overview of the field campaign, since a detailed description of the site, experimental strategy, and LiDAR data processing can be found in previous related publications \cite{El-Asha2017,Zhan2020,Maulik2021}.

A meteorological (met) mast is located at the South-East corner of the farm and recorded for the whole duration of the field campaign mean and standard deviation of horizontal wind speed and direction at 36-m, 60-m, and 80-m heights, barometric pressure at 2-m and 75-m heights, and temperature at 3-m and 75-m heights. Furthermore, SCADA data in the form of mean and standard deviation of hub-height wind speed, active power, and other operational parameters were acquired.
\begin{figure*}
\centerline{\includegraphics[width=\textwidth]{Panhandle_map_v2.pdf}}
\caption{Characterization of the wind farm under investigation: (a) topographic map of the site; (b) directional histogram of the normalized wind speed at hub height (Eq. (\ref{eq:U_norm})); (c) directional histogram of the incoming turbulence intensity at hub height. The green sectors in (a) represent the unwaked wind sectors, while the grey ones in (b) and (c) correspond to wind sectors where the met mast is likely affected by turbine wakes.}\label{fig:Panhandle_map_v2}
\end{figure*}

The wind LiDAR was set up with a range gate of 50 m and an accumulation time of 500 ms. The scanning wind LiDAR performed narrow Plan Position Indicator (PPI) scans with low elevation angles (typically $3^{\circ}$), azimuthal range of $20^{\circ}$, and rotation speed of the scanning head of $2^{\circ}$ s$^{-1}$, leading to a typical scanning time for a single PPI scan of 10 s. PPI scans were performed for wind directions included within the sector 145$^\circ$-235$^\circ$, when the wakes of turbines 1-6 are advected towards the LiDAR location. Overall, 16765 quality-controlled PPI scans are available for the present study.

Individual LiDAR samples undergo a quality control process, including data rejection based on the LiDAR signal intensity, which is performed based on the dynamic filtering \cite{Beck2017}. After realigning the wake according to the 10-minute-averaged estimated wind direction, the horizontal equivalent velocity field, $u_\text{eq}$, is made non-dimensional through the 10-minute-averaged incoming vertical wind profile, according to the procedure outlined in Zhan et al. \cite{Zhan2020}. Finally, the LiDAR measurements are reported in an axisymmetric-equivalent reference frame $(x,r)$, where the radial position, $r$, represents the distance from the expected wake center \citep{Zhan2020}. The re-aligned and non-dimensional horizontal-equivalent wake velocity data are clustered in bins based on operative conditions (i.e. aerodynamic rotor thrust) and incoming wind turbulence intensity. The selection of the bin edges is guided by a preliminary characterization of the wind farm performance as a function of the inflow conditions, which is the object of the following subsection.

\subsection{Wind turbine and wind farm performance}
The operation of a pitch-regulated turbine can be characterized as a function of the hub-height density-corrected wind speed \cite{IEC61400_12_1}, as follows:
\begin{equation}\label{eq:U_norm}
    U_\text{norm}=U_\text{hub}\left(\frac{\rho}{\rho_0}\right)^{1/3},
\end{equation}
where $U_\text{hub}$ is the 10-minute-averaged horizontal wind speed at hub height, $\rho$ is the measured air density, and $\rho_0=1.225$ kg m$^{-3}$ is the reference air density. Wake turbulent mixing and recovery can be conveniently investigated through the turbulence intensity, TI, defined as the ratio between the standard deviation of the hub-height wind speed over the mean of the same quantity. Specifically, Zhan et al \cite{Zhan2020} identified through an atmospheric stability analysis values of $7\%$ and $13.5\%$ as thresholds of the undisturbed turbulence intensity bounding stable (TI$_\infty<7\%$), near-neutral (TI$_\infty\in[7,13.5)\%$), and convective (unstable) (TI$_\infty\geq13.5\%$) conditions for the site under investigation. In general, the subscript $\infty$ refers to quantities recorded only by the SCADA data of turbines unaffected by wakes and, thus, representative of the incoming wind field. In this work, wind sectors affected by wakes are defined according to the 61400-12, Annex B of the IEC standards \cite{IEC61400_12_2}.

Wind speed at hub height, $U_\text{hub}$, is retrieved based on the measurements of the nacelle anemometer, which undergoes quality control and correction processes. Firstly, occurrences of unrealistic velocity and turbulence intensity, extreme veer, yaw misalignment, idling sensors, and zero or negative active power are rejected (see Table \ref{tab:Quality_check}). Then the nacelle transfer function is calculated based on met-tower normalized wind speed collected at 80-m height and the nacelle anemometer of Turbine 06 (namely the closest turbine to the met tower) only for unwaked wind sectors, in compliance with the International IEC Standards \cite{IEC61400_12_2}. Such transfer function cannot be applied directly to the whole park since biases of up to 5\% between different SCADA anemometers are present. The mentioned biases are estimated by comparing the raw individual power curves \cite{Sebastiani2020}. Therefore, all the wind speed data are re-scaled by a multiplicative factor using the data of Turbine 06 as a reference before the correction for the nacelle flow distortion.
\begin{table}[]
\centering
\caption{Breakdown of the quality check on met and SCADA data.}\label{tab:Quality_check}
\begin{tabular}{lll}
 & Formulation   & Rejection rate {[}\%{]} \\
 \hline
Veer                  & $\frac{d\theta_w}{dz}>0.25 ^\circ$ m$^{-1}$ & 16.3                    \\
Idling wind direction & $\frac{d \theta_w}{dt}=0$                  & 2.1                     \\
Yaw misalignment      & $|\theta_\text{yaw}-\theta_w|>30 ^\circ$   & 1.6                     \\
Idling power          & $\frac{d P_\text{norm}}{dt}=0$             & 0.04                    \\
No power              & $P_\text{norm}\leq0$                        & 9.9                     \\
Zero velocity         & $U_\text{norm}=0$                          & 0.02                    \\
Off-range TI          & TI$\leq0$ or TI$>100 \%$                    & 1.1\\
\end{tabular}
\end{table}

The experimental turbine power curve is calculated by leveraging the data from all the turbines following the IEC Standards \cite{IEC61400_12_2}. To flag power curtailments associated with off-design operations, an initial nominal power curve is generated based on the power of reference Turbine 06 only and the met-tower normalized wind speed and only for unwaked sectors. Power reads exhibiting a deficit larger than 300 kW from the power estimated through the nominal power curve ($\sim 1\%$ of the cases) are then rejected \citep{Hamilton2020}. The final power curves are generated using all the available SCADA data (August 2015 to April 2017) and differentiated as a function of TI \cite{Teng2020}. It is noteworthy that power in region III can exceed the nameplate value by virtue of the power-boost mode, which enables the generator to produce additional power under favorable operative conditions.
\begin{figure*}
\centerline{\includegraphics[width=0.9\textwidth]{Power_curve.pdf}}
 \caption{(a) Power curves for different stability classes. The grey dots represent quality-controlled 10-minute-averaged data. (b) Normalized Mean Absolute Error (NMAE) of the SCADA 10-minute-averaged power and the power obtained through several power curves and using different inputs as in incoming velocity.}\label{fig:Power_curve}
\end{figure*}
The TI-based normalized power ($P_\text{norm}$) curves are reported in Fig. \ref{fig:Power_curve}(a). The accuracy of this method is assessed quantitatively by evaluating the error between the measured power and the power obtained from the different power curves and using either met or SCADA (corrected) hub-height wind speeds as input variables. The error metric used here is the normalized mean absolute error \cite{Moriarty2014}, which is defined as:
\begin{equation}
    \text{NMAE}(x,x_\text{ref})=\frac{\sum_i|x_i-x_{i,\text{ref}}|}{\sum_i x_{i,\text{ref}}},
\end{equation}
where $x$ and $x_\text{ref}$ represent the modeled and reference value of the target quantity, respectively.
This error analysis indicates that the best agreement with the SCADA data is achieved by using the TI-based power curves and the corrected SCADA wind speed as input. Slightly worse accuracy is obtained by using the nominal power curve. The highest error associated with the met tower wind speed is due to the presence of spatial heterogeneity of the flow over the farm \cite{Hamilton2020}. Further, higher uncertainty is observed for the power output of waked turbines, likely due to the  more heterogeneous velocity distribution over the rotor area not captured by the nacelle anemometer.

Similar procedures are used to evaluate the power coefficient, $c_p$, and total farm power curve, which guide the selection of the bin edges of the undisturbed normalized hub-height velocity, $U_{\text{norm},\infty}$, used for the clustering of the LiDAR data. The bin selection is summarized in Fig. \ref{fig:Cluster_stats}. Specifically, the heat map in Fig. \ref{fig:Cluster_stats}(a) shows the occurrence of a $[U_{\text{norm},\infty},$ TI$_\infty]$ pairs based on the all available SCADA data and averaged over the set of unwaked turbines.
The green dots refer to the SCADA data paired with the wake LiDAR scans, which can be seen to cover satisfactorily most of the typical wind conditions, except for the $U_{\text{norm},\infty}<0.5$, TI$_\infty<13.5\%$ region. For this study, bin edges in $U_{\text{norm},\infty}$ are selected as $[0.25, 0.5, 0.7, 0.85, 1, 1.15, 2]$ (see Fig. \ref{fig:Cluster_stats}(b)), and are purposely refined in the proximity of the turbine rated power (i.e. $U_{\text{norm},\infty}=1$), where a stronger variation of $c_p$ and, thus, thrust coefficient, are expected. Finally, as highlighted by the farm power curve calculated including all the turbines regardless of their waked/unwaked state (Fig. \ref{fig:Cluster_stats}(c)), for $U_{\text{norm},\infty}>1.15$, the average total farm power is practically equal to the nominal capacity of the wind plant. The numerical modeling of such cases is irrelevant for the assessment of the proposed wind farm model in terms of power capture, and further bins for very high $U_{\text{norm},\infty}$ are thus not defined. Table \ref{tab:Cluster_points} reports the number of raw and quality-controlled experimental LiDAR samples available for each cluster.

\begin{figure*}
\centerline{\includegraphics[width=\textwidth]{Cluster_stats.pdf}}
\caption{Selection of the bin edges: (a) two-dimensional pdf of $U_{\text{norm},\infty}$ and TI$_\infty$ (green dots correspond to available LiDAR scans, while crosses indicate the bin centroids); (b) average $c_p$ experimental curves; (c) overall normalized farm power as a function of $U_{\text{norm},\infty}$, where the shaded areas represent standard deviation. The dashed lines separate the different bins in $U_{\text{norm},\infty}$.}\label{fig:Cluster_stats}
 \end{figure*}

\begin{table}
\caption{Number (in millions) and percentage (in brackets) of quality-controlled LiDAR samples available for each cluster.}\label{tab:Cluster_points}
\centering
\begin{tabular}{lccc}
& TI$_\infty <7\%$ & TI$_\infty \in [7,13.5)\%$ & TI$_\infty\geq13.5\%$\\
\hline
$U_{\text{norm},\infty}\in[0,0.25)$ & 0    & 0.04 (0.3)  &       0.24 (1.8)\\
$U_{\text{norm},\infty}\in[0.5,0.7)$ &0.29 (2.2) &0.88 (6.6) &0.52 (3.9) \\
$U_{\text{norm},\infty}\in[0.7,0.85)$ &0.77 (5.8) &1.32 (9.9) &0.42 (3.1) \\
$U_{\text{norm},\infty}\in[0.85,1)$ &1.77 (20.7) &2.87 (21.5) &0.3 (2.2)\\
$U_{\text{norm},\infty}\in[1,1.15)$ &0.58 (4.4) &1.46 (11.0) &0.16 (1.2) \\
$U_{\text{norm},\infty}\in[1.15,2)$ &0.21 (1.6) &0.48 (3.6) &0.03 (0.2) \\
\end{tabular}
\end{table}

\subsection{Data-driven calibration of the RANS model}
The mean streamwise velocity for each cluster is reconstructed using the LiDAR Statistical Barnes Objective Analysis (LiSBOA) tool \cite{Letizia2021,Letizia2021b} on a Cartesian grid with fundamental half-wavelengths $\Delta n_{0,x}=1.5 D,\Delta n_{0,r}= 0.5 D$ and region with local data spacing $\Delta \tilde{d}>0.25$ are rejected (refer to \cite{Letizia2021} for more theoretical details on LiSBOA). Data voids are filled through a bi-harmonic interpolation algorithm \cite{inpaintnans}. The velocity fields in axisymmetric coordinate are finally depth-averaged to build the reference dataset of $\hat{u}(x,y,U_{\text{norm},\infty},\text{TI}_\infty)$, which are reported in Fig. \ref{fig:Cal_u_LiDAR}. Results are not shown for two bins, which fail to satisfy the statistical constraints imposed by the LiSBOA due to the limited availability of data points. The used methodology for the characterization of the wake velocity field captures satisfactorily the reduced velocity deficit occurring for cases with lower $c_t$ (i.e. for $U_{\text{norm},\infty} \gtrapprox 0.85$), as well as the faster wake recovery occurring with increasing turbulence intensity. The quantification of the wake variability with incoming wind speed and turbulence intensity is instrumental for the calibration of the actuator disk and the turbulence closure model of the P2D-RANS and, more in general, for accurate modeling of the wind farm flow.
\begin{figure*}[b!]
\centerline{\includegraphics[width=\textwidth]{Cal_u_lidar.pdf}}
 \caption{Mean streamwise, depth-averaged velocity fields ($\hat{u}$) calculated from the LiDAR data for the various bins defined based on hub-height wind speed, $U_{\text{norm},\infty}$, and, turbulence intensity, $TI_\infty$.}\label{fig:Cal_u_LiDAR}
 \end{figure*}

The tuning of the turbulent eddy-viscosity and axial force distribution are performed following the methodology introduced by Iungo et al. \cite{Iungo2018}, but adopting depth-averaged velocity fields and the P2D-RANS model instead of the axisymmetric flow and RANS model of the original publication. The present calibration procedure includes the following steps:
\begin{itemize}
     \item the near wake length, $x_\text{NW}$, namely the extension of the region where the pressure gradients are relevant, is estimated as the streamwise location where the  minimum rotor-averaged velocity occurs \cite{Vermeulen1980};
     \item a constant turbulent eddy-viscosity is estimated by minimizing the mean absolute percentage error (MAPE) between the average streamwise wake velocity predicted with the P2D-RANS and that obtained from the LiDAR data, where:
     \begin{equation}
         \text{MAPE}=\bigg\langle\frac{|\hat{u}_\text{P2D-RANS}(x,|y|<1)-\hat{u}_\text{LiDAR}(x,|y|<1)|}{\hat{u}_\text{LiDAR}(x,|y|<1)}\bigg\rangle,
     \end{equation}
     with $\langle \rangle$ indicating space average. In this case, the modified P2D-RANS solver does not include the rotor region, instead the $\hat{u}$ velocity profile calculated from the LiDAR data at $x_\text{NW}$ is injected as inlet boundary condition, while assuming constant pressure to calculate $\hat{v}$. The optimization algorithm used to estimate turbulent eddy-viscosity is described in Iungo et al.\cite{Iungo2018};
     \item once the optimal $\nu_T$ is known, the axial thrust force on the turbine rotor is estimated by correcting iteratively the streamwise forcing, $\hat{f}_x$, in the full version of the P2D-RANS model to match the experimental velocity profile at $x_\text{NW}$ as:
     \begin{equation}
     \begin{array}{l}
     \int_{-\infty}^{+\infty} \hat{f_x}^{n+1}(x,y)~dx=\int_{-\infty}^{+\infty} \hat{f_x}^{n}(x,y)~dx
     +\hat{u}_\text{P2D-RANS}(x_{NW},y)^2-\hat{u}_\text{LiDAR}(x_\text{NW},y)^2,
     \end{array}
     \end{equation}
     where the index $n$ is the iteration counter. This method represents a simplified version of the technique based on the local momentum budget \cite{Iungo2018}.
     It is noteworthy that the integrals in the $x$-direction are included only to compensate for the streamwise smoothing of the forces applied in the P2D-RANS to prevent numerical instabilities (further details on the formulation of the aerodynamic loads are provided in Appendix \ref{app:Loads}). The algorithm stops when variation in the above-defined MAPE lower than 0.1\% are attained (generally with less than 5-6 iterations).
\end{itemize}

The results of the calibration procedure are summarized in Figs. \ref{fig:Cal_EV_Ct} and \ref{fig:Cal_loads}.
\begin{figure*}[b]
\centerline{\includegraphics[width=\textwidth]{Cal_EV_Ct.pdf}}
 \caption{Calibration of the P2D-RANS: (a) optimal turbulent eddy-viscosity; (b) optimal thrust coefficient, $c_t$; (c) mean absolute percentage error (MAPE) (empty bars represent the MAPE for the far wake velocity field after the optimization of $\nu_T$ only, and the full ones represent the final MAPE after the estimation of the axial load).}\label{fig:Cal_EV_Ct}
 \end{figure*}
 \begin{figure*}[b!]
\centerline{\includegraphics[width=\textwidth]{Cal_loads.pdf}}
 \caption{Optimal axial forcing distribution after calibration: (a) TI$_\infty<7\%$; (b) TI$_\infty$ $\in[7,13.5)\%$; (c) TI$_\infty\geq13.5\%$.}\label{fig:Cal_loads}
 \end{figure*}
The behavior of the optimal turbulent eddy-viscosity (Fig. \ref{fig:Cal_EV_Ct}(a)) lends interesting insight into the wake recovery process. The most consistent pattern is the increase of $\nu_T$ as a function of TI$_\infty$, which is expected since TI$_\infty$ is linked to the magnitude of the turbulent momentum fluxes in the wake. The turbulent eddy viscosity increases by nearly tenfold in convective conditions compared to stable conditions, which confirms the necessity of including effects of atmospheric stability in wind farm models. Furthermore, the turbulent eddy-viscosity for low and moderate TI$_\infty$ decreases as a function of $U_{\text{norm},\infty}$, except for the last bin $U_{\text{norm},\infty}\in [1.15,2)$. Our interpretation is that the $\nu_T$ optimizer identifies effects related to the wake-added turbulence, which is indeed stronger for higher $c_t$ values (i.e. low $U_{\text{norm},\infty}$) conditions due to the enhanced shear-generated turbulence. The proportionality of the turbulent eddy viscosity to the velocity deficit, $\Delta u$, (i.e. $c_t$) is predicted by the shear flow theory, which produces the well-known scaling $\nu_T \sim r_w \Delta u$, where $r_w$ is the wake radius, for a self-similar axisymmetric wake \cite{Sforza1981,Ainslie1988,Lange2003}.

It is noteworthy that $\nu_T$ for convective conditions increases between $U_{\text{norm},\infty}=0.5$ to $1.15$. This trend, which is reversed compared to what is observed for low to moderate incoming TI$_\infty$ conditions, is likely due to a different wake recovery mechanism occurring for unstable atmospheric conditions, where the incoming turbulence significantly contributes to the momentum diffusion in the wake, especially through the flow dynamics connected with wake meandering \cite{Larsen2008}. Furthermore, the inverse proportionality between $\nu_T$ and $c_t$ observed in convective conditions may be due to the damping of the largest incoming coherent structures operated by the rotor already documented by Chamorro et al. \cite{Chamorro2012}, which is expected to be more severe for high $c_t$. A more thorough investigation of this interesting behavior will be the object of future work. Finally, the reduction of $\nu_T$ observed for the last bin $U_{\text{norm},\infty}\in[1.15,2)$ is ascribed to the very shallow wake deficit, which makes the definition of an eddy viscosity quite elusive

Fig. \ref{fig:Cal_EV_Ct}(b) and \ref{fig:Cal_loads} show $c_t$ and the spanwise distribution of depth-averaged thrust force. The thrust coefficient is practically constant for $0.25 \leq U_{\text{norm},\infty}<0.7$ (i.e. region II of the power curve), then it decreases due to the blade pitching and tip-speed ratio reduction operated by the controller to limit the power output to the nameplate capacity. The only available thrust-force profile at very low $U_{\text{norm},\infty}$ corresponds to convective conditions, which shows an overall higher aerodynamic force. The force distribution is less sensitive to the atmospheric stability, although slightly higher values are estimated in high TI conditions \citep{Iungo2018}

Fig. \ref{fig:Cal_EV_Ct}(c) shows a MAPE that lies below the 5\% threshold for all the wake after the calibration, which is considered satisfactory. To conclude, the calibrated numerical fields are shown in Fig. \ref{fig:Cal_u_P2D}, while the difference fields (P2D-RANS minus LiDAR) are reported in Fig. \ref{fig:Cal_u_diff}. Finally, Tables \ref{tab:EV} and \ref{tab:Ct} provide the $\nu_T$ and $c_t$ resulting from the calibration procedure.
 \begin{figure*}
\centerline{\includegraphics[width=\textwidth]{Cal_u_P2D.pdf}}
 \caption{Mean streamwise, depth-averaged velocity fields ($\hat{u}$) of the P2D-RANS after calibration.}\label{fig:Cal_u_P2D}
 \end{figure*}
  \begin{figure*}
\centerline{\includegraphics[width=\textwidth]{Cal_u_diff.pdf}}
 \caption{Difference $\hat{u}$ velocity field between P2D-RANS and LiDAR after calibration.}\label{fig:Cal_u_diff}
 \end{figure*}

\begin{table}
\caption{Optimal turbulent eddy viscosity of the P2D-RANS model. The value is made non-dimensional by $U_\infty D$.}\label{tab:EV}
\centering
\begin{tabular}{lccc} & TI$_\infty <7\%$ & TI$_\infty \in [7,13.5)\%$ & TI$_\infty\geq13.5\%$\\
\hline
$U_{\text{norm},\infty}\in[0,0.25)$ & N/A    & N/A  &       0.0575\\
$U_{\text{norm},\infty}\in[0.5,0.7)$ &0.0082 &0.0139 &0.0318 \\
$U_{\text{norm},\infty}\in[0.7,0.85)$ &0.0095 &0.0133 &0.0356 \\
$U_{\text{norm},\infty}\in[0.85,1)$ &0.0057 &0.0097 &0.0499 \\
$U_{\text{norm},\infty}\in[1,1.15)$ &0.003 &0.0051 &0.0663 \\
$U_{\text{norm},\infty}\in[1.15,2)$ &0.0081 &0.0176 &0.0323 \\
\end{tabular}
\end{table}

\begin{table}
\caption{Coefficient of thrust calibrated on LiDAR data.}\label{tab:Ct}
\centering
\begin{tabular}{lccc}& TI$_\infty <7\%$ & TI$_\infty \in [7,13.5)\%$ & TI$_\infty\geq13.5\%$\\
\hline
$U_{\text{norm},\infty}\in[0,0.25)$ & N/A    & N/A  &       1.83\\
$U_{\text{norm},\infty}\in[0.5,0.7)$ &1.21 &1.28 &1.51 \\
$U_{\text{norm},\infty}\in[0.7,0.85)$ &1.3 &1.24 &1.52 \\
$U_{\text{norm},\infty}\in[0.85,1)$ &0.81 &0.8 &1.11 \\
$U_{\text{norm},\infty}\in[1,1.15)$ &0.46 &0.48 &0.86 \\
$U_{\text{norm},\infty}\in[1.15,2)$ &0.23 &0.37 &0.64 \\
\end{tabular}
\end{table}
For the implementation of the P2D-RANS, the values of the calibrated rotor thrust coefficient and turbulent eddy-viscosity are assigned to the bin centroids displayed in Fig. \ref{fig:Cluster_stats}(a) and interrogated through a linear/nearest-neighbor Delauney interpolation to perform simulations with an arbitrary inflow.

\section{Verification}\label{sec:Verification}
According to an AIAA technical report \cite{AIAA2002}, the process of determining if a model implementation accurately represents the conceptual description and the solution of a model is referred to as verification. The purpose of this section is indeed to verify the validity of the proposed shallow-water approach corrected for vertical fluxes and dispersive terms, the pressure correction algorithm of the P2D-RANS, and the numerical method. To this aim, the P2D-RANS solution of a single-wake velocity field is compared to the equivalent result calculated through an elliptic axisymmetric RANS solver. The latter uses the Newton algorithm to solve the full non-linear set of RANS equations in cylindrical coordinates discretized through a spectral collocation method based on Chebyshev polynomials (for more details, please refer to \citep{Viola2014,Iungo2018}, while the code is publicly available \cite{GRANS}. The simulations with the axisymmetric code are carried out using equivalent boundary conditions, turbine load, and turbulence model. The inflow conditions are arbitrarily chosen as $U_\infty=10$ m s$^{-1}$, TI$_\infty=10\%$, which in turn define the axial aerodynamic thrust and the turbulent eddy viscosity (Sec.  \ref{sec:Calibration}).

Fig. \ref{fig:Verification_prof} shows the spanwise profiles of $\hat{u}, \hat{v}, \hat{p}$ extracted at several downstream locations calculated by the axisymmetric solver (dots) and the P2D-RANS (solid lines), as well as the same quantities for a purely 2D solver, i.e. without 3D corrections (dashed lines). The agreement between the axisymmetric solver (i.e. the benchmark) and the P2D-RANS is excellent, with negligible discrepancies, which are due to the slightly different implementation of the boundary conditions for the two models \cite{Iungo2018}. The almost perfect overlapping between the pressure obtained by the axisymmetric RANS solving the fully elliptic Navier-Stokes equations and through the iterative marching scheme of the P2D-RANS described in Sec.  \ref{sec:Computational} serves as verification of the proposed pressure-velocity coupling scheme. Conversely, the profiles obtained with the merely 2D-RANS code exhibit significantly higher speedups, hampered wake recovery, excessive spanwise velocity magnitude, and pressure gradients compared to the benchmark simulations, which justifies the utilization of the proposed 3D corrections to relax the vertical confinement, which is a typical limitation of 2D wind farm flow simulators.
\begin{figure*}
\centerline{\includegraphics[width=\textwidth]{Verification_prof.pdf}}
 \caption{Comparison between the depth-averaged fields of the elliptic axisymmetric RANS solver, a purely 2D-RANS code, and the P2D-RANS with 3D corrections: (a) streamwise velocity component; (b) spanwise velocity component; (c) pressure.}\label{fig:Verification_prof}
\end{figure*}

The 3D correction terms introduced in the P2D-RANS are assessed by plotting the individual contributions of the vertical fluxes and dispersive terms (overbarred and tilded terms, respectively, in Eqs. (\ref{eq:DepthAverage})) to the mass (Fig. \ref{fig:Verification_3Dcorr}(a)), $x$-momentum (Fig. \ref{fig:Verification_3Dcorr} (b, d, f)) and $y$-momentum (Fig. \ref{fig:Verification_3Dcorr} (c, e, g)) budgets. For an isolated turbine, these terms peak in the rotor area and drop sharply outside of the rotor span, decay further downstream as the wake recovers, and cross-stream components dissipate.
\begin{figure*}
\centerline{\includegraphics[width=\textwidth]{Verification_3Dcorr.pdf}}
 \caption{3D corrections for the single turbine case: (a) vertical mass flux; (b) streamwise vertical advection; (c) spanwise vertical advection; (d) streamwise vertical turbulent flux; (e) spanwise vertical turbulent flux; (f) streamwise dispersive stresses; (g) spanwise dispersive stresses.}\label{fig:Verification_3Dcorr}
\end{figure*}
It is noteworthy that the results of the present verification hold valid for an arbitrary turbine thrust and eddy viscosity, namely for the full range of inflow conditions defined in the calibration phase (Sec.  \ref{sec:Calibration}).

For the simulation of a wind farm flow, the assumption of axially symmetrical and non-swirling flow underpinning the  method for the estimation of the 3D effects described in Appendix \ref{app:3D_corr} fails due to the merging of multiple wakes. Therefore, it is necessary to formulate an appropriate superposition method of the vertical fluxes and dispersive terms of the single turbines. To this aim, the flow resulting from the interaction of two in-tandem turbines with a streamwise spacing of $5D$ is calculated through the axisymmetric elliptic RANS and the P2D-RANS. The latter estimates the overall vertical mass and momentum fluxes through a superposition of those obtained for single-wake simulations and stored in a look-up table as a function of $U_{\text{norm},\infty}$ and TI$_\infty$ and interrogated through a nearest-neighbor method. Namely, the 3D corrections are evaluated for 15 values of $U_{\text{norm},\infty}$ from 0.2 to 2 and 12 values of TI$_\infty$ from 4\% to 50\% (i.e. 180 total single-wake cases), which cover the 98.5\% of the climatology of the Panhandle wind farm. The domain considered for single-wake fluxes extends from $-4D$ to $20D$ in the streamwise direction and $\pm 2D$ in the spanwise direction, since vertical fluxes and dispersive terms beyond such region are negligible for a wide range of operative conditions (see Fig. \ref{fig:Verification_3Dcorr}). Different methods have been considered for the estimate of the total fluxes resulting from wake overlapping, such as linear sum, root-squared signed sum, and the maximum absolute value of the single-wake contributions. The error with respect to the axisymmetric RANS, which solves directly the full flow equations without any simplifying assumption, is quantified to identify the best performing superposition principle. A total of 18 simulations of two in-tandem wind turbines with $5D$ streamwise spacing are performed for each method, each one having as input the centroid of the clusters in terms of $U_{\text{norm},\infty}$ and TI$_\infty$ previously defined. This analysis aims to characterize the flow-dependent error for various wake conditions.

The MAPE of $\hat{u}$ for all the simulations is provided in Fig. \ref{fig:Verification_VFS} and indicates the root squared method with sign as that exhibiting the overall best agreement, with error exceeding 2\% for only one case (stable, region II).
\begin{figure*}
\centerline{\includegraphics[width=\textwidth]{Verification_VFS.pdf}}
 \caption{Mean absolute relative error of $\hat{u}$ for between the axisymmetric RANS and the P2D-RANS using different superposition methods for the 3D corrections: (a) linear sum; (b) root-squared signed sum; (c) maximum absolute value.}\label{fig:Verification_VFS}
\end{figure*}
It is noteworthy that the MAPE in $\hat{u}$ due to the superposition of the 3D corrections is always lower than 3\%, which is quite smaller than the errors reported for traditional engineering wake models, for which these types of superposition methods are applied directly on the velocity field \cite{Gunn2016,Machefaux2015}.

Figure \ref{fig:Verification_example} shows three selected cases with two in-tandem turbines with $5D$ streamwise spacing for the centroid of the cluster $U_{\text{norm},\infty}\in (0.75, 0.8]$ (optimal $c_t$) and different TI$_\infty$ simulated through the reference elliptic RANS solver and the P2D-RANS with the root-squared signed sum method. The difference fields (Fig. \ref{fig:Verification_example} (g-i)) show that a systematic velocity underestimation takes place in the near wake region for the P2D-RANS, where vertical fluxes and dispersive terms peak, while the error decreases in magnitude and switches sign moving downstream. These features are observed in most of the simulations performed for this analysis and indicate that the empirical superposition technique of the 3D corrections may lead to substantial errors only in case of very close wake interaction (e.g. with streamwise spacing smaller than $3D$), which is however unlikely for well-designed wind farms. A more in-depth investigation on the accuracy of the 3D corrections superposition might be the object of future work.

\begin{figure*}
\centerline{\includegraphics[width=\textwidth]{Verification_example.pdf}}
 \caption{Streamwise velocity simulated for two in-tandem wind turbines with a streamwise spacing of $5D$, $U_{\text{norm},\infty}=0.775$, and different turbulence intensity levels. By row: (a-c) axisymmetric RANS; (d-f) P2D-RANS using the root-squared superposition method; (g-i) difference (P2D-axisymmetric RANS). By columns: (a, d, g) TI$_\infty<7\%$; (b, e, h) TI$_\infty\in(7,13.5]\%$; (c, f, i) TI$_\infty \geq13.5\%$. }\label{fig:Verification_example}
\end{figure*}

\section{Validation}\label{sec:Validation}
The validation of a computational model is formally defined as "the process of determining the degree to which a model is an accurate representation of the real world from the perspective of the intended uses of the model"\cite{AIAA2002}. In our case, real-world data are available in the form of 10-minute-averaged power from the SCADA data of the wind farm under investigation. The applications of the P2D-RANS model, considering its potential in accuracy and computational costs, can be broadly classified into two categories: $i)$ real-time diagnostic and control of wind turbines; $ii)$ wind farm performance analysis and optimal design. The former focuses on the estimation of the time variability of power production (typically in the form of 10-minute-averaged quantities) and accuracy of the P2D-RANS is assessed by comparing model predictions with the SCADA data for specific time series. The second application mainly entails the evaluation of wind turbine and wind farm performance in terms of power capture and quantification of wake losses as a function of atmospheric and turbine operational parameters. In this case, the following performance indicators are evaluated:
\begin{itemize}
\item percentage power losses \cite{El-Asha2017}, which are more relevant from the single-turbine modeling standpoint, and are defined as follows:
\begin{equation}\label{eq:Power_loss}
\Delta P_i(t_j)=1-\frac{P_i(t_j)}{P_{\infty,i}(t_j)},
\end{equation}
where $P_i(t_j)$ and $P_{\infty,i} (t_j)$ are the actual and unwaked power of the $i^\text{th}$ turbine at time $t_j$, respectively;
\item the wind farm efficiency, which is analogous to the previous parameter, yet for the entire turbine array \cite{Walker2016}:
\begin{equation}\label{eq:WF_eff}
    \eta=\frac{\sum_{i=1}^{N_t} P_i(t_j)}{\sum_i^{N_t} {P_{\infty,i}(t_j)}},
\end{equation}
with $N_t=25$ turbines for this study;
\item the annual (or cumulative) energy loss \cite{El-Asha2017}, which is more relevant from a financial standpoint since it considers the occurrence of each flow scenario within the local climatology, which can be recast as:
\begin{equation}\label{eq:Energy_loss}
   \Delta E=\frac{1}{N}\sum_{j=1}^N \sum_{i=1}^{N_t}(P_{\infty,i}(t_j)-P_i(t_j)) \cdot 8760,
\end{equation}
where $N$ is the number of available quality controlled data for the flow conditions under consideration and the unit is kWh/y.
\end{itemize}
The unwaked power for each turbine, $P_{\infty,i}(t_j)$, is defined as the ideal power capture in absence of wake interactions, and it is estimated through a method similar to that proposed by Farrel et al. \cite{Farrell2021} to map a heterogeneous wind velocity field over a farm from pointwise sparse measurements. In this case, the same technique is applied to the power, using the power of unwaked turbines as input. Namely, the power from unwaked turbines (see Fig. \ref{fig:Panhandle_map_v2}(a)) is interpolated (extrapolated) at the location of waked turbines through a linear (nearest-neighbor) interpolation approach after a Delaunay triangulation. This ensures that the gross power takes into account possible heterogeneity in the wind field over the farm, whose relevance was already confirmed through the power curve analysis (section \ref{sec:Calibration}).

It is noteworthy that the occurrence of overpower conditions, i.e. $P_i>P_{\infty,i}$, for a wind farm in flat terrain can be due either due to speedups caused by the induction pressure field of neighboring turbines or to an incorrect estimation of $P_{\infty,i}$, which might be associated with low availability of unwaked turbines for a certain wind direction. To exclude this last type of outliers yet retain cases with speedups, cases with overpower higher than 10\% are discarded only if the P2D-RANS model indicates that no pressure-induced speedup occurs at that turbine. In other words, the data point is rejected if SCADA data indicate $P_i/P_{\infty,i}>1.1$ and the P2D-RANS shows $P_i/P_{\infty,i}\leq1$ (1.3\% rejection rate).

Similarly to power, the undisturbed hub-height wind speed, $U_{\text{norm},\infty}$, turbulence intensity, TI$_\infty$, and wind direction, $\theta_w$, are estimated as averages of the respective quantities recorded by all the nacelle anemometers ($U_\text{norm}$ and TI$_\infty$) and yaw encoders ($\theta_w$) over the set of unwaked turbines at a specific time, $t_j$. Cases exhibiting error on the mean over the unwaked turbines at time $t_j$ based on the Student's $t$ distribution \cite{Wheeler2004} higher than 0.25 for $U_{\text{norm},\infty}$ (0.5\% of the overall datasets), 25\% of the average value for TI$_\infty$ (16\%), and 5$^\circ$ for $\theta_w$ (1.9\%) consequent to high flow variability within the park are discarded. Furthermore, periods when more than two turbines are either curtailed (based on the criteria described in Sec.  \ref{sec:Calibration}) or likely affected by wakes from neighboring turbine arrays that are not modeled are entirely excluded (35\%). This slightly relaxed criterion for the rejection of off-design farm operation is beneficial in terms of data availability and has minimal impact on the final statistics \cite{Nygaard2015}.

To showcase the capabilities of the P2D-RANS and highlight the challenges connected with the predictions of power production for individual onshore wind turbines, the simulations of a full day of turbine operation are first analyzed. The simulations are carried out for operations on August 26$^\text{th}$, 2016, and use the previously defined undisturbed hub-height wind speed, direction, and turbulence intensity as input with a 10-minute resolution. The selected day exhibits the typical stability-driven diurnal cycle of the wind resource \cite{El-Asha2017}, as shown in Fig. \ref{fig:Time_series_Inflow}.
\begin{figure*}
\centerline{\includegraphics[width=\textwidth]{Time_series_Inflow.pdf}}
 \caption{Inflow characterization for August 26$^\text{th}$, 2016: (a) hub-height normalized wind speed; (b) wind direction; (c) turbulence intensity; (d) total normalized farm power capture. The background colors in (a) and (c) facilitate the identification of the $U_{\text{norm},\infty}$ and TI$_\infty$ clusters used in this work.}\label{fig:Time_series_Inflow}
 \end{figure*}
Specifically, the SCADA-recorded and P2D-RANS simulated power capture indicate a generally good agreement, except around 06:00 and 21:00 when the power is overestimated with the P2D-RANS, which is probably a consequence of the inherently complex transition in the atmospheric stability regime. Further insight is gained by selecting four illustrative cases indicated with vertical dashed lines in Fig. \ref{fig:Time_series_Inflow}. For these cases, experimental and numerical wind speed and power capture are reported in Figs. \ref{fig:Time_series_LiDAR} and \ref{fig:Time_series_Power}, respectively.
\begin{figure*}
\centerline{\includegraphics[width=\textwidth]{Time_series_LiDAR.pdf}}
 \caption{Horizontal equivalent LiDAR (top row) and depth-averaged P2D-RANS (bottom rows) streamwise velocity for selected 10-minute periods on August 26$^\text{th}$, 2016: (a, e) Case 1, 1:10-1:20, $U_{\text{norm},\infty}=1.07$, $\theta_w=194^\circ$, TI$_\infty=5.3\%$;  (b, f) Case 2, 7:10-7:20, $U_{\text{norm},\infty}=0.85$, $\theta_w=262^\circ$, TI$_\infty=8.5\%$; (c, g) Case 3, 14:40-14:50, $U_{\text{norm},\infty}=0.55$, $\theta_w=202^\circ$, TI$_\infty=21.2\%$;  (d, h) Case 4, 18:50-19:00, $U_{\text{norm},\infty}=0.74$, $\theta_w=156^\circ$, TI$_\infty=6.8\%$;}\label{fig:Time_series_LiDAR}
 \end{figure*}
 \begin{figure*}
\centerline{\includegraphics[width=\textwidth]{Time_series_Power.pdf}}
 \caption{As in Fig. \ref{fig:Time_series_LiDAR}, but for SCADA (top row) and P2D-RANS simulated (bottom row) power capture.}\label{fig:Time_series_Power}
 \end{figure*}

Case 1 represents a typical wind farm operation in region III of the power curve under stable atmospheric conditions and southerly wind direction.  Relatively long wakes are observed with a weak velocity deficit (Fig. \ref{fig:Time_series_LiDAR} (a, e)) leading to significant and localized power losses (Fig. \ref{fig:Time_series_Power} (a, e)). The main challenges in modeling this type of wind farm operations arise from the uncertainty in the quantification of the local wind direction throughout the farm, which is a major source of error in case of persistent wakes with limited lateral expansion.

Case 2 shows the occurrence of strong wake interactions with relatively small streamwise spacing ($\sim 3D$) associated with moderate westerly winds for operations in region II of the power curve. The proximity of the waked rotors to the complex near-wake turbulent flow of multiple upstream turbines is arguably the most challenging (yet unlikely for a well-designed farm layout) case to simulate. Nevertheless, the P2D-RANS can reproduce satisfactorily the pattern of power losses across the farm (Fig. \ref{fig:Time_series_Power}(f)), albeit an evident overestimation in power capture is observed for specific turbines, which may be attributed to the detrimental effect on power capture of partial and/or intermittent wake interaction not reproduced by the steady RANS model.

Case 3 corresponds to highly convective daytime conditions, as portrayed by the utterly complex instantaneous flow field probed by the LiDAR (Fig. \ref{fig:Time_series_LiDAR}(c)). Thanks to the optimally calibrated $\nu_T$, the P2D-RANS model is capable in reproducing the enhanced turbulent mixing of the wakes, which results in moderate yet widespread power losses (Fig. \ref{fig:Time_series_Power}(g)), in agreement with the respective SCADA data (Fig. \ref{fig:Time_series_Power}(c)). The latter plots also show the occurrence of significant spanwise wind gradient at the southern-most array facing the incoming wind field, which is a quite common phenomenon at this site that justified the formulation of the turbine-specific unwaked power, $P_{\infty,i}$, adopted in Eqs. (\ref{eq:Power_loss}, \ref{eq:WF_eff}, \ref{eq:Energy_loss}).

Finally, case 4 represents an interesting south-east wind event, which results in $\sim 80\%$ power loss at turbine 07, which is accurately simulated by the P2D-RANS (Fig. \ref{fig:Time_series_Power} (d, h)).

This qualitative analysis is meant to remark the considerable complexity of a real wind farm flow and the challenges inherent to the simulation of the 10-minute-averaged power at the single turbine level, but also to stress the ability of the P2D-RANS model to reproduce the physical behavior of individual turbines experiencing a real-world turbulent ABL flow, yet keeping the required computational costs and resources very limited.

More systematic quantification of the uncertainty of the P2D-RANS model is carried out through a statistical analysis encompassing all the available SCADA data. To this aim, a look-up table of simulated single-turbine power is built for a wide range of inflow conditions covering the overall climatology of the site. The input inflow matrix is the combination of several $U_{\text{norm},\infty}$ values between 0.2 and 2 (15 cases), TI$_\infty$ between 4\% and 50\% (12 cases), and wind directions from $0^\circ$ to $360^\circ$ with a resolution of $5^\circ$, resulting in 12960 simulations and a 98.5\% coverage of the climatic conditions of the wind farm. A further refinement of the input matrix did not produce any significant differences in the results. The simulated power at the generic time $t_j$ is obtained by interrogating independently for the $i^\text{th}$ turbine the power look-up table, $P_i(U_{\text{norm},\infty}, \text{TI}_\infty,\theta_w)$, through a tri-linear interpolation with the experimental values of $U_{\text{norm},\infty}(t_j)$, TI$_\infty(t_j)$ and $\theta_w(t_j)$ as inputs.

The characterization of the accuracy of the P2D-RANS in simulating the time-resolved 10-minute-averaged power is carried out through a linear regression of the power, as reported in Fig. \ref{fig:Validation_Power}.
\begin{figure*}
\centerline{\includegraphics[width=0.9\textwidth]{Validation_Power.pdf}}
 \caption{Linear regression of SCADA vs. P2D-RANS 10-minute-averaged power: (a) normalized wind farm power; (b) normalized turbine power; (c) normalized power of unwaked turbines; (d) normalized power of waked turbines.}\label{fig:Validation_Power}
 \end{figure*}
There is an outstanding agreement for the overall power of the entire farm (Fig. \ref{fig:Validation_Power}(a)) with a negligible bias and $R^2=0.998$, which drops when considering the regression between the power of individual turbines (Fig. \ref{fig:Validation_Power}(b), $R^2=0.956$). This feature indicates that a noticeable error cancellation occurs at the farm level, thus making the actual accuracy in power predictions for single turbine power a better indicator to assess the validity of the proposed flow model. For a better understanding of the error source, Figs. \ref{fig:Validation_Power} (c, d) report the linear regression of single-turbine power for unwaked and waked wind turbines, respectively. It is noteworthy that, although the prediction of the power of waked turbines exhibits larger discrepancies compared to that of unwaked turbines ($R^2=0.945$ vs. $R^2=0.974$), the latter still is encumbered with a non-negligible scattering. This mismatch between the predicted and measured power for unwaked turbines is due to spatial wind gradients, unsteady conditions within the 10-minute periods, power curve uncertainty, and undetected wakes, which represent a baseline error that is inevitably present regardless of the accuracy of the used wind farm model. The linear regression of power capture under waked conditions shows an almost unitary slope and nearly zero intercept, which indicates a negligible contribution of systematic biases in the prediction of the wake interactions to the overall uncertainty.

To provide a more focused statistical analysis of power and energy losses due solely to wake interactions, further validation is performed also including results obtained with the Jensen \cite{Jensen1983}, Multizone \cite{Gebraad2016}, and Gaussian \cite{Bastankhah2014,Niayifar2016} wake models implemented in the FLORIS package \cite{FLORIS_2021}. All these models use the nominal power curve (see Fig. \ref{fig:Power_curve}) and the optimally calibrated $c_t$ reported in Sec.  \ref{sec:Calibration}, while other wake parameters have default FLORIS settings (specific documentation and references are available at \cite{FLORIS_2021}).

Firstly, the ability of the P2D-RANS model to predict the directional wake losses is tested by calculating the bin-average of the percentage power losses for each turbine (Eq. \ref{eq:Power_loss}) in $5^\circ$-wide sectors, as shown in Fig. \ref{fig:Validation_PL}. Statistics for bins having an error on the mean larger than $25\%$ with $95\%$ confidence are excluded.
\begin{figure*}[b!]
\centerline{\includegraphics[width=\textwidth]{Validation_PL.pdf}}
 \caption{Polar chart of the directional percentage power losses. The blue bars represent the SCADA data, while the red lines are the P2D-RANS results and the black (zoomed region only) the Gauss model of FLORIS. The grey dots represent the turbines, with the radius corresponding also the $\Delta P_i=0\%$ losses. The dashed circles indicate $\Delta P_i=100\%$.}\label{fig:Validation_PL}
 \end{figure*}
The model can reproduce the directional pattern of power losses of each turbine satisfactorily with a NMAE $=30.5\%$ with respect to the SCADA data. In Fig. \ref{fig:Validation_PL}, the zoomed-in frame shows how the model can accurately reproduce the wake losses resulting from both near (between turbine 07 and 08, with spacing $4.7D$) and far (South direction for turbine 08, with spacing $\sim 15 D$) wake interactions, thus confirming the validity of the wake model for a wide range of streamwise distances, even exceeding the domain of LiDAR data used for the calibration. The zoomed region also indicates that the P2D-RANS provides accurate quantification of the physically consistent overpower occurring at turbine 08 for wind directions at the edge of a heavily-waked eastern sector (overpower is visualized as a negative power loss, i.e. the part of the polar plot that lies inside the grey circle corresponding to zero losses). For these cases the P2D over-performs other engineering wake models, which lack modeling of the pressure field, as the Gaussian wake model (black line in the same plot), and, thus, any effect on the incoming wind field induced by the thrust force of a turbine rotor.

The analysis of the wind farm performance as a function of the incoming flow is further refined by bin-averaging the wind farm efficiency, $\eta$, for all the wind directions and separately for each $U_{\text{norm},\infty}$ and TI$_\infty$ cluster (Fig. \ref{fig:Validation_WF_eff}).
\begin{figure*}[t!]
\centerline{\includegraphics[width=\textwidth]{Validation_WF_eff.pdf}}
 \caption{Directional and clustered wind farm efficiency for all the $U_{\text{norm},\infty}$ and TI$_\infty$ clusters defined in Sec.  \ref{sec:Calibration}). The shaded area represents the error on the mean with 95\% confidence level. Dashed lines refer to points rejected due to an error on the mean that exceeds 0.25.}\label{fig:Validation_WF_eff}
 \end{figure*}
The resulting directional and clustered wind farm efficiency exhibits, as expected, significant variability with lower values occurring for low $U_{\text{norm},\infty}$ (i.e. high $c_t$), low TI$_\infty$ (i.e. slower wake recovery), and winds aligned with the turbine arrays (i.e. W-E direction).
This analysis corroborates the high accuracy of the model in reproducing the wind farm flow at the turbine level. The Gauss model shows comparatively a poorer agreement with the experimental wind farm efficiency, especially for heavily waked sectors ($\theta_w \sim 90^\circ$ and $\theta_w \sim 270^\circ$), a general underestimation of losses for southerly winds, and a complete lack of wind farm overpower (i.e. $\eta>1$), which on the contrary is observed for the SCADA and the P2D-RANS for narrow wind sectors around $\theta_w=65^\circ$ and $\theta_w=245^\circ$.

In this regard, two occurrences of wind farm efficiency greater than one,
which are circled in green in Fig. \ref{fig:Validation_WF_eff}, are further investigated by extracting the individual turbine power gain (i.e. $-\Delta P_i$, which is more suitable to investigate speed-up conditions). Fig. \ref{fig:Validation_Speedup} displays $-\Delta P_i$ and the non-dimensional wind speed anomaly, $\hat{u}-1$, obtained from the P2D-RANS simulation corresponding to the bin centroids. The visualizations of the wind speed reveal the presence of channeling effects leading to power increases up to $\sim 10\%$, which are also observed in the SCADA data. This is further evidence of the importance of the rotor-induced pressure field on the performance of utility-scale wind turbines.
 \begin{figure*}[t!]
\centerline{\includegraphics[width=0.85\textwidth]{Validation_Speedup.pdf}}
 \caption{Power gain for selected wind direction bins within the cluster $U_{\text{norm},\infty}\in [0.7, 0.85)$, TI$_\infty \in[7,13.5)\%$: (a) SCADA power gain for $\theta_w\in [62.5,67.5)^\circ$; (b) P2D-RANS power gain for $\theta_w\in (62.5,67.5)^\circ$; (c) SCADA power gain for $\theta_w\in [242.5,247.5)^\circ$; (d) P2D-RANS power gain for $\theta_w\in (242.5,247.5)^\circ$. The velocity fields in (b) and (d) corresponds to the P2D-RANS solution for the bin centroid.} \label{fig:Validation_Speedup}
 \end{figure*}

The influence of the incoming flow on the wind farm performance is further characterized by disregarding power capture variability with the wind direction and calculating clustered wind farm efficiency and annual energy losses associated with the different $U_{\text{norm},\infty}$ and TI$_\infty$ bins (Fig. \ref{fig:Validation_Energy}).
 \begin{figure*}[t!]
\centerline{\includegraphics[width=\textwidth]{Validation_Energy.pdf}}
 \caption{Clustered wind farm efficiency (a-d) and annual energy losses (e-f) for difference models: (a) and (e) P2D-RANS; (b) and (f) Jensen; (c) and (g) Multizone; (d) and (h) Gauss. The dots indicate the SCADA statistics, the empty circles refer to the models. All points displayed satisfy the requirement of less than 25\% of error on the mean with a 95\% confidence.}\label{fig:Validation_Energy}
 \end{figure*}
The omnidirectional wind farm efficiency exhibits a neat variability for the different TI$_\infty$ clusters, with higher efficiency achieved for convective conditions, and the expected trend as a function of $U_{\text{norm},\infty}$, as higher losses are observed in region II of the power curve. From the modeling standpoint, the P2D-RANS shows a general slight underestimation of wake losses in stable atmospheric conditions and overestimation for neutral and convective clusters. Nevertheless, the P2D-RANS achieves better agreement in clustered wind farm efficiency (Fig. \ref{fig:Validation_Energy}(a)) and annual energy losses (Fig. \ref{fig:Validation_Energy}(e)) than the other wake models considered for this study.

The Jensen model (Fig. \ref{fig:Validation_Energy} (b, f)) generally overestimates the wake losses and does not capture the change in wind farm efficiency due to atmospheric stability. The difference in wind farm efficiency among different TI$_\infty$ values is solely due to intra-cluster variability of wind direction, being the analytical wake-field insensitive to the turbulence intensity due to the constant wake expansion coefficient $k_e=0.05$. The Multizone model (Fig. \ref{fig:Validation_Energy} (c, g)) also fails to fully capture the effect of the atmospheric stability, since it still uses a fixed wake expansion coefficient $k_e=0.05$. The overestimation of wake losses is however assuaged thanks to the improved formulation of the velocity field, compared to Jensen. The Gaussian model with the TI-dependent wake expansion \cite{Niayifar2015} $k_e=0.004 + 0.38$TI$_\infty$ shows a significant improvement, although it generally underestimates the wake losses in region II of the power curve. The performance of the Gaussian model is remarkable, considering the extremely low computational time (0.1 s per simulation) and the lack of a site-specific calibration of the wake expansion coefficient.

\begin{figure*}[t!]
\centerline{\includegraphics[width=\textwidth]{Validation_NMAE.pdf}}
 \caption{Normalized mean absolute error (NMAE) for the different performance indicator of the P2D-RANS and FLORIS wake models: (a) power capture; (b) single-turbine directional percentage power losses; (c) directional and clustered wind farm efficiency; (d) clustered annual energy losses.}\label{fig:Validation_NMAE}
\end{figure*}
The NMAE for all the performance indicators presented in this section are summarized in Fig. \ref{fig:Validation_NMAE} and for all the models. Although the P2D-RANS shows the smallest error for all the statistics, there are important distinctions that it is worth pointing out. Table \ref{tab:DNMAE} provides the percentage improvement attained by the P2D-RANS compared to the other wake models to assist the interpretation of the results. The P2D-RANS achieves a significant reduction of error for the farm and single-turbine power (especially waked turbines) compared to Jensen and Multizone, while the improvement compared to the more accurate TI-dependent Gaussian model is more limited (8\% for the single turbine power). Better relative performances are observed for the directional single-turbine power losses, with a 22\% error reduction achieved by the P2D-RANS compared to the Gaussian model, a consequence of the more accurate description of the wake velocity field and the modeling of the rotor-induced pressure field and the associated speed-up conditions for specific wind directions. A further significant margin of improvement is obtained for the statistics taking into account atmospheric stability through the clustering in TI$_\infty$, which is a clear indication that the site-specific calibration of $\nu_T$ is beneficial. The NMAE reduction with respect to the Gaussian model is 24\% for the directional and clustered wind farm efficiency and 40\%  for the annual energy losses. This last parameter is particularly relevant for the quantification of the financial losses due to wake interactions at a specific site. Therefore, the significantly enhanced accuracy in the prediction of the energy losses achieved by the P2D-RANS compared to other models is actually required to justify the investment in such LiDAR-driven CFD approach, which is computationally and logistically more expensive than the other engineering wake models.
\begin{table}[]
\centering
\caption{Percentage improvement in NMAE of the P2D-RANS model compared to the FLORIS wake models.}\label{tab:DNMAE}
\begin{tabular}{lccc}
 & Jensen & Multizone & Gauss \\ \hline
$P_\text{norm}$ (farm)                       & 37     & 27        & 16    \\
$P_\text{norm}$ (single-turbine)             & 23     & 15        & 8     \\
$P_\text{norm}$ (single-turbine, unwaked)    & 8      & 8         & 7     \\
$P_\text{norm}$ (single-turbine, waked)      & 29     & 18        & 8     \\
$\Delta P_i(\theta_w)$                       & 51     & 33        & 22     \\
$\eta(U_{\text{norm},\infty}$,TI$_\infty,\theta_w)$ & 57     & 39        & 24    \\
$\Delta E(U_{\text{norm},\infty}$,TI$_\infty)$      & 82     & 61        & 40   \\
\end{tabular}
\end{table}

For the sake of completeness, the gross and net capacity factors and the overall wake losses for the full SCADA dataset are reported in Table \ref{tab:CapacityFactor}. The P2D-RANS is still the model giving the best prediction, nevertheless its capability in simulating the wake interaction under different atmospheric stability and operative conditions is not clearly identifiable using these global statistical parameters.
\begin{table}[]
\centering
\caption{Gross capacity factor, net capacity factor, and wake loss factor estimated with the various wind farm models.}\label{tab:CapacityFactor}
\begin{tabular}{lccccc}
&SCADA  & P2D-RANS & Jensen & Multizone & Gauss  \\ \hline
Gross capacity factor & 0.801 & 0.800   & 0.799 & 0.799    & 0.800 \\
Net capacity factor   & 0.77 & 0.768   & 0.746 & 0.763    & 0.774 \\
Wake losses {[}\%{]}  & 3.87   & 4.00     & 6.63   & 4.51      & 3.25   \\
\end{tabular}
\end{table}

\section{Conclusions}\label{sec:Conclusions}
The theoretical fundamentals, numerical scheme, verification, and validation process of the Pseudo 2D-RANS (P2D-RANS) model have been discussed in this manuscript. The P2D-RANS uses a shallow-water approximation of the RANS equations to achieve an efficient solution of a wind farm flow vertically averaged over the rotor heights by inserting into the RANS solver appropriate corrective source terms, which take into account the contribution to the mass and momentum budgets of vertical fluxes and spatial dispersive stresses.

The set of depth-averaged and corrected flow equations is solved through an iterative marching scheme, which alternates sweeps in the streamwise directions for the solution of the momentum equations and global pressure corrections to ensure a solenoidal velocity field upon the convergence of the algorithm. The detailed derivation of the pressure correction, which may represent a generalizable and convenient numerical tool for the solution of partially-parabolized flows with source terms is also provided.

Clustered and statistically analyzed LiDAR measurements of wakes generated by utility-scale wind turbines have been leveraged to calibrate the turbulence closure and the actuator disk models implemented into the P2D-RANS. The optimally tuned turbulent eddy-viscosity exhibits physically reasonable trends as a function of the turbulence intensity of the incoming wind with a striking wake recovery enhancement detected for convective conditions, which might be connected with the occurrence of wake meandering. The thrust-force distribution over the rotor radial direction has been estimated for the various bins of incoming wind speed and turbulence intensity  at hub height. This data-driven approach has enabled the detection of a general reduction of the thrust force when turbine operations transition from region II to region III of the power curve

The correctness of the conceptual and computational framework of the P2D-RANS and its numerical implementation has been verified versus the solution of an elliptic axisymmetric RANS solver for single and in-tandem wind turbines. The validation process of the P2D-RANS has been broken down into three phases. First, effects of the daily cycle of the atmospheric stability on operations of an onshore farm have been simulated with a 10-minute resolution, then assessed against the SCADA data at the single-turbine level. Subsequently, the 10-minute-averaged power has been simulated for the whole set of available SCADA data to perform statistically meaningful error analysis. The power of single turbines shows an $R^2=0.956$, a normalized mean absolute error of $5.6\%$, and negligible bias, which represents an 8\% error reduction compared to the Gaussian wake model. Finally, wake losses have been quantified through different indicators and extensively analyzed as a function of the inflow conditions. The P2D-RANS reproduces satisfactorily the single-turbine directional wake losses, including the effect of a local increase of the freestream wind speed, referred to as flow speed-up conditions, which are induced by the rotor pressure field for specific wind directions. Such speed-up conditions have also repercussions on the wind farm efficiency, which is predicted with a normalized mean absolute error of 3\%, 24\% smaller than for the Gaussian model.

The P2D-RANS reproduces convincingly the annual energy wake losses (12\% error) for different stability and operative conditions, which indicates the financial relevance of LiDAR-driven CFD modeling in the context of the optimal design of wind farms. It is also proven the necessity for correct modeling of the variability of wind farm efficiency with atmospheric stability.

The current implementation of the P2D-RANS can solve the steady flow over a wind farm with 25 turbines and extending $30D$ by $30D$ in an averaged time of $\sim80$ s on a commercial laptop equipped with a single i5 1.70GHz core, which can enable performing real-time wind farm monitoring and control. The required computational time has been observed to be roughly half when operating an Intel Xeon Platinum 8160 48-core 2.1GHz cluster, which implies that the full characterization of the performances of a mid-sized farm can be completed in $\sim170$ CPU-hours.

\section*{Acknowledgments}
This material is based upon work supported by the National Science Foundation (grant no. IIP-1362022; Collaborative Research – I/UCRC for Wind Energy, Science, Technology, and Research) and from the WindSTAR I/UCRC Members of Aquanis, Inc., EDP Renewables, Bachmann Electronic Corp., GE Energy, Huntsman, Hexion, Leeward Asset Management, LLC, Pattern Energy, EPRI, LMWind, Texas Wind Tower, and TPI Composites. This research has been funded by the National Science Foundation CBET Fluid Dynamics (grant no. 1705837).

\section*{Conflict of interest declaration}
The authors have no conflicts to disclose.

redpajama_set_name: RedPajamaArXiv


Increased government scrutiny of the safety of the nation's trucking industry is contributing to a shortage of long-haul truck drivers. But it's not resulting in more jobs for some of the nearly 13 million unemployed Americans who can't find employment, USA Today reports.
"It's getting harder to get drivers," Mike Card, president of Combined Transport of Central Point, Ore., and incoming chairman of the American Trucking Associations, tells the newspaper. "I could hire 50 guys right now," says Card, who currently employs nearly 400 truck drivers.
A rise in the number of retiring baby boomers is contributing to the dwindling numbers of commercial truck drivers, the report says.
Young people who aren't interested in long-haul-trucking careers, which often require them be away from home for weeks at a time.
The training course is expensive for unemployed workers. Many unemployed construction and factory workers simply can't afford the $4,000 to $6,000 cost of a six-week driver-training course to obtain a commercial driver's license and meet individual states' regulatory requirements.
Truck drivers must be at least 21, leading many 18-year-olds who might consider trucking as a profession to pursue other trades, such as plumbing.
Trucking companies raised their standards. Recent publication of trucking companies' safety records by the federal government have led employers to hire only drivers with unblemished records, further reducing the applicant pool.
CareerBuilder says it has more than 35,000 job openings for truck drivers in its database. The job-search portal notes that the average national salary for a licensed commercial trucker is $44,500, and the top cities hiring for these positions are Chicago, Dallas and Columbus, Ohio.
The shortage of qualified drivers has prompted some trucking firms to poach truckers from other companies.
Competition among trucking companies has become so fierce that drivers are continually lured by better offers, the Los Angeles Times reports.
Job turnover rates for drivers at large, interstate fleets rose 2 percent to 90 percent during the first three months of the year, the highest in four years, according to the Trucking Activity Report, published by the American Trucking Associations.
The report notes that the rate of increase is even higher among smaller fleet operators -- up 16 percent to 71 percent.
Still, some within the industry question whether there is indeed a shortage of drivers. Rather, they say, truckers are being pickier about the jobs they are taking.
Writing recently on the Truck News Blog, a contributor named Stephen Large says he has to pay good drivers $35 to 40 an hour. "And it has to be local work, or they are not interested."
He further blames "too many ridiculous rules and regulations in the industry, and shippers and receivers who treat truck drivers poorly."
The solution to the problem is common sense, wrote another commenter, named Desiree Wood. In part, she said, "Pay a living wage [and] do not treat drivers like second class citizens."


How can you see wild horses in Montana or Wyoming? A trip with Pryor Wild ensures a fascinating visit with people who love the wild mustangs.
Las Vegas can be expensive, but there are also lots of things to do in Vegas for free. From free Vegas shows, to fun free things to do on the Las Vegas Strip, to free walking tours, and more.


Devote the vow you made during your wedding ceremony with this stunning wedding band brought to you by the prestigious brand, CrownRing. 14k white gold decorate this ring beautifully as a symbol of true love. This wedding band with model number WB-9604W features Round-cut stones intimately carved on a Pave setting.
Diamonds on Wabash is an authorized distributor of CrownRing in Chicago, IL. If you need more details about this CrownRing wedding band, feel free call us by phone or email, and our professional jewelers will be more than happy to guide you.


text:
Theodore Maxant Conservation Restriction
In 2016, Theodore "Ted" Maxant approached the Harvard Conservation Trust (HCT) with an idea for how he could help preserve the rural character of Still River Village. Born and raised in Still River, it is place Ted knows by heart. In 2007, his father, William T. Maxant, conserved the family homestead through the gift of a conservation restriction (CR) to HCT. William also took steps to ensure the property would someday benefit the Town's Conservation Commission and Historical Society. Ted's idea was to build on his family's conservation legacy by purchasing a key abutting parcel of land, and then extinguishing the development rights to that land through the gift of a new CR to HCT. As Ted explained, this would preserve the scenic and open character of this part of Still River Village, as well as create the possibility for a having the Historical Society property serve as a gateway into the extensive network of conservation lands that lie beyond (see map). With trail linkages being an important part of both HCT's strategic land protection efforts as well as the Town's recently updated Open Space & Recreation Plan, HCT requested that a public access clause be included in the CR. Thankfully, Ted agreed. Given that land put under a CR by a private landowner generally remains private property, the right to public trail access was an exceptional addition to this already generous conservation gift. The CR gift was completed in December, and HCT is looking forward to working with Ted to realize his vision for a connecting trail.

redpajama_set_name: RedPajamaCommonCrawl


i have a problem in running my sap system.when i start it ,the message server and gateway are starting,but there is some problem with my disp+work .its not getting started up.
Could you plz paste more info from the dev_disp.
Can youp lease post the log file: dev_w0, dev_disp is a higher level log and generally doesn't give enough information about why the dispatcher died.
Sorry couldn't figure out the cause for the ERROR message in dev_disp . Did you check the dev_rd ( the trace file for gateway) ? do you see any error there ?