1. Introduction
Wind turbine wakes are central to power loss and flow interactions across wind farms. Yet accurately and efficiently capturing their unsteady behaviour remains a major challenge in wind energy research. Conventional engineering wake models, such as the widely used Jensen model (Jensen Reference Jensen1983), and its extensions (Katic, Højstrup & Jensen Reference Katic, Højstrup and Jensen1987; Frandsen et al. Reference Frandsen, Barthelmie, Pryor, Rathmann, Larsen, Hojstrup and Thogersen2006), typically assume simple top-hat velocity deficit profiles for steady-state power estimation. More recent models have adopted more realistic profiles, such as Gaussian (Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2013) or cosine (Tian et al. Reference Tian, Zhu, Shen, Zhao and Shen2015) distributions, but they still overlook key unsteady features such as meandering, intermittency and coherent structures. These dynamics significantly affect power fluctuations and structural loads, and are critical for instantaneous power estimation, dynamic control design, and grid frequency regulation (Milan, Wächter & Peinke Reference Milan, Wächter and Peinke2013; Stevens & Meneveau Reference Stevens and Meneveau2017). Efficient models capable of resolving the full spatiotemporal evolution of wakes are thus vital for next-generation wind farm design and operation.
In this context, the dynamic wake meandering (DWM) model was proposed (Larsen et al. Reference Larsen, Madsen, Thomsen and Larsen2008), offering a framework to incorporate unsteady effects in engineering-focused simulations. This model treats the wake as a passive tracer whose large-scale motion is driven by atmospheric boundary layer (ABL) turbulence, thus capturing the temporal displacement of the wake. The physical foundation of DWM traces back to the early hypothesis by Lissaman, Zambrano & Gyatt (Reference Lissaman, Zambrano and Gyatt1983), which identified large-scale atmospheric eddies as the dominant driver of wake meandering – a view supported by wind tunnel experiments and field measurements (Espana et al. Reference Espana, Aubrun, Loyer and Devinant2011; Trujillo et al. Reference Trujillo, Bingöl, Larsen, Mann and Kühn2011; Brugger, Markfort & Porté-Agel Reference Brugger, Markfort and Porté-Agel2022). Over time, the DWM model has been extended to incorporate more physical effects, including ABL shear (Keck et al. Reference Keck, De Maré, Churchfield, Lee, Larsen and Madsen2015), yaw misalignment (Branlard, Martínez-Tossas & Jonkman Reference Branlard, Martínez-Tossas and Jonkman2023) and wake-added turbulence (WAT) modelling (Branlard et al. Reference Branlard, Jonkman, Platt, Thedin, Martínez-Tossas and Kretschmer2024). Model parameters have also been systematically calibrated using experimental and numerical approaches (Madsen et al. Reference Madsen, Larsen, Larsen, Troldborg and Mikkelsen2010; Doubrawa Moreira et al. Reference Doubrawa Moreira, Annoni, Jonkman and Ghate2018; Yu et al. Reference Yu, Euzenat, Maheshwari, Martinez, Eneau and Pallud2024), and different DWM variants have been compared and evaluated (Rivera-Arreba et al. Reference Rivera-Arreba, Eliassen, Bachynski and Panjwani2021; Hanssen-Bauer et al. Reference Hanssen-Bauer, Doubrawa, Madsen, Asmuth, Jonkman, Larsen, Ivanell and Stenbro2023). To avoid the complexity of coupling with external aero-elastic solvers, statistical wake meandering models based on wake centreline probability density functions (PDFs) have been proposed (Thøgersen et al. Reference Thøgersen, Tranberg, Herp and Greiner2017; Braunbehrens & Segalini Reference Braunbehrens and Segalini2019; Brugger et al. Reference Brugger, Markfort and Porté-Agel2024). These mid-fidelity models have demonstrated reliable accuracy in far-wake predictions when validated against large-eddy simulations (LES) results (Jonkman et al. Reference Jonkman, Doubrawa, Hamilton, Annoni and Fleming2018; Shaler & Jonkman Reference Shaler and Jonkman2021; Rivera-Arreba et al. Reference Rivera-Arreba, Li, Yang and Bachynski-Polić2024) and field data (Shaler, Debnath & Jonkman Reference Shaler, Debnath and Jonkman2020; Kretschmer et al. Reference Kretschmer, Jonkman, Pettas and Cheng2021; Shaler et al. Reference Shaler, Quon, Ivanov and Jonkman2024), making them suitable for load estimation and transient response studies (Wise & Bachynski Reference Wise and Bachynski2020).
Despite these advances, the physical origin of wake meandering remains debated. While the classic DWM formulation attributes wake displacement to large-scale turbulence, an alternative perspective highlights the role of intrinsic shear-layer instabilities, as evidenced by laboratory experiments (Medici & Alfredsson Reference Medici and Alfredsson2006; Okulov et al. Reference Okulov, Naumov, Mikkelsen, Kabardin and Sørensen2014; Howard et al. Reference Howard, Singh, Sotiropoulos and Guala2015; Foti et al. Reference Foti, Yang, Campagnolo, Maniaci and Sotiropoulos2018) and numerical studies (Kang, Yang & Sotiropoulos Reference Kang, Yang and Sotiropoulos2014; Li, Dong & Yang Reference Li, Dong and Yang2022). Recent work suggests a coexistence of both mechanisms (Heisel, Hong & Guala Reference Heisel, Hong and Guala2018; Yang & Sotiropoulos Reference Yang and Sotiropoulos2019b ), where large-scale turbulence amplifies inherent wake instabilities (Coudou, Buckingham & van Beeck Reference Coudou, Buckingham and van Beeck2017). Building on this dual-mechanism hypothesis, subsequent extensions to the DWM framework have incorporated both effects by selectively amplifying turbulent fluctuations, leading to improved predictions of WAT intensity (Feng et al. Reference Feng, Gupta, Li and Wan2024a ).
Nevertheless, most existing DWM models still lack accurate treatment of near-wake dynamics. Like traditional steady-state models, they often prioritise far-wake accuracy at the expense of physical consistency in the near wake. Major deficiencies include the omission of pressure recovery and the artificial prescription of wake expansion onset, leading to physically incomplete calculations. The role of pressure recovery is in fact substantial in the near wake (Liew, Heck & Howland Reference Liew, Heck and Howland2024), and has recently begun to receive attention in wake modelling (Ge et al. Reference Ge, Wu, Liu and Yang2019). For simplicity, steady-state models have adopted self-similar double-Gaussian (DG) velocity deficit profiling (Keane et al. Reference Keane, Aguirre, Ferchland, Clive and Gallacher2016) to improve near-wake performance, with extensions that enhance mass conservation (Schreiber, Balbaa & Bottasso Reference Schreiber, Balbaa and Bottasso2020; Keane Reference Keane2021) and validate performance against LES (Soesanto, Yoshinaga & Iida Reference Soesanto, Yoshinaga and Iida2022). More recently, data-driven variants of DG models (Wang et al. Reference Wang, Chen and Chen2024a ) have achieved improved agreement with high-fidelity data through refined profile fitting. However, these approaches remain fundamentally limited: they treat near-wake modelling as a curve-fitting task, without addressing the deeper issue of how the underlying governing equations should be closed. Consequently, while they may reproduce mean velocity profiles, they lack dynamical consistency – particularly under unsteady or complex inflow conditions. More importantly, the fitted profiles offer limited physical insight. Therefore, a salient open question is how to achieve equation-level closure of DWM that not only improves wake flow prediction but also preserves physical interpretability.
Symbolic regression (SR) offers a promising solution to this modelling gap. As an interpretable machine learning technique, SR can automatically discover concise and accurate mathematical expressions that describe complex nonlinear relationships in data (Angelis, Sofos & Karakasidis Reference Angelis, Sofos and Karakasidis2023). Compared to black-box data-driven methods such as neural networks, SR is better suited for explicitly uncovering governing physical laws (Makke & Chawla Reference Makke and Chawla2024). Mainstream approaches include sparse identification of nonlinear dynamics (SINDy) (Brunton, Proctor & Kutz Reference Brunton, Proctor and Kutz2016) and genetic programming with expression trees (Chen et al. Reference Chen, Luo, Liu, Xu and Zhang2022; Cranmer Reference Cranmer2023), which have shown strong performance across applications (La Cava et al. Reference La Cava, Orzechowski, Burlacu, de Francca, Virgolin, Jin, Kommenda and Moore2021). Recent advances in SR algorithms (Kamienny et al. Reference Kamienny, d’Ascoli, Lample and Charton2022; Tenachi, Ibata & Diakogiannis Reference Tenachi, Ibata and Diakogiannis2023; Du, Chen & Zhang Reference Du, Chen and Zhang2024) have further accelerated its adoption in data-driven physical modelling (Gamella, Peters & Bühlmann Reference Gamella, Peters and Bühlmann2025). In fluid dynamics, SR has recently demonstrated success in modelling flow behaviour with both accuracy and interpretability (Ma et al. Reference Ma, Zhang, Feng and Wen2024; Wang et al. Reference Wang, Chen and Chen2024a ).
To address the limitations of current models, we propose a symbolic regression-enhanced dynamic wake meandering (SRDWM) framework. The first contribution is the integration of symbolic expressions, extracted from LES data, into the DWM governing equations to reconstruct volumetric forcing and boundary terms – achieving equation-level closure and restoring near-wake dynamics. Second, we introduce a hierarchical SR strategy guided by domain knowledge to reduce dimensionality and improve generalisability. Additionally, SRDWM incorporates a revised WAT formulation to enhance wake turbulence predictions under the dual-mechanism hypothesis of wake meandering. Together, these developments yield a model that recovers physical consistency in both space and time, capturing unsteady wake evolution with high accuracy and computational efficiency.
The rest of the paper is structured as follows. Section 2 details the inflow generation, wake simulation, and the proposed SRDWM framework. Section 3 presents the discovered symbolic terms, and evaluates wake dynamics across inflow conditions. Finally, we draw a conclusion in § 4.
2. Methods
This section introduces the SRDWM framework (see figure 1), which enhances DWM by embedding symbolic expressions – reconstructed from LES data via hierarchical SR – into the governing equations. These replace empirical terms for forcing and boundary conditions, improving physical consistency. A revised WAT model further refines wake turbulence predictions. Together, these components yield an accurate, efficient and interpretable wake model with full spatiotemporal resolution.
Figure 1. Schematic of the SRDWM modelling framework. The proposed method (black) augments the DWM model (blue) with hierarchical SR (red) to reconstruct volumetric forcing and wake boundary conditions from data of LES (green), and improves predictions of wake dynamics with the WAT model (yellow).
2.1. Inflow turbulence generation
For LES and DWM-type models, turbulence characteristics in the incoming flow strongly influence wake recovery and meandering behaviour (Yang & Sotiropoulos Reference Yang and Sotiropoulos2019a
). In this study, the inflow turbulence, serving as the inlet boundary condition for turbine wakes, is generated using two approaches: precursor simulation (PS) and synthetic turbulence (ST). Both methods produce ABL turbulence with identical logarithmic mean streamwise velocity profiles, but differ in turbulence intensity. Specifically, mean velocity at hub height is set to
$U_{\infty }=8 $
m s−1 in all inflow cases. The PS case exhibits hub-height turbulence intensity 6.1 %, while the ST cases span a broader range: 1.8 % (ST1), 4.1 % (ST2), 6.4 % (ST3), 8.4 % (ST4) and 10.3 % (ST5).
In the PS approach, a fully developed turbulent ABL flow is simulated using LES with periodic horizontal boundary conditions. Once statistical steady state is reached, the velocity fields at the inlet are recorded at each time step for use in subsequent turbine wake simulations.
Alternatively, the ST method generates inflow turbulence by modelling a stationary stochastic process governed by prescribed cross-spectral density. Initially, the Kaimal spectrum (Kaimal et al. Reference Kaimal, Wyngaard, Izumi and Coté1972) is applied to define the turbulence energy distribution in the frequency domain, ensuring that the resulting velocity time series adheres to the target power spectral density (PSD). The turbulence spectrum, following the IEC 61400-1 standard, is expressed as
\begin{equation} \frac {f S_{\chi }(f)}{\sigma _{\chi }^{2}} = \frac {4 f L_{\chi } / \bar {U}_{\textit{hub}}}{(1 + 6 f L_{\chi } / \bar {U}_{\textit{hub}})^{5/3}}, \end{equation}
where
$S_{\chi }$
denotes the auto-spectral density function,
$\chi \in \{x,y,z\}$
represents directions,
$\sigma _{\chi }$
is the standard deviation of velocity components,
$L_{\chi }$
is the integral length scale,
$f$
is the cyclic frequency, and
$\bar {U}_{\textit{hub}}$
is the mean hub-height velocity. The velocity fluctuation components are related by
$\sigma _x = 1.25\sigma _y = 2\sigma _z$
. To further account for spatial coherence, a model is incorporated to characterise the correlation properties of turbulent fluctuations:
\begin{equation} \varTheta (r, f) = \exp \left ( -12 \sqrt {\left (\frac {f r}{\bar {U}_{\textit{hub}}}\right )^2 + \left (\frac {0.12 r}{L_c}\right )^2 } \right )\!, \end{equation}
where
$\varTheta (r,f)$
quantifies coherence between two nodes separated by distance
$r$
(perpendicular to wind direction), and
$L_c$
is the coherence length scale. Synthesised turbulent fluctuations are subsequently superimposed onto the mean velocity profile, forming spatiotemporally correlated wind fields. The resulting turbulence data undergo inverse Fourier transformation and are structured as a time series along a two-dimensional boundary, establishing the inflow turbulence conditions.
Figure 2 summarises the statistics of all inflows by showing the time-averaged streamwise velocity and the turbulent fluctuation profiles as functions of height. The mean velocity profiles in all cases follow the logarithmic law and remain nearly identical. For the turbulent fluctuations, the PS case, derived from a fully developed ABL flow, exhibits good agreement with higher-resolution high-fidelity simulations. In contrast, for all ST cases, the synthetic turbulent fluctuations are prescribed to be independent of height. These profiles define the inflow conditions and provide necessary context for interpreting case-to-case differences in subsequent wake behaviour.
Figure 2. The normalised horizontal- and time-averaged (a) streamwise velocity and velocity variances of the (b) streamwise, (c) spanwise, and (d) vertical components for all turbulent inflows. Cyan dashed and brown solid lines denote high-resolution LES results from Stevens, Wilczek & Meneveau (Reference Stevens, Wilczek and Meneveau2014) and Feng et al. (Reference Feng, Gupta, Li and Wan2024b ), respectively.
2.2. Turbine wake simulation
Utilising the numerical solver SOWFA (Churchfield et al. Reference Churchfield, Lee, Michalakes and Moriarty2012), LES are performed to simulate wind turbine wakes under varying inflow turbulence intensities. The turbine effect is represented using the actuator disk method (ADM) with rotation (ADM-R) without nacelle and tower influence. Prior studies (Porté-Agel et al. Reference Porté-Agel, Wu, Lu and Conzemius2011; Wu & Porté-Agel Reference Wu and Porté-Agel2011) demonstrated that ADM-R and the more complex actuator line method (ALM) show good agreement with wind tunnel data for both time-averaged velocity and turbulence intensity, even in the near wake at
$x = 2D$
, and significantly outperform the standard non-rotating ADM. Additionally, Qian, Song & Ishihara (Reference Qian, Song and Ishihara2022) reported minimal differences in the power spectrum of streamwise velocity between ADM-R and ALM in the near-wake region (
$x \lt 1D$
), except for high-frequency content at
$x=0.1D$
. The wake dynamics is governed by the filtered incompressible Navier–Stokes (NS) equations:
\begin{align} \frac {\partial {\widetilde {U}_j}}{\partial {x_j}}=0,\\[-28pt] \nonumber \end{align}
\begin{align} \frac {\partial \widetilde {U}_i}{\partial t}+\frac {\partial (\widetilde {U}_i\widetilde {U}_j)}{\partial {x_j}}=-\frac {\partial \widetilde {p}^*}{\partial x_i}-\frac {\partial \tau _{i j}^D}{\partial x_j}+\frac {1}{\rho } F_i^T. \\[2pt] \nonumber \end{align}
Here,
$\widetilde {U}_{i}$
denotes the filtered velocity component, where
$i = 1, 2, 3$
correspond to the streamwise (
$x$
), spanwise (
$y$
) and vertical (
$z$
) directions, respectively. The term
$\widetilde {p}^*$
represents the modified filtered kinematic pressure, and
$\tau _{\textit{ij}}^D$
is the deviatoric component of the sub-grid scale (SGS) stress. The SGS effects are parametrised using the Lagrangian-averaged scale-dependent dynamic model (Bou-Zeid, Meneveau & Parlange Reference Bou-Zeid, Meneveau and Parlange2005), which dynamically computes anisotropic eddy viscosity with scale-adaptive Smagorinsky coefficients. The turbine force term
$F_i^T$
is derived by spatially smoothing the aerodynamic blade forces via a convolution operation:
\begin{equation} F_i^T = \frac {\boldsymbol{f}}{\epsilon ^3 \pi ^{3/2} \,{\rm d}V} \exp \left (-\left (\frac {d}{\epsilon }\right )^2 \right )\!, \end{equation}
where
${\rm d}V$
is the volume element,
$\epsilon = 20$
m controls the smoothing effect, and
$d$
represents the distance between the actuator point and the computational grid centre. This regularisation ensures numerical stability and a smoothly distributed turbine force. The aerodynamic force
$\boldsymbol{f}$
, acting on the blade elements, is computed using blade element momentum theory:
\begin{equation} \boldsymbol{f}=\frac {1}{2} \rho U_{{rel}}^2 N_b c\, {\rm d}r\, (C_L \boldsymbol e_L+C_D \boldsymbol e_D), \end{equation}
where
${U}_{{rel}}$
represents the relative velocity,
$N_b$
is the number of blades,
$c$
is the blade-chord length, and
${\rm d}r$
is the radial width. The lift (
$C_L$
) and drag (
$C_D$
) coefficients are interpolated from aerodynamic look-up tables, while
$\boldsymbol e_L$
and
$\boldsymbol e_D$
denote unit vectors. The wind turbine modelled in this study is the NREL 5 MW reference turbine, characterised by a
$126$
m rotor diameter (D) and
$90$
m hub height. It operates at
$9$
rotations per minute, maintaining an optimal tip-speed ratio. Following prior studies (Liu & Stevens Reference Liu and Stevens2020; Yang et al. Reference Yang, Lang, Du, Jin, Li and Ge2022), the computational domain, spanning
$1800\times 1000\times 600$
m
$^3$
, is meshed using a structured grid comprising
$180\times 100\times 80$
nodes in the
$x$
-,
$y$
- and
$z$
-directions. Periodic boundary conditions are enforced in the spanwise directions, while a wall model resolves near-ground dynamics. Velocity at the fourth grid point is adjusted to correspond to the wall shear stress, reducing discrepancies in the logarithmic layer (Kawai & Larsson Reference Kawai and Larsson2012). The accuracy of the simulated flow fields has been validated against high-fidelity data (Chen et al. Reference Chen, Wang, Feng, Tian, Gupta, Cao, Wan and Chen2025), with further computational set-up and solver implementation detailed in Feng et al. (Reference Feng, Li, Gupta and Wan2022) and Wang et al. (Reference Wang, Feng, Peng, Mao, Doranehgard, Gupta, Li and Wan2023).
2.3. The SRDWM model
2.3.1. The DWM framework
The DWM model hypothesises that wake behaviour resembles that of a passive tracer advected by large-scale atmospheric turbulence (Larsen et al. Reference Larsen, Madsen, Thomsen and Larsen2008). Under this assumption, the model applies a thin shear layer approximation of the NS equations (Ainslie Reference Ainslie1988) to compute the velocity deficit in a meandering frame, where the deficit is assumed quasi-steady relative to the instantaneous wake centre. Wake meandering is then introduced by advecting the deficit pattern along the pathlines determined by the low-pass filtered velocity field.
The original governing equations neglect pressure gradient effects and assume dominance of radial over axial velocity gradients in the axisymmetric formulation. However, as shown in figure 3(a), LES results reveal that both mean pressure gradient (term VI) and mean turbine-induced forcing (term VII) make significant contributions, particularly in the near-wake region. We therefore restore these terms, leading to the following momentum and continuity equations:
\begin{align} \underbrace {r U_x \frac {\partial U_x}{\partial x}}_{\substack {\text{I} \\[2pt] \text{axial advection}}} \!+\! \underbrace {r U_r \frac {\partial U_x}{\partial r}}_{\substack {\text{II} \\[2pt] \text{radial advection}}} &\! = \underbrace {\nu _T \frac {\partial U_x}{\partial r} + r \nu _T \frac {\partial ^2 U_x}{\partial r^2} + r \frac {\partial \nu _T}{\partial r} \frac {\partial U_x}{\partial r}}_{\substack {\text{III}{-}\text{V} \\[2pt] \text{radial diffusion}}}\!-\!\! \underbrace {r \frac {\partial P}{\partial x}}_{\substack {\text{VI} \\[2pt] \text{pressure gradient}}}\! +\!\! \underbrace {r f_x}_{\substack {\text{VII} \\[2pt] \text{turbine forcing}}}\!\!\!, \\[-28pt] \nonumber \end{align}
\begin{align} U_r + r \frac {\partial U_r}{\partial r} + r \frac {\partial U_x}{\partial x} & = 0.\\[2pt] \nonumber \end{align}
Here,
$U_x$
and
$U_r$
denote the axial and radial velocity components,
$r$
stands for the radial position, and
$\nu _T$
represents the eddy viscosity. Following Shaler & Jonkman (Reference Shaler and Jonkman2021),
$\nu _T$
is modelled as the sum of ambient velocity and wake shear contributions:
\begin{align} \nu _T(x, r) = & {\mathcal{F}}_{\nu , \textit{Amb}}(x)\, k_{\nu , \textit{Amb}}\, {\widehat{TI}}_{\textit{Amb}}\, R\, {\widehat{U}}_{x,\textit{Amb}}\nonumber \\ & + {\mathcal{F}}_{\nu , {\textit{Shr}}}(x)\, k_{\nu , {\textit{Shr}}}\max \left ( R^2 \left | \frac {\partial U_x}{\partial r} \right |, R {{\min}_r} U_x \right )\!, \end{align}
where a hat denotes low-pass temporal filtering,
$R$
is the rotor radius, and subscript
$\textit{Amb}$
denotes quantities in ambient velocity field. The filter functions
${\mathcal{F}}$
account for the initial imbalance between velocity and turbulence fields:
\begin{equation} {\mathcal{F}}(\dot {x}) = \begin{cases} C_{\textit{min}}, & \dot {x} \leq \xi _{\textit{min}}, \\ C_{\textit{min}} + \left ( 1 - C_{\textit{min}} \right ) \left ( \dfrac {\dot {x} - \xi _{\textit{min}}} {\xi _{\textit{max}} - \xi_{\textit{min}}} \right )^{\lambda }, & \xi _{\textit{min}} \lt \dot {x} \lt \xi _{\textit{max}}, \\ 1, & \dot {x} \geq \xi _{\textit{max}} , \end{cases} \end{equation}
with parameters calibrated in a prior study (Doubrawa Moreira et al. Reference Doubrawa Moreira, Annoni, Jonkman and Ghate2018). The filter function thus represents the delayed build-up of turbulent stresses, with a fixed minimum value near the rotor, a unity value in the far wake, and a smooth transition in between. The coefficients
$k_{\nu ,\textit{Amb}}$
and
$k_{\nu ,{\textit{Shr}}}$
associated with dissipation are linearly regressed from the turbulence intensity
$TI$
, as shown in figure 3(b).
Figure 3. (a) Budget of the mean momentum equation normalised by
$U_\infty ^2$
for the LES PS case. (b) Linear fitting of coefficients
$k_{\nu , \textit{Amb}}$
and
$k_{\nu , {\textit{Shr}}}$
versus turbulence intensity
$TI$
.
The SRDWM model developed in this study is implemented within the FAST.Farm platform of the National Renewable Energy Laboratory (Shaler & Jonkman Reference Shaler and Jonkman2021). The wake field is discretised on a polar grid with radial spacing
$5$
m (Rivera-Arreba et al. Reference Rivera-Arreba, Li, Yang and Bachynski-Polić2024), and the governing equations are solved using a second-order finite difference scheme with the implicit Crank–Nicolson method. The instantaneous wake dynamics follow the passive wake meandering concept, in which large-scale lateral and vertical displacements of the wake plane are resolved. The displacement velocity of each wake plane is obtained by spatial averaging of the wake velocity field, with recommended characteristic length
$3D$
(Feng et al. Reference Feng, Gupta, Li and Wan2024a
), corresponding to a meandering cut-off frequency at Strouhal number
$St=0.3$
. The spatial averaging employs a windowed jinc function (first-order Bessel function), with weights determined by the radial distance from the wake centreline.
To improve the underestimation of turbulence intensity in the original DWM formulation, the WAT model (Madsen et al. Reference Madsen, Larsen, Larsen, Troldborg and Mikkelsen2010) is commonly employed. Building upon this, Feng et al. (Reference Feng, Gupta, Li and Wan2024a ) proposed a scale-dependent amplification mechanism based on shear instability, realising the dual-mechanism hypothesis. While this WAT model enhances spectral predictions, it also indiscriminately amplifies velocity fluctuations across the entire wake. As a result, it leads to overpredictions of near-wake turbulence intensity, which has been shown to exhibit localised suppression (Heisel et al. Reference Heisel, Hong and Guala2018; Li et al. Reference Li, Zhang, Li and Yang2024).
To resolve this issue, we propose a revised amplification scheme using a wake-scaling factor
$K_{\textit{WAT}}(x, r)$
that selectively enhances band-pass filtered fluctuations. These filtered components, denoted by
$\boldsymbol{u}^{\prime }_{{bp}}$
, are extracted through recursive exponential smoothing, combining low-pass
$\boldsymbol{u}_{{lp}}$
and high-pass
$\boldsymbol{u}^{\prime }_{{hp}}$
signals as follows:
\begin{equation} \begin{aligned} \boldsymbol{u}^{\prime }_{{bp}}[n+1] &= \beta \, \boldsymbol{u}^{\prime }_{{bp}}[n] + (1 - \beta )\, \boldsymbol{u}^{\prime }_{{hp}}[n+1], \\ \boldsymbol{u}^{\prime }_{{hp}}[n+1] &= \boldsymbol{u}[n] - \boldsymbol{u}_{{lp}}[n], \\ \boldsymbol{u}_{{lp}}[n+1] &= \alpha \, \boldsymbol{u}_{{lp}}[n] + (1 - \alpha )\, \boldsymbol{u}[n+1], \end{aligned} \end{equation}
where
$\boldsymbol{u}$
denotes instantaneous velocity, and
$\boldsymbol{u}^{\prime }$
denotes its fluctuation. The filtering coefficients are defined as
$\alpha = \textrm{e}^{-2\pi\, \Delta t\, f_1}$
and
$\beta = \textrm{e}^{-2\pi\, \Delta t\, f_2}$
, with cut-off frequencies
$f_1$
and
$f_2$
spanning the energy-dominant Strouhal (
$St=fD/U_{\infty }$
) band
$0.1 \lt St \lt 1.0$
(Feng et al. Reference Feng, Li, Gupta and Wan2022).
The refined wake-scaling factor is defined as
\begin{equation} K_{\textit{WAT}}(x, r) = \frac{\mathcal{F}_{s,\textit{Shr}}(x) \, k_{s,\textit{Shr}}}{U_{x,\textit{Amb}}} \, \frac{\partial U_x}{\partial r} - \mathcal{F}_{s,\textit{Amb}}(x) \, k_{s,\textit{Amb}} \left| 1 - \frac{U_x}{U_{x,\textit{Amb}}} \right| \end{equation}
with
$k_{s,\textit{Amb}}=3.2$
and
$k_{s,{\textit{Shr}}}=2.7$
. Unlike previous formulations, this structure – accounting for negative velocity gradients near the centreline – enhances shear-driven scaling while suppressing unphysical growth in the near wake, where the velocity deficit peaks. By amplifying coherent band-pass fluctuations at each time step, expressed as
$\boldsymbol{u}=\boldsymbol{u}+K_{\textit{WAT}}\boldsymbol{u}^{\prime }_{{bp}}$
, this approach yields more accurate turbulence intensity distributions, as will be demonstrated in § 3.
2.3.2. Symbolic regression
Symbolic regression is a data-driven approach for discovering interpretable equations that explicitly capture complex nonlinear relationships between variables, providing concise mathematical descriptions of physical system behaviour. It systematically explores combinations of variables, constants and mathematical operators, while balancing descriptive accuracy with model simplicity (Cranmer Reference Cranmer2023). In this study, we adopt the tree-based SR method, which evolves expression trees – composed of operators (nodes) and variables or constants (leaves) – using a multi-population genetic algorithm. The evolutionary process includes selection, crossover and mutation, as illustrated schematically in figure 4. Selection is governed by asynchronous tournament strategies (Goldberg & Deb Reference Goldberg and Deb1991), while inter-population migration enhances diversity and prevents premature convergence. Simulated annealing (Kirkpatrick, Gelatt & Mario Reference Kirkpatrick, Gelatt and Mario1983) is integrated to improve exploration, and the final symbolic expressions are refined through BFGS-based constant optimisation (Broyden Reference Broyden1970).
Figure 4. Schematic of evolution operations in the genetic algorithm-based SR.
Figure 5. Schematic of the hierarchical wake modelling strategy, where spanwise profiles are prescribed a priori (e.g. Gaussian) and their parameters (amplitude
$a$
, standard deviation
$\sigma$
) evolve along the streamwise direction.
This evolutionary SR approach employs a population-based search to identify the most suitable analytical expressions. The core of this process is the fitness criterion, which is designed to balance accuracy and simplicity, reflecting a key principle of scientific discovery. Accuracy is quantified by the mean squared error (MSE), while simplicity is measured by the expression’s complexity, defined as the total number of nodes in its expression tree. The fitness of an expression is evaluated as
\begin{equation} F=\frac {\sum _i^N \left (\mathscr{h}(x_i)-{\mathscr{Y}}_i \right )^2}{N} \times \exp {(f_c)}, \end{equation}
where
$\mathscr{h}(x)$
is the discovered expression,
${\mathscr{Y}}$
are the corresponding true values, and
$N$
is the number of data points. The term
$f_c$
is the frequency of expressions at a given complexity within a moving window, which ensures a balanced exploration of solutions with varying complexities across the population. The search ultimately yields a Pareto front of optimal expressions that represents the trade-off between accuracy and simplicity.
Given the limited availability of computational fluid dynamics (CFD) data in wake simulations, we adopt a hierarchical regression strategy guided by domain knowledge. Specifically, we assume that pressure and turbine forcing terms exhibit separable forms: the spanwise distribution is prescribed a priori using physically motivated functions (see (2.14) and (2.15) below), while the streamwise evolution is inferred through SR. This strategy, schematically illustrated in figure 5, is inspired by the well-established paradigm of traditional wake modelling (Jensen Reference Jensen1983; Katic et al. Reference Katic, Højstrup and Jensen1987; Frandsen et al. Reference Frandsen, Barthelmie, Pryor, Rathmann, Larsen, Hojstrup and Thogersen2006; Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2013; Keane et al. Reference Keane, Aguirre, Ferchland, Clive and Gallacher2016), which similarly assumes a self-similar spanwise profile. However, instead of using simplified conservation-based parameters to govern streamwise evolution, our method leverages SR to discover the governing relationships directly from data. To validate this foundational assumption, the LES results shown in figure 6 reveal that both pressure and turbine forcing display strong self-similarity in the spanwise direction, which can be compactly represented by Gaussian and DG profiles, respectively:
\begin{align} 1-\frac {P}{\rho U_{\infty }^2} &= a \times \exp \left (-\frac {\dot {y}^2}{2\sigma ^2}\right )\!,\\[-28pt] \nonumber \end{align}
\begin{align} -\frac {F_x D}{\rho U_{\infty }^2} & = a \times \left (\exp \left (-\frac {(\dot {y} - \mu )^2}{2\sigma ^2}\right ) + \exp \left (-\frac {(\dot {y} + \mu )^2}{2\sigma ^2}\right )\right )\!,\\[2pt] \nonumber \end{align}
where
$\dot {x} = x/D$
and
$\dot {y} = y/D$
denote normalised coordinates. The parameters
$a(\dot {x})$
,
$\mu (\dot {x})$
and
$\sigma (\dot {x})$
describe the streamwise variation and are treated as independent scalar functions to be identified. Notably, the turbine forcing is highly localised in space, acting only within a narrow region immediately downstream of the rotor, and it decays rapidly and becomes negligible within approximately
$0.3D$
. To reflect this physical behaviour, we adopt a piecewise formulation in which (2.15) applies only within the localised influence region, while the forcing is set to zero elsewhere. These formulations transform the original two-dimensional regression task into one-dimensional sub-problems, forming a hierarchical model structure. By decoupling the spatial dimensions based on physical insight, the SR process becomes more tractable, with significantly reduced searching space of SR, and the resulting expressions more interpretable.
Figure 6. Spanwise distribution of (a) mean pressure and (b) mean turbine force at various downstream positions in the PS case. Circles show LES results, while dash-dotted lines denote fitted Gaussian/DG profiles.
The normalised streamwise position
$\dot {x}$
from the CFD data serves as the SR input variable. For Gaussian/DG fitting, the dataset comprises
$125 \times 63$
samples spanning both the streamwise and spanwise directions. Accordingly, the dataset used for SR reduces to
$125 \times 3$
, focusing solely on the fitted parameters at
$125$
streamwise locations. To perform SR, we construct a candidate operator library consisting of basic arithmetic operations, powers, exponentials and cosine functions – operators commonly encountered in wake modelling. To ensure parsimony and prevent overfitting, structural constraints are imposed on operator nesting (e.g. disallowing nested exponentials). Each of these streamwise parameters (
$a$
,
$\mu$
,
$\sigma$
) is regressed independently using the same input variable
$\dot {x}$
, with 20 000 iterations to yield concise symbolic expressions. This process generates individual Pareto fronts, which are then combined to form a new, comprehensive Pareto front
${\mathbb{P}}$
of volumetric forcing. The final loss
$L$
and total complexity
$C$
are calculated by combining the individual components:
$L=\lVert {\mathcal{G}}({\mathcal{H}}(x),y)-{\mathbb{Y}} \rVert ^2_2 /N$
and
$C=c_{a}+c_{\mu }+c_{\sigma }$
, respectively, where
${\mathcal{G}}$
represents (2.14) or (2.15),
${\mathcal{H}}$
denotes SR-discovered parameters, and
${\mathbb{Y}}$
are the true volumetric force values.
Unlike synthetic benchmarks with known target expressions, turbine wake modelling as a real-world problem lacks a definitive ground truth, requiring an additional criterion to select the most physically meaningful expression from the Pareto front
${\mathbb{P}}$
. The objective is to find the optimal balance between a high regression quality and the interpretability of a simple expression. To achieve this, a scoring system is used to quantitatively assess whether the benefit of increasing complexity outweighs the cost. The score for each expression is defined as
\begin{equation} S_i=-\frac {\log (L_i/L_{i-1})}{C_i-C_{i-1}}, \end{equation}
where
$i$
is the expression’s position in
${\mathbb{P}}$
. A high score indicates a substantial reduction in loss for only a small increase in complexity. The final expression is selected by finding the highest score among all expressions whose loss is within
$1.2$
times the minimum loss observed in
${\mathbb{P}}$
. To maintain interpretability and prevent overly complex expressions, a maximum complexity of
$20$
is set for each of the parameters.
Figure 7. Streamwise evolution of parameters for (a,b) pressure and (c,d) turbine force terms. Blue, orange and green lines represent the PS case ground truth (GT) results, SR-discovered expressions, and mean values of ST cases, respectively.
All expressions are derived exclusively from the PS case, which serves as the training dataset. The remaining ST cases are used for verification to assess the generalisability of the discovered optimal expressions. These expressions explicitly encode the streamwise evolution of the volumetric forcing, and are directly embedded into the DWM numerical solver, forming the physical core of the SRDWM framework.
3. Results and discussion
This section presents a comprehensive evaluation of SRDWM performance against LES, treated as ground truth, under various inflow conditions. For reference, excluding inflow generation, LES require on average approximately
$1.6\times 10^{-1}$
core-hours per physical second of simulated time per turbine wake case on a computing cluster, whereas SRDWM achieves approximately
$3.8\times 10^{-5}$
core-hours per physical second on a personal computer virtual machine. This corresponds to an efficiency gain over three orders of magnitude. The resulting low runtime opens the possibility of edge-computing deployment and real-time wind farm control, where rapid yet physically consistent predictions are required. This section first presents the extracted symbolic expressions (§ 3.1), then analyses key features of wake evolution (§ 3.2), and finally further examines turbulence statistics and wake dynamics (§ 3.3). Results demonstrate that SRDWM delivers high accuracy and generalisability with full spatiotemporal resolution.
3.1. The SR-discovered volumetric terms and boundary conditions
Figure 7 displays the optimal symbolic expressions discovered from the PS case for parameters
$a$
,
$\mu$
and
$\sigma$
of volumetric terms in (2.14) and (2.15). Remarkably, the amplitude function of the pressure term includes an exponential decay factor, capturing the localised pressure drop region behind the rotor. This indicates that pressure recovery occurs downstream of the actuator region rather than immediately behind the turbine, in agreement with physical expectations. The SR expression for the pressure standard deviation includes a cosine term modulated by a polynomial. This concise form captures the observed spanwise non-monotonic variation – initial broadening, subsequent narrowing, and a mild re-broadening – with a downstream weakening of amplitude, consistently across all cases. Although approximate, it improves interpretability and mitigates overfitting relative to more complex alternatives. Because the expression is selected by minimising mean error over all downstream stations, it should be viewed as a compact surrogate for nonlinear behaviour rather than evidence of true periodicity, with validity restricted to the training domain (
$x \lt 10D$
). For the turbine force distribution, the extracted DG profile shows a near-linear streamwise decay of amplitude, with constant mean and standard deviation. This behaviour reflects how momentum loss from turbine extraction rapidly weakens downstream while maintaining a consistent spanwise pattern. These findings validate the accuracy and physical interpretability of the SR-derived expressions.
Figure 8. Distribution of (a) mean pressure and (b) mean turbine force: (i) LES colour maps and isocontours for the PS case; (ii) SR predictions for the PS case (colour map) with isocontours from all LES ST cases.
Figure 8 further assesses generalisation performance by comparing the reconstructed pressure and turbine forcing in both the PS case (training data) and ST cases (test data). In figure 8(i), background colour maps show LES data for the PS case, while figure 8(ii) overlays the SR predictions with LES isocontours for all ST cases. The close alignment across all cases demonstrates that near-wake pressure and turbine forcing are governed primarily by turbine aerodynamics and remain largely invariant with respect to inflow turbulence. Far-wake variations, while present, are minor relative to the dominant near-wake features. This result highlights the robustness and strong generalisation capacity of the SR-discovered equations for volumetric forcing.
For the wake velocity boundary conditions used in SRDWM, we similarly decompose the velocity profile into fixed spanwise functions, with streamwise position set to the rotor location (
$\dot {x} = 0$
). The streamwise and radial velocity components are discovered respectively as
\begin{equation} \frac {\Delta U_x(\dot {x}=0)}{U_{\infty }} = 0.27 \times \left (\exp \left (-\frac {(\dot {r} -0.27)^2}{2\times 0.19^2}\right )+\exp \left (-\frac {(\dot {r}+0.27)^2}{2\times 0.19^2}\right ) \right )\!, \end{equation}
\begin{equation} \frac {\Delta U_r(\dot {x}=0)}{U_{\infty }} = \frac {0.11\dot {r}}{e^{2.0\dot {r}} - 4.4\dot {r}}. \end{equation}
As shown in figure 9, the SR-derived expressions accurately reconstruct the initial mean wake velocity distribution across varying inflow turbulence intensities. The spanwise location for the peak of streamwise velocity deficit aligns with that of the turbine forcing, and both variables remain largely insensitive to inflow variations. This consistency reflects the underlying blade aerodynamics (Keane Reference Keane2021), particularly the lift distribution peaking near mid-span (Magnusson Reference Magnusson1999). These forces are projected onto the mesh via Gaussian smoothing (Churchfield et al. Reference Churchfield, Lee, Michalakes and Moriarty2012), and iteratively generate the observed velocity profile through the momentum equations. Importantly, these discovered symbolic expressions not only maintain physical consistency, but also provide new interpretable insights into turbine wakes.
Figure 9. Comparison of mean wake velocity boundary conditions derived from SR and LES under different inflow turbulence intensities.
3.2. Mean and instantaneous wake evolution
Figure 10 presents the streamwise time-averaged velocity deficit and turbulence intensity increments in turbine wakes across inflow cases, computed as the differences between fields with and without the turbine. The statistics are obtained over the last
$1600$
s of each wake simulation after reaching statistically steady state. The SRDWM framework accurately reconstructs wake evolution with full spatiotemporal resolution, including both time-averaged velocity and turbulence intensity – quantities critical to wind energy research. It shows strong robustness and generalisation, automatically adapting to different inflow conditions. In contrast, the baseline DWM substantially overestimates the near-wake velocity deficit, and markedly underestimates turbulence intensity increments, particularly under low ambient turbulence. The quantitative comparison of errors across all cases is provided in table 1, highlighting the significant improvements of SRDWM over traditional DWM in both time-averaged and fluctuating predictions within the wake widths. A notable advantage of SRDWM lies in its accurate representation of the near-wake region, which has remained a persistent challenge for conventional models. This improvement is further illustrated in figure 11, where SRDWM reproduces instantaneous wake velocity more faithfully to LES results than the baseline method.
Table 1. Normalised MSE of streamwise time-averaged velocity deficit (
$\Delta U/U_{\infty }$
) and turbulence intensity increments (
$\Delta \sigma _u/U_{\infty }$
), evaluated within the wake widths, with DWM and SRDWM compared against LES ground truth.
Figure 10. Comparison of (a,b,c) streamwise time-averaged velocity deficit and (d,e,f) turbulence intensity increments from (a,d) LES, (b,e) SRDWM, and (c,f) baseline DWM, under varying inflow turbulence conditions: (i) PS case; (ii)–(vi) ST cases with increasing turbulence intensity.
Figure 11. Instantaneous wake velocity fields from (a) LES, (b) SRDWM, and (c) conventional DWM.
Interestingly, figure 10 shows a consistent asymmetry in the time-averaged velocity deficit, with larger deficits on the
$\dot {y} \lt 0$
side. The asymmetry in time-averaged velocity and turbulence intensity can be attributed to rotor rotation-enhanced mixing between higher-momentum flow above the hub and lower-momentum flow below (Placidi, Hancock & Hayden Reference Placidi, Hancock and Hayden2023; Bastankhah et al. Reference Bastankhah, Zunder, Hydon, Deebank and Placidi2024). Stronger ambient turbulence further enhances this mixing, resulting in more pronounced asymmetries for turbulence intensity in ST cases (Vinod, Han & Banerjee Reference Vinod, Han and Banerjee2021). Such spanwise asymmetries under high turbulence intensity have also been documented recently (Wang et al. Reference Wang, Guan, Wang and Lu2024b
). However, this asymmetry in the PS case differs from the ST cases, which can be linked to spanwise non-uniformity in the inflow itself (Chamorro & Porté-Agel Reference Chamorro and Porté-Agel2010). Specifically, the lower velocity on the
$\dot {y} \lt 0$
side is observed in figure 12(i), which alters wake development. Despite these complexities, SRDWM effectively captures wake dynamics, highlighting its adaptability and robustness under diverse inflow conditions. Moreover, increasing turbulence intensity accelerates wake recovery, a well-established observation. Notably, ST3 shows a more pronounced turbulence intensity increase than PS, despite similar incoming turbulence intensity, indicating additional influencing factors to be explored in § 3.3.
Figure 12. (a) Time-averaged streamwise velocity contours, and (b) spanwise profiles of its statistics under varying inflows: (i) PS case; (ii)–(vi) ST cases with increasing turbulence intensity. The blue, green and orange lines denote the maximum, mean and minimum values across the streamwise extent, while the shaded grey band indicates one standard deviation around the mean.
Figure 13 presents global quantitative verification using evenly spaced downstream probes from the PS case and from the cases with the largest MSE in table 1. For the PS inflow, results from the more accurate but computationally expensive ALM are also provided as a benchmark. In line with previous studies (Porté-Agel et al. Reference Porté-Agel, Wu, Lu and Conzemius2011; Wu & Porté-Agel Reference Wu and Porté-Agel2011; Qian et al. Reference Qian, Song and Ishihara2022), the velocity deficit and turbulence intensity predicted by ADM and ALM are generally consistent, with ADM slightly underestimating peak turbulent fluctuations in the near wake. The SRDWM predictions show close agreement with LES data for both time-averaged velocity and turbulence intensity across the full wake, including the near-rotor region. These results suggest that SRDWM can deliver engineering-focused wake simulations with accuracy comparable to LES, while requiring substantially lower computational cost.
Figure 13. Profiles of (i) time-averaged streamwise velocity deficit and (ii) turbulence intensity increment at various downstream locations. Results are shown for (a) the PS case, and (b) the cases with the largest MSE in table 1, namely (bi) ST1 and (bii) ST4. The LES with ADM (ground truth) are shown as circles, SRDWM as solid lines, and PS-ALM as dotted lines. Grey vertical dashed lines indicate selected downstream positions.
Our proposed SRDWM framework integrates two main improvements to the baseline DWM model: the SR closures and the revised WAT formulation. To systematically disentangle the respective contributions of these two components, we conducted a sensitivity analysis, presented in figure 14. This figure contrasts LES results with four model configurations: the SR closures alone, where SR-discovered equations replace DWM’s original empirical terms; the WAT correction alone, where the revised WAT formulation is added to the baseline DWM; their combined form (SRDWM); and the baseline DWM. Within the SRDWM workflow, the SR closures first refine the quasi-steady wake velocity through a thin shear layer approximation of the NS equations. Subsequently, the revised WAT formulation utilises this improved wake field to selectively scale velocity fluctuations when computing the instantaneous wake velocity in the meandering frame. Accordingly, the SR closures are primarily responsible for improvements in time-averaged velocity, while the revised WAT formula specifically corrects turbulent fluctuations, and their synergy underpins the overall performance gain achieved by SRDWM.
Figure 14. Contours of (a) time-averaged streamwise velocity deficit, and (b) turbulence intensity increments for the PS case, comparing results from (i) LES, (ii) SRDWM, (iii) DWM, (iv) SR closures only, and (v) revised WAT only.
To further quantify wake recovery, figure 15(a) presents the time-averaged velocity deficit along the wake centreline. Enhanced inflow turbulence clearly accelerates wake recovery, whereas ST1 – with very low inflow turbulence – exhibits little recovery even at
$10D$
downstream. The SRDWM accurately captures these trends across all conditions, with minor deviations in ST1. As shown in figure 15(b), turbulence intensity is suppressed in the near wake and enhanced downstream. Maximum rotor-extracted energy consistently occurs at
$x = 0.7D$
along the centreline, while peak wake-added energy appears near the rotor edge (
$r = 0.5D$
) in the far wake. As inflow turbulence increases, this location shifts upstream, stabilising at approximately
$x = 5.7D$
. These two characteristic locations, denoted as
$P_1$
and
$P_2$
, are therefore selected as representative monitoring points for the near- and far-wake regions in subsequent full-domain flow analysis.
Figure 15. (a) Time-averaged streamwise velocity deficit along the wake centreline from LES (circles) and SRDWM (lines). (b) Locations of peak rotor-extracted energy (LES denoted
$\times$
, SRDWM denoted
$\square$
) and wake-added energy (LES denoted
$+$
, SRDWM denoted
$\vartriangle$
). Dashed lines at
$x = 0.7D$
,
$x = 5.7D$
,
$y = 0$
and
$y = 0.5D$
indicate monitoring references. Red points
$P_1$
and
$P_2$
indicate the selected near-wake and far-wake monitoring locations for further analysis.
Another major limitation of conventional DWM models lies in their incorrect representation of spanwise velocity. Figure 16 compares spanwise velocity from LES, SRDWM and the traditional DWM model. In the latter, the spanwise velocity at the boundary is fixed at zero, while a large streamwise velocity deficit is imposed. As a result, solving the continuity equation yields radial velocity
$U_r\lt 0$
in the near wake, artificially compensating for this imbalance. This leads to a misrepresentation of the wake flow’s outward expansion as inward contraction (Rivera-Arreba et al. Reference Rivera-Arreba, Li, Yang and Bachynski-Polić2024), violating physical laws. In contrast, SRDWM accurately captures spanwise wake dynamics in close agreement with LES, via physically consistent boundary and forcing inputs. Identity plots in figure 16(d,e) confirm that the conventional model severely underperforms, even scoring below zero in R-squared metrics, indicating systemic directional errors. These results demonstrate the proposed framework’s success in overcoming inherent deficiencies of conventional approaches.
Figure 16. Time-averaged spanwise velocity deficit in the PS case from (a) LES, (b) SRDWM, and (c) conventional DWM. Identity plots comparing LES (
$Y$
,
$x$
-axis) and predicted (
$\widetilde {Y}$
,
$y$
-axis) values of normalised radial velocity from (d) SRDWM and (e) conventional DWM across all cases. Colours represent radial position. Dotted diagonal shows perfect match; dashed lines are zero baselines.
Figure 17. (a) Wake centre deflection along the streamwise direction, (b) wake width, and premultiplied PSD of wake centre deflection at (c)
$x=0.7D$
and (d)
$x=9.5D$
under varying turbulence intensities. The LES are shown as solid lines and circles, SRDWM as dashed lines, and baseline DWM as dash-dotted lines.
Wake centre deflection and wake width across cases are analysed in figure 17, where deflection statistics are derived from the instantaneous momentum-deficit centre (Vahidi & Porté-Agel Reference Vahidi and Porté-Agel2024). A minimum wake-centre displacement occurs near
$x=0.7D$
, indicating that both turbulent velocity fluctuations and wake deflection reach their lowest levels at this location. This observation supports the DWM hypothesis that wake meandering responds linearly to velocity fluctuations. The SRDWM successfully reproduces these near-wake dynamics, as further confirmed by the premultiplied PSD at this position (figure 17
c). However, the traditional DWM substantially underestimates low-frequency wake deflection from the near- to far-wake regions, primarily due to its underprediction of turbulent fluctuations. In contrast, the proposed WAT module within the SRDWM framework selectively scales wake velocity fluctuations, thereby reducing this bias even at the downstream location where the SRDWM error is largest (figure 17
d). Although discrepancies in magnitude remain, both the passive tracer-based DWM and SRDWM accurately capture the dominant frequencies of wake deflection, which reinforces the validity of the DWM hypothesis. Increased incoming turbulence intensity results in enhanced wake deflection and expansion, a trend well recovered by our model. However, wake deflection and widths are slightly overestimated, consistent with previous findings (Feng et al. Reference Feng, Gupta, Li and Wan2024a
), which attributed such discrepancies to intrinsic limitations of DWM-type models. Because the correlation between high-frequency fluctuations and wake meandering is generally considered negligible (Muller, Aubrun & Masson Reference Muller, Aubrun and Masson2015), small-scale wake deflection is filtered out in these models. Consequently, SRDWM performs less accurately for
$St\gt 0.3$
, tending to overestimate wake meandering and, subsequently, the wake width.
3.3. Turbulence statistics and wake dynamics
Figure 18 presents the premultiplied PSD of turbulence kinetic energy (TKE), defined as
$f\phi _k = f (\phi _{u^{\prime }} + \phi _{v^{\prime }} + \phi _{w^{\prime }})$
, where
$u^{\prime }, v^{\prime }, w^{\prime }$
are velocity fluctuation components. The spectra are evaluated along the rotor edge at two downstream locations. Immediately behind the rotor (
$x=0$
), turbine operation enhances TKE primarily in the low-frequency range (
$St \lt 1$
), indicating the presence of large-scale coherent structures. As the flow progresses downstream (
$x = 5.7D$
), the influence of the incoming turbulence intensity on the low-frequency TKE gradually diminishes. In particular, lower inflow turbulence results in higher wake-added TKE contributions from the turbine, consistent with observations in figure 10. Across all cases, SRDWM reproduces turbulence energy observed in LES, confirming its ability to resolve spectral wake characteristics.
Figure 18. Premultiplied PSD of TKE along the rotor edge at (a)
$x=0$
and (b)
$x=5.7D$
. Solid and dashed lines represent LES and SRDWM results, respectively. Dotted lines indicate cases without the turbine.
Notably, despite similar inflow turbulence intensities, the PS case exhibits lower wake TKE compared to the ST cases. Following the interpretation of Li et al. (Reference Li, Zhang, Li and Yang2024), this phenomenon is attributed to the smaller turbulence integral length scale present in the PS condition. Their findings suggest that wake-added TKE is more pronounced when the inflow has either larger integral length scales or lower turbulence intensities. To further investigate the spatial structure of wake turbulence, we conduct an integral length scale analysis using two-point velocity correlation functions. The streamwise velocity autocorrelation is defined as (Vahidi & Porté-Agel Reference Vahidi and Porté-Agel2024)
\begin{equation} R_{uu}(\delta \chi ) = \frac {\langle u^{\prime }(\chi ) \, u^{\prime }(\chi + \delta \chi ) \rangle }{\sigma _{u^{\prime }(\chi )} \, \sigma _{u^{\prime }(\chi + \delta \chi )}}, \end{equation}
where
$\delta \chi$
denotes the spatial separation, and
$\chi$
represents either the streamwise or spanwise direction. The corresponding integral length scale is calculated as
\begin{equation} L^{\chi }_{uu} = \int _0^{\infty } R_{uu}(\delta \chi ) \, {\rm d}\delta \chi . \end{equation}
As shown in figure 19, turbine operation reduces spatial coherence relative to no-turbine cases, reflecting large eddies breakdown into small-scale structures. The SRDWM closely tracks LES across directions and inflow conditions, confirming its ability to faithfully reproduce the turbine-induced structural modifications in wakes.
Figure 19. (a) Spanwise and (b) streamwise integral length scales of the streamwise velocity component. Circles and dashed lines represent LES and SRDWM results, respectively. Dotted lines indicate cases without the turbine.
Figure 20. Spanwise profiles at
$x=0.7D$
of (a) skewness, (b) kurtosis, and (c) PDF of streamwise velocity fluctuations at point
$P_1$
. Circles and dashed lines represent LES and SRDWM, respectively.
High-order turbulence statistics are shown in figure 20. At monitoring point
$P_1$
, the streamwise velocity fluctuations display a near-Gaussian PDF, with skewness and kurtosis close to Gaussian reference values
$S_u = 0$
and
$K_u = 3$
. These metrics are defined as
\begin{equation} \begin{aligned} S_u &= \frac {\langle u^{\prime 3} \rangle }{\langle u^{\prime 2} \rangle ^{3/2}}, \\ K_u &= \frac {\langle u^{\prime 4} \rangle }{\langle u^{\prime 2} \rangle ^{2}}. \end{aligned} \end{equation}
This observation suggests that turbulence at this TKE-suppressed location is relatively homogeneous, with energy primarily distributed among continuous small-scale fluctuations rather than intermittent large-scale bursts. The SRDWM accurately replicates these statistical properties, demonstrating its fidelity in capturing not only first- and second-order but also higher-order turbulence characteristics.
Figure 21. (a) The Q2 and Q4 contributions at
$x=0.7D$
; (b) self-similarity of Q2/Q4 profiles for
$x \geq 3D$
in the ST1 case; (c) streamwise evolution of Q2 and Q4 at the rotor edge; (d) normalised difference in contribution to Reynolds stress between Q2 and Q4,
$\Delta C_{\eta }=(C_{2,\eta }-C_{4,\eta })/(C_{2,\eta }+C_{4,\eta })$
. Markers, dashed lines and dotted lines represent LES, SRDWM and DWM, respectively; vertical dotted lines mark
$x=0.7D$
and
$x=5.7D$
; horizontal line is zero baseline.
To investigate momentum transport in turbine wakes and characterise the dynamics at monitoring points
$P_1$
and
$P_2$
, we employ quadrant analysis to quantify the contributions of coherent turbulence structures to the Reynolds shear stress (Wallace Reference Wallace2016). The analysis is conducted in a coordinate system where the horizontal and vertical axes represent
$u^{\prime }$
and
$v^{\prime }$
, respectively, dividing the space into four quadrants
${\rm Q}i$
that correspond to different turbulence interactions. The relative contribution of each quadrant is defined as follows:
\begin{equation} \begin{aligned} -\langle u^{\prime }v^{\prime } \rangle _i &= \langle -u^{\prime }v^{\prime } I_i(u^{\prime }, v^{\prime })\rangle , \\[3pt] I_i(u^{\prime }, v^{\prime }) &= \begin{cases} 1 & \text{if } (u^{\prime }, v^{\prime }) \text{ lies in quadrant } {\rm Q}i, \\ 0 & \text{otherwise}. \end{cases} \end{aligned} \end{equation}
This conditional sampling method reveals the flow directionality relative to the mean flow. In figure 21, we focus on the contribution of ejection events (Q2) and sweep events (Q4), which are the primary mechanisms driving wake flow transport (Li et al. Reference Li, Zhang, Li and Yang2024; Xiong et al. Reference Xiong, Laima, Li and Zhou2024). The evolution of these events strongly influences turbine power output and unsteady loads (Kadum, Knowles & Cal Reference Kadum, Knowles and Cal2019). Moreover, characterising the size, organisation and relative contribution of these turbulent structures provides insights for improved wind farm design and for understanding wind farm–atmosphere interactions (Hamilton et al. Reference Hamilton, Kang, Meneveau and Cal R>B>2012).
In the near-wake region, Q2 and Q4 contribute comparably in the spanwise direction, reaching their maximum at the shear layer. At the rotor edge, sweeps dominate slightly inside the rotor disk, whereas ejections dominate slightly outside (Ali et al. Reference Ali, Hamilton, Calaf and Cal2019; Mouchref et al. Reference Mouchref, Viggiano, Ferčák, Bossuyt, Ali, Meneveau, Gayme and Cal2024). The baseline DWM exhibits small overall errors but overestimates momentum transport near the wake centre within the rotor energy-extraction region. By contrast, SRDWM accurately replicates these patterns, including the spanwise self-similarity of Q2 and Q4 events in low-turbulence conditions, consistent with prior findings (Xiong et al. Reference Xiong, Laima, Li and Zhou2024) and illustrated in figure 21(b). This improvement mainly arises from the revised WAT formulation, which locally suppresses turbulence within the rotor energy-extraction region while enhancing velocity fluctuations in the outer and downstream regions, yielding a more accurate representation of momentum transport. The self-similarity is observed after normalising the peak event values and their corresponding spanwise locations at downstream positions (
$x \geq 3D$
). Farther downstream, wake expansion drives outdraft-related ejections, which become the dominant mechanism of momentum transport. While SRDWM effectively characterises sweep events, it slightly underestimates the contribution of outdraft-driven ejections in the far wake, particularly under high-turbulence conditions. In contrast, the baseline DWM always markedly underpredicts wake turbulent fluctuations, resulting in a substantial underestimation of quadrant events – a deficiency that becomes even more pronounced when normalised to evaluate self-similarity under low-turbulence inflow.
To further assess the influence of intense turbulent events, we conduct quadrant-hole analysis (Zhou, Liu & Wan Reference Zhou, Liu and Wan2023), introducing a thresholded indicator function
$I_{i,\eta }$
:
\begin{equation} I_{i,\eta }(u^{\prime }, v^{\prime }) = \begin{cases} 1 & \text{if } (u^{\prime }, v^{\prime }) \text{ lies in quadrant } {\rm Q}i \text{ and } |u^{\prime }v^{\prime }| \geq \eta \sigma _u \sigma _v, \\ 0 & \text{otherwise}. \end{cases} \end{equation}
Here,
$\sigma _u$
and
$\sigma _v$
denote the standard deviation of streamwise and spanwise velocity, respectively. This filtering enables isolation of high-intensity events above a prescribed threshold
$\eta$
. The relative contribution of such events to the total momentum flux is then expressed as
\begin{equation} C_{i,\eta } = \frac {\langle u^{\prime }v^{\prime } I_{i,\eta }(u^{\prime }, v^{\prime }) \rangle }{\langle u^{\prime }v^{\prime } \rangle }. \end{equation}
As shown in figure 21(d), sweep events dominate the rotor-extracted energy region, where limited outdrafts restrict momentum mixing (Viestenz & Cal Reference Viestenz and Cal2016). In contrast, wake-added energy in the far wake is primarily associated with outdrafts induced by ambient turbulence, which enhance momentum exchange between high- and low-velocity regions. This effect becomes increasingly pronounced with rising hole thresholds
$\eta$
, emphasising the key role of extreme coherent structures in governing wake dynamics.
Turbulent coherent structures dominate wind velocity fluctuations, influence aerodynamic performance, and shape the wake behaviour of wind turbines within the ABL (Wang et al. Reference Wang, Guan, Wang and Lu2024b ). To analyse these structures, we employ proper orthogonal decomposition (POD) (Berkooz, Holmes & Lumley Reference Berkooz, Holmes and Lumley1993), which extracts flow structures that maximise variance, making it an effective tool for identifying meandering-induced features in wakes (Hamilton et al. Reference Hamilton, Doubrawa, Moriarty, Letizia and Thedin2025). Beyond its use in reduced-order modelling, POD also provides physical insight into dominant coherent structures in turbulence (Moreno et al. Reference Céspedes Moreno, Leon and Andersen2025), while the associated eigenvalues quantify the relative TKE of each mode (Hamilton et al. Reference Hamilton, Viggiano, Calaf, Tutkun and Cal2018). In this study, POD is implemented by performing singular value decomposition on the full-field velocity fluctuation matrix, with the domain spanning both near- and far-wake regions to ensure comprehensive coverage. The decomposition yields a set of orthogonal modes (spatial structures) and their singular values, whose squares represent the energy content of each mode. All modes are retained in the analysis, with the total number equal to the number of instantaneous snapshots. Bastine et al. (Reference Bastine, Witha, Wächter and Peinke2015) showed that in their case, the first two POD modes are sufficient to represent a dynamic wake model when wake meandering is considered as the only dynamic feature. Therefore, figure 22 presents the first two POD modes of the streamwise and radial velocity components. The SRDWM successfully reproduces the dominant large-scale wake structures observed in LES, where these oscillations are typically linked to wake meandering (Sørensen et al. Reference Sørensen, Mikkelsen, Henningson, Ivanell, Sarmast and Andersen2015; De Cillis et al. Reference De Cillis, Cherubini, Semeraro, Leonardi and De Palma2020). In contrast, the baseline DWM exhibits clear discrepancies, whose POD modes appear overly smoothed and misaligned relative to LES, indicating its limited ability to capture the correct spatial distribution of wake meandering.
Figure 22. Comparison of the first two POD modes of (i,ii) streamwise and (iii,iv) radial velocity in the ST1 case: (a) LES, (b) SRDWM, and (c) baseline DWM.
Figure 23. Variation with POD rank of (a) individual energy contribution, (b) cumulative energy, and (c) mean frequency
$St_m$
of POD modes. The LES are shown as circles, SRDWM as dashed lines, and baseline DWM as dotted lines.
The relative energy contributions of the first
$200$
modes are expressed as percentages, with both individual and cumulative values normalised by the total energy of these modes, as shown in figure 23. The resulting distribution shows little sensitivity to inflow turbulence intensity, consistent with field measurements reporting similar modal ordering under different atmospheric conditions (Hamilton et al. Reference Hamilton, Doubrawa, Moriarty, Letizia and Thedin2025). To further characterise frequency content, a mean Strouhal number
$St_m$
is defined as the geometric mean of the mode frequency obtained from the PSD of temporal coefficients. As expected, lower-order modes correspond to lower frequencies and represent the large-scale, energetic structures (Cherubini et al. Reference Cherubini, De Cillis, Semeraro, Leonardi and De Palma2022; Wang et al. Reference Wang, Guan, Wang and Lu2024b
). Differences across cases are most evident at low frequencies, where large-scale atmospheric structures strongly influence wake dynamics (Ali et al. Reference Ali, Cortina, Hamilton, Calaf and Cal2017; Zhang & Stevens Reference Zhang and Stevens2020). Most modes fall within
$0.1 \lt St_m \lt 1$
, aligning with the model’s filtering bandwidth. Compared with traditional DWM, which struggles to capture the most energetic POD modes, SRDWM more effectively resolves the multiscale dynamics of turbine wakes owing to the revised WAT model. This leads to a closer agreement with LES, particularly in the dominant-mode
$St_m$
, which the baseline DWM markedly underestimates under low-turbulence inflow.
4. Conclusion
We have developed a symbolic regression-enhanced dynamic wake meandering (SRDWM) model that significantly advances the capabilities of engineering wake simulations. Unlike prior DWM variants that rely on empirical tuning, this framework embeds symbolic expressions discovered from high-fidelity LES data, directly into the DWM model to reconstruct both volumetric forcing and wake boundary conditions. By incorporating these interpretable expressions, SRDWM achieves equation-level closure – enhancing physical consistency, particularly in the near wake, which has long posed a challenge for conventional approaches.
A key innovation of this work lies in the use of a hierarchical symbolic regression (SR) strategy informed by domain knowledge. By decomposing the expression space and imposing physics-informed structural constraints, we reduce the dimensionality of the learning problem and improve the interpretability and generalisation of the extracted expressions. This approach supports robust equation discovery from limited data and reinforces the physical fidelity of the resulting model.
To further improve wake turbulence predictions, we introduce a revised wake-added turbulence formulation that selectively amplifies energy-dominant fluctuations while suppressing unphysical growth near the rotor. This enables accurate reproduction of both wake spectra and turbulence statistics across a range of inflow conditions.
The SRDWM framework retains the temporal resolution capabilities of DWM while incorporating enhanced spatial integrity. It accurately captures the evolution of time-averaged profiles and turbulent properties from the rotor to the far wake. Compared to LES, SRDWM offers comparable accuracy with runtime speed-ups of over three orders of magnitude. This makes it well suited for wind farm layout optimisation, control co-design, and operational forecasting. Additionally, its high computational efficiency enables the generation of large-scale wake data for training and benchmarking modern machine learning models.
These findings underscore the potential of SR as a physics-informed modelling paradigm that bridges data-driven discovery and mechanistic understanding. Unlike black-box data-driven models, SR emphasises transparency and physical interpretability, providing deeper understanding into turbine wake dynamics. Moreover, the proposed SR framework can serve as a general tool for constructing closure models that explicitly correct empirical formulas derived under simplified assumptions, a common need across many physical domains.
Despite its demonstrated advantages under the tested cases, SRDWM is developed from LES data at a single hub-height wind speed under neutral stability, and its fidelity in markedly different regimes – such as strongly stratified atmospheres, large wind veer, or substantially different inflow speeds – remains to be quantified. These, however, represent natural directions for extension within the same physics-informed SR framework. Stratification can be incorporated by adding buoyancy forcing term in the symbolic search, while veer and inflow-speed effects can be addressed through multi-condition training, inclusion of additional SR candidate operands, or informed priors and transfer-learning strategies. The framework’s fidelity also depends on the reference LES data; although SOWFA and the NREL 5 MW model are well validated, newer solvers and turbine models could provide richer datasets for modern turbines. Future work will therefore emphasise multi-condition training, targeted sensitivity tests (including stratification and veer), transfer-learning approaches, and advanced datasets to enhance robustness and broaden applicability.
Acknowledgements
The authors would like to thank Professor V. Gupta for his valuable feedback, which helped to improve this paper. The authors also acknowledge support from the High Performance Computing Centres at Eastern Institute of Technology, Ningbo, and Ningbo Institute of Digital Twin.
Funding
This work was supported by the National Natural Science Foundation of China (no. 12572266), the National Key Research and Development Programme of China (2024YFF1500600), the Innovation Capability Support Programme of Shaanxi (programme no. 2023-CX-TD-30), and the Yongjiang Talent Programme of Ningbo (2022A-242-G).
Declaration of interests
The authors report no conflict of interest.
Code availability
The symbolic regression tool used in this study is publicly available at https://github.com/MilesCranmer/PySR and https://github.com/Scientific-Artificial-Intelligence-Lab/kd.
Open access
