JHTDB-wind: a web-accessible large-eddy simulation database of a wind farm with virtual sensor querying

This paper introduces JHTDB-wind (https://turbulence.idies.jhu.edu/datasets/windfarms, last access: 11 November 2025), a publicly accessible database containing large-eddy simulation (LES) data from wind farms. Building on the framework of the Johns Hopkins Turbulence Database (JHTDB), which hosts direct numerical simulation (DNS) and some LES datasets of canonical turbulent flows, JHTDB-wind stores the 4D space–time history of the flow and provides users the ability to access and query the data via a web-based virtual sensor interface. The initial dataset comprises LES results from a large wind farm with 10×6 turbines, modeled using a filtered actuator line method, under conventionally neutral atmospheric conditions. These data comprise 1 h (hour) of flow field data (velocity, pressure, potential temperature deviation, subgrid-scale (SGS) eddy viscosity, and turbine forces, approximately 15 TB (terabytes) and wind turbine data – including both turbine-level operational quantities and blade-level aerodynamic quantities (approximately 1.3 TB) – stored in Zarr and Parquet formats, respectively. Data retrieval is facilitated by the giverny Python package, allowing remote users to query the database in Python or MATLAB (C and Fortran support are available for flow field data). This paper details the simulation setup and demonstrates data access through examples that analyze wind farm flow structures and turbine performance. The framework is extensible to future datasets, including the JHTDB-wind diurnal cycle simulation analyzed in Xiao et al. (2025).

1 Introduction

Eddy-resolving simulations of atmospheric boundary layer (ABL) phenomena (Porté-Agel et al., 2000; Bou-Zeid et al., 2004; Kumar et al., 2006) and of wind farms in particular (Calaf et al., 2010; Meyers and Meneveau, 2012; Gebraad et al., 2016; Stevens and Meneveau, 2017; Zhang et al., 2023) have significantly advanced our understanding of the complex, multi-scale, and multi-physics processes involved. Large-eddy simulations (LESs) offer high spatial and temporal resolution, capturing the dynamics of relatively small and fast turbulent eddies (Churchfield et al., 2012; Chatelain et al., 2013; Yang et al., 2021; Li et al., 2022). While the range of resolved scales in LES is constrained by computational resources, the number of LES grid points in typical simulations continues to increase. However, data handling and post-processing capabilities have not kept pace with the resulting rapid increase in data volumes. For instance, a single LES of turbulent flow outputting five field variables (e.g., the three velocity components, potential temperature, and pressure) on 2048³ spatial grid points and integrated over, say, 10⁴ time steps (McWilliams et al., 1994; Alexakis et al., 2024) can generate petabytes (PB) of data. As a result, most studies store only a few selected snapshots and rely heavily on pre-defined run-time diagnostics when time-resolved analysis is required. This approach reduces storage requirements but limits the ability to revisit data when new questions and concepts arise, often necessitating costly recomputation. Furthermore, certain analyses – such as backward-in-time particle tracking from an extreme dissipation event – cannot be performed without the full temporal data.

To address these challenges, modern database technologies have increasingly been applied to preserve and store data from simulation-based turbulence research (Perlman et al., 2007; Zhang et al., 2018; Chung et al., 2022; Duraisamy et al., 2019). One example is the Johns Hopkins Turbulence Database (JHTDB; https://turbulence.idies.jhu.edu, last access: 11 November 2025), an open-access platform supported by the National Science Foundation (Perlman et al., 2007; Li et al., 2008). JHTDB enables researchers to interact with easily accessible large-scale simulation data. The system currently hosts more than 1 PB of direct numerical simulation (DNS) data for canonical, turbulent flows of fundamental interest (over 2 PB if counting warm backup copies), including six space–time-resolved datasets and several others with a few snapshots available. Some LES datasets of stably stratified atmospheric turbulence are also included in JHTDB. Through web-service-based tools, users can query the database using a “virtual sensors” interface, specifying spatial and temporal locations for which the system returns properly interpolated field or derivative values (Li et al., 2008; Yu et al., 2012). A hallmark of the platform is that it allows users to access only the specific subsets of the data they require, eliminating the need to download massive datasets or manage complex file formats. This approach has significantly broadened access to high-fidelity eddy-resolving simulation data and has contributed to democratizing high-performance computational turbulence research. To date, JHTDB data have been used in research reported in over 400 peer-reviewed journal articles.

At the same time, with the growing global demand for renewable energy, enhancing wind energy efficiency has become a key priority. As wind turbines grow larger and wind farms expand in scale, their interactions with the ABL become increasingly complex – particularly with respect to wake dynamics, energy extraction, and the redistribution of momentum within the flow. LESs of large wind turbines have emerged as a crucial complement to field measurements, enabling researchers to explore flow–turbine interactions in detail and to develop engineering models that inform turbine placement strategies and improve wind farm efficiency. For example, Calaf et al. (2010) used LES with periodic boundary conditions to study the performance of “infinite” arrays of wind turbines under neutrally stratified conditions. Abkar and Porté-Agel (2013, 2014) examined how wind farm density and free-atmosphere stability influence kinetic energy fluxes in a conventionally neutral boundary layer (CNBL) – defined as neutrally stratified surface layers capped by stably stratified free atmospheres (Zilitinkevich et al., 2002). Allaerts and Meyers (2015) explored the effect of capping inversion profile on wind farm performance. Numerous additional LES-based studies have further advanced the field (Yang et al., 2014; Aitken et al., 2014; Martínez-Tossas et al., 2015; Stevens et al., 2018; Gharaati et al., 2022, 2024; Aiyer et al., 2024), highlighting the continued value of high-resolution simulation tools for understanding and optimizing wind energy systems.

These simulations, like many previous numerical studies of large-scale wind farms, generate extensive datasets. However, access to these data often remains restricted to the original researchers who conducted the simulations. The data (typically 4D space–time fields of velocity, temperature, etc.) are ephemeral: they must be analyzed in real time during the simulation, or, at best, a limited number of snapshots are stored for post-processing, while the large majority of the data is discarded. As demonstrated in the case of the JHTDB database, providing access to the 4D space–time history of a simulation could provide substantial benefits for the broader research community. The value of open access to time-resolved numerical datasets is now being recognized beyond fluid dynamics, particularly in the field of geosciences. For example, the recently released NOW-23 dataset (Bodini et al., 2024) comprises a full year of Weather Research and Forecasting (WRF) model simulations of off-shore wind conditions over several expansive (hundreds of km) US coastal regions, offering valuable data for wind farm developers. However, no equivalent open-access LES datasets currently exist at smaller scales that explicitly include wind turbine effects – datasets that would be highly valuable for researchers focused on wake interactions, turbine siting, and wind farm optimization. More in general, the lack of data sharing in the wind energy sector has been recognized to hinder technical progress and leads to missed opportunities for improving the efficiency of energy markets (Kusiak, 2016).

To begin addressing the need for open access to LES wind farm data, we construct JHTDB-wind (see https://turbulence.idies.jhu.edu/datasets/windfarms, last access: 11 November 2025; Zhu et al., 2025), a publicly accessible turbulence database built on the JHTDB framework. This paper presents the dataset by detailing the simulation framework (Sect. 2) and flow configuration – specifically, a CNBL interacting with a 60-turbine wind farm using National Renewable Energy Laboratory (NREL) 5 MW reference turbines. Here, CNBL is chosen because it is a less complicated atmospheric state, observed in nature (Liu and Stevens, 2022), for example, during the transition period after sunset or on cloudy days with powerful winds (Allaerts and Meyers, 2017; Liu et al., 2024). Simulation parameters are described in Sect. 3. The construction of the database system is described in Sect. 4, followed by an overview of representative data access methods based on the JHTDB virtual sensor method, illustrated here via Python examples (Sect. 5). Conclusions are summarized in Sect. 6. Further documentation is available directly on the database website.

2 Large-eddy simulation framework

In this study, we use the open-source LES code LESGO (https://lesgo.me.jhu.edu, last access: 11 November 2025) as a numerical solver to simulate ABL flows and its interactions with wind turbines (Calaf et al., 2010; Stevens and Meneveau, 2017; Martinez, 2017; Stevens et al., 2018; Shapiro et al., 2018, 2020; Gharaati et al., 2022; Narasimhan et al., 2022, 2024a, b, 2025; Gharaati et al., 2024; Ayala et al., 2024). The model represents all variables on a 3D Cartesian grid, with x, y, and z denoting the streamwise, spanwise, and vertical directions, respectively. In index notation, these are expressed as x_i, where i=1, 2, 3. The corresponding velocities are denoted by u_i or also with u, v, and w for its x-, y-, and z-direction components, respectively.

2.1 Governing equations and numerical methods

The turbulent flow is simulated by solving the filtered Navier–Stokes equations in their rotational form with Boussinesq thermal forcing and Coriolis effects, along with the transport equation for the potential temperature field. The governing equations include the filtered mass conservation,

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i}{\partial x_i}}=\mathrm{0};\end{array}$$

the filtered momentum conservation,

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i}{\partial t}}& +{\stackrel{\mathrm{̃}}{u}}_j\left({\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i}{\partial x_j}}-{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_j}{\partial x_i}}\right)=-{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{p}}^{*}}{\partial x_i}}+{\displaystyle \frac{g}{{\mathit{\theta }}_{\mathrm{0}}}}(\stackrel{\mathrm{̃}}{\mathit{\theta }}-{\mathit{\theta }}_{\mathrm{0}}){\mathit{\delta }}_{i\mathrm{3}}\\ & -{\displaystyle \frac{\partial {\mathit{\tau }}_{ij}^{\mathrm{SGS},d}}{\partial x_j}}-f_i+f_{\mathrm{c}}({\stackrel{\mathrm{̃}}{u}}_{\mathrm{2}}-V_{\mathrm{g}}){\mathit{\delta }}_{i\mathrm{1}}\\ & -f_{\mathrm{c}}({\stackrel{\mathrm{̃}}{u}}_{\mathrm{1}}-U_{\mathrm{g}}){\mathit{\delta }}_{i\mathrm{2}};\end{array}\end{array}$$

and the filtered heat conservation,

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{\mathit{\theta }}}{\partial t}}+{\stackrel{\mathrm{̃}}{u}}_j{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{\mathit{\theta }}}{\partial x_j}}=-{\displaystyle \frac{\partial {\mathrm{\Pi }}_j}{\partial x_j}}.\end{array}$$

Here, the tilde indicates filtering at the LES grid scale $\stackrel{\mathrm{̃}}{\mathrm{\Delta }}=\sqrt[\mathrm{3}]{\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z}$; ρ is the density of air; ${\mathit{\tau }}_{ij}^{\text{SGS}}=\stackrel{\mathrm{̃}}{u_i{u}_j}-{\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{u}}_j$ is the subgrid-scale (SGS) stress tensor, and ${\mathit{\tau }}_{ij}^{\text{SGS,d}}={\mathit{\tau }}_{ij}^{\text{SGS}}-{\mathit{\delta }}_{ij}{\mathit{\tau }}_{kk}^{\text{SGS}}/\mathrm{3}$ is the deviatoric (trace-free) part of ${\mathit{\tau }}_{ij}^{\text{SGS}}$, where δ_ij is the Kronecker delta; $\stackrel{\mathrm{̃}}{p^{*}}=\stackrel{\mathrm{̃}}{p}/\mathit{\rho }+{\stackrel{\mathrm{̃}}{u}}_k{\stackrel{\mathrm{̃}}{u}}_k/\mathrm{2}+{\mathit{\tau }}_{kk}^{\text{SGS}}/\mathrm{3}$ is the pseudo-pressure, where $\stackrel{\mathrm{̃}}{p}$ is the resolved pressure; $g=\mathrm{9.81}\mathrm{m}{\mathrm{s}}^{-\mathrm{2}}$ is the gravitational acceleration; θ₀ is the reference potential temperature scale; and f_i is the distributed body force for modeling the turbine-induced aerodynamic forces on the air flow (see Sect. 2.3). In the present study, ${\mathit{\tau }}_{ij}^{\text{SGS,d}}$ is parameterized using the Lilly–Smagorinsky eddy-viscosity-type model (Smagorinsky, 1963; Lilly, 1966), i.e., ${\mathit{\tau }}_{ij}^{\text{SGS,d}}=-\mathrm{2}{\mathit{\nu }}_{\text{SGS}}{\stackrel{\mathrm{̃}}{S}}_{ij}=-\mathrm{2}\left(C_{\mathrm{s}}\stackrel{\mathrm{̃}}{\mathrm{\Delta }}{)}^{\mathrm{2}}\right|\stackrel{\mathrm{̃}}{S}|{\stackrel{\mathrm{̃}}{S}}_{ij}$, where ${\stackrel{\mathrm{̃}}{S}}_{ij}=\mathrm{0.5}(\partial {\stackrel{\mathrm{̃}}{u}}_i/\partial x_j+\partial {\stackrel{\mathrm{̃}}{u}}_j/\partial x_i)$ is the resolved strain-rate tensor, $\left|\stackrel{\mathrm{̃}}{S}\right|=\sqrt{\mathrm{2}{\stackrel{\mathrm{̃}}{S}}_{ij}{\stackrel{\mathrm{̃}}{S}}_{ij}}$ is the strain-rate magnitude, and ${\mathit{\nu }}_{\text{SGS}}=\left(C_{\mathrm{s}}\stackrel{\mathrm{̃}}{\mathrm{\Delta }}{)}^{\mathrm{2}}\right|\stackrel{\mathrm{̃}}{S}|$ is the modeled SGS eddy viscosity. The coefficient C_s is dynamically determined using the Lagrangian-averaged scale-dependent dynamic model (Bou-Zeid et al., 2005), which has been successfully applied in several prior LES studies of wind turbine wake flows (Calaf et al., 2010; Stevens and Meneveau, 2017; Martinez, 2017; Stevens et al., 2018; Narasimhan et al., 2022; Gharaati et al., 2022; Narasimhan et al., 2024a; Gharaati et al., 2024). In Eq. (3), the term ${\mathrm{\Pi }}_j=\stackrel{\mathrm{̃}}{u_j\mathit{\theta }}-{\stackrel{\mathrm{̃}}{u}}_j\stackrel{\mathrm{̃}}{\mathit{\theta }}$ is the SGS heat flux whose eddy diffusivity (κ_SGS) is determined from ${\mathit{\kappa }}_{\text{SGS}}={\mathit{Pr}}_{\text{SGS}}^{-\mathrm{1}}{\mathit{\nu }}_{\text{SGS}}$, where the SGS Prandtl number of Pr_SGS=1 (Narasimhan et al., 2022) is prescribed.

The atmospheric boundary layer flow is driven by a geostrophic wind whose pressure gradient is given by $-\mathbf{\nabla }{P}_{\mathrm{\infty }}/\mathit{\rho }=(f_{\mathrm{c}}{V}_{\mathrm{g}},-f_{\mathrm{c}}{U}_{\mathrm{g}})$. Here, $f_{\mathrm{c}}=\mathrm{2}\mathrm{\Omega }\sin \mathit{\varphi }={\mathrm{10}}^{-\mathrm{4}}{\text{s}}^{-\mathrm{1}}$ is the Coriolis parameter corresponding to a mid-latitude position (specifically to ϕ=43.44° with Earth's rotation rate $\mathrm{\Omega }=\mathrm{7.27}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{rad}{\mathrm{s}}^{-\mathrm{1}}$). The quantities U_g,V_g are the geostrophic wind velocity components along the x and y directions, respectively, with magnitude $G=\sqrt{U_{\mathrm{g}}^{\mathrm{2}}+V_{\mathrm{g}}^{\mathrm{2}}}$, and directed at an angle of α_G relative to the x direction such that U_g=Gcos α_G, V_g=Gsin α_G. At each time step, a proportional–integral (PI) controller is utilized to control the direction of the geostrophic wind such that the wind flows in the streamwise direction with zero wind veer at the hub height (Sescu and Meneveau, 2014; Narasimhan et al., 2022).

LESGO uses a Fourier-series-based pseudo-spectral method based on collocated grids for the spatial discretizations in the horizontal (x and y) directions and a second-order central difference method based on staggered grids in the vertical (z) direction. The $\mathrm{3}/\mathrm{2}$ rule is used to eliminate the aliasing error associated with the pseudo-spectral discretization of the nonlinear convective terms. The simulation is advanced in time using a fractional-step method. First, the velocity field is advanced in time by integrating Eq. (2) using the second-order Adams–Bashforth scheme to obtain a predicted velocity field. Then, a pressure Poisson equation is constructed based on the divergence-free constraint Eq. (1) for the new time step and is solved to obtain the pseudo-pressure field. Lastly, the predicted velocity field is projected to the divergence-free space using the gradient of the pseudo-pressure to obtain the velocity field for the new time step. The above fractional steps are repeated at every time step in LES to advance the flow field in time. More details of the numerical schemes used in the LESGO solver can be found in the original references (Albertson, 1996; Albertson and Parlange, 1999).

2.2 Boundary conditions

In the streamwise (x) direction, inflow–outflow boundary conditions are applied using the concurrent precursor simulation approach (Stevens et al., 2014). Specifically, a separate precursor domain without wind turbines is simulated to generate realistic turbulent inflow conditions, which are then imposed at the inlet of the wind farm domain. To ensure periodicity, a fringe region is introduced at the end of the wind farm domain where the outflow is gradually forced to match the inflow from the mapped region in the precursor domain. More details of the inflow–outflow conditions implemented in the current pseudo-spectral solver are provided in Stevens et al. (2014). Additionally, the simulation in the precursor domain uses a shifted periodic boundary condition where the flow field in a spanwise shifting region is shifted to prevent persistent spanwise locking of large-scale turbulent structures (Munters et al., 2016). Following the recommendation in Munters et al. (2016) a shift of $L_{y-\text{shift}}=\mathrm{0.25}{L}_z$ is used in this study, where L_z is the domain height. In the spanwise (y) direction, periodic boundary conditions are used. In the vertical (z) direction, the ground surface boundary condition is specified in both the precursor and wind turbine domains using the Monin–Obukov Similarity Theory (MOST)-based equilibrium surface flux modeling (Monin and Obukhov, 1954). The components of local surface shear stress are computed as a function of the prescribed roughness length according to

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\mathit{\tau }}_{\text{i,3|surf}}=-u_{*}^{\mathrm{2}}{\displaystyle \frac{\widehat{\stackrel{\mathrm{̃}}{u_i}}}{\sqrt{{\widehat{\stackrel{\mathrm{̃}}{u}}}^{\mathrm{2}}+\widehat{{\stackrel{\mathrm{̃}}{v}}^{\mathrm{2}}}}}},i=\mathrm{1},\mathrm{2};\\ & \text{and }{u}_{*}=\mathit{\kappa }{\displaystyle \frac{\sqrt{{\widehat{\stackrel{\mathrm{̃}}{u}}}^{\mathrm{2}}\left(\mathrm{0.5}\mathrm{\Delta }z\right)+{\widehat{\stackrel{\mathrm{̃}}{v}}}^{\mathrm{2}}\left(\mathrm{0.5}\mathrm{\Delta }z\right)}}{\ln (\mathrm{0.5}\mathrm{\Delta }z/z_{\mathrm{0}})}}.\end{array}\end{array}$$

Here, κ=0.41 is the von Kármán constant, and z₀ is the prescribed roughness length. The friction velocity u_* is expressed in terms of the horizontal velocity $(\widehat{\stackrel{\mathrm{̃}}{u}},\widehat{\stackrel{\mathrm{̃}}{v}})$ at the first grid point (z₁=0.5Δz), filtered at twice the grid resolution, $\widehat{\stackrel{\mathrm{̃}}{\mathrm{\Delta }}}=\mathrm{2}\stackrel{\mathrm{̃}}{\mathrm{\Delta }}$ (Bou-Zeid et al., 2005). Since we simulate conventionally neutral conditions, the surface heat flux is set to zero; thus no stability correction terms (as used in Xiao et al., 2025) are included. At the top of the domain, a stress-free boundary condition is imposed. A sponge or Rayleigh-damping layer (Durran and Klemp, 1983) is included approaching the top boundary, ranging from 0.75L_z to L_z, with a sponge inverse relaxation timescale (frequency) parameter of $\mathrm{3.9}\times {\mathrm{10}}^{-\mathrm{3}}{\text{s}}^{-\mathrm{1}}$. In this layer, a damping body force with a cosine profile is applied to suppress the reflection of gravity waves.

Henceforth, the $\stackrel{\mathrm{̃}}{(\cdot )}$ notation for LES-filtered field variables (e.g., velocity ${\stackrel{\mathrm{̃}}{u}}_i$, temperature $\stackrel{\mathrm{̃}}{\mathit{\theta }}$) will be omitted for brevity. All subsequent variables should be interpreted as implicitly filtered quantities obtained from the LES solution, governed by the equations presented in this section.

2.3 Wind turbine representation

The aerodynamic forces exerted by wind turbines on the airflow are modeled through the distributed body force term f_i in the momentum transport equations (Eq. 2). During the initial spin-up phase (i.e., Phase 1), we employ an actuator disk model (ADM) on a coarse grid for computational efficiency, with the thrust force magnitude calculated as $f=\frac{\mathit{\pi }}{\mathrm{8}}\mathit{\rho }{C}_T^{\prime }\langle u_T{\rangle }_{\mathrm{d}}^{\mathrm{2}}{D}^{\mathrm{2}}$ (Calaf et al., 2010; Howland et al., 2016). Here, ρ is the air density, 〈u_T〉_d is the local wind velocity averaged over the rotor disk, D is the diameter of the wind turbine, and $C_T^{\prime }$ is the local thrust coefficient (set to a common value $C_T^{\prime }=\mathrm{1.33}$). We recall that $C_T^{\prime }$ is based on the disk-averaged velocity 〈u_T〉_d, which, unlike the far-upstream velocity U_∞, is immediately available in LES (Calaf et al., 2010).

After the spin-up simulation converges to quasi-steady behavior, the grid is refined to its final resolution, and the actuator line model (ALM) is adopted (Sørensen and Shen, 2002; Troldborg, 2009; Jha et al., 2014; Martínez-Tossas et al., 2015). In ALM, each turbine blade is represented by a collection of actuator points along a line, where forces are applied according to the velocity field and the angle of attack. The forces per unit width at every actuator point are computed as

$$\begin{array}{}\text{(5)}& {\mathbf{f}}_{\text{alm}}=\mathrm{0.5}\mathit{\rho }c\left|{\mathbf{V}}_{\text{rel}}{|}^{\mathrm{2}}\right(C_{\mathrm{L}}{\mathit{e}}_{\mathrm{L}}+C_{\mathrm{D}}{\mathit{e}}_{\mathrm{D}}),\end{array}$$

where c is the airfoil chord length, $\left|{\mathbf{V}}_{\text{rel}}\right|$ is the magnitude of the relative velocity of the upwind flow to the turbine blade, C_L and C_D are lift and drag coefficients obtained from tabulated airfoil data, and e_L and e_D are unit vectors along the direction of the lift and drag forces at each actuator point, respectively. These forces are then smeared using a Gaussian kernel to project them into the computational LES grid:

$$\begin{array}{}\text{(6)}& {\mathit{\eta }}_{\mathit{ϵ}}={\displaystyle \frac{\mathrm{1}}{{\mathit{ϵ}}^{\mathrm{3}}{\mathit{\pi }}^{\mathrm{3}/\mathrm{2}}}}{e}^{-r^{\mathrm{2}}/{\mathit{ϵ}}^{\mathrm{2}}},\end{array}$$

where r is the distance from the grid point to the actuator point and ϵ denotes the width of the kernel. The kernel width is chosen to be at least $\mathit{ϵ}=\mathrm{2}({\mathrm{\Delta }}_x{\mathrm{\Delta }}_y{\mathrm{\Delta }}_z{)}^{\mathrm{1}/\mathrm{3}}$, as recommended to avoid numerical instabilities (Troldborg, 2009; Martínez-Tossas et al., 2015).

The accuracy of the ALM can be sensitive to grid resolution and the choice of ϵ. The optimal ϵ_opt needed to resolve the induced velocities is typically much smaller than the ϵ used to avoid numerical instabilities (Martínez-Tossas et al., 2017). To address this challenge, we use the generalized filtered lifting line theory correction to accurately represent the blade aerodynamics (Martínez-Tossas and Meneveau, 2019; Martínez-Tossas et al., 2024), including the shedding of unresolved vorticity leading to missing induced velocities at the blade. The correction accounts for subgrid-scale induced velocity that would be obtained by using an optimal ϵ_opt by estimating its contribution and adding it to the resolved velocity in the LES. With the correction, the ALM provides consistent blade loading predictions across varying grid resolutions.

The NREL-5MW baseline wind turbine (Jonkman et al., 2009) is adopted as our reference model. It is a widely used benchmark model developed by NREL to standardize research on wind technologies. The turbine has a diameter of D=126 m, three blades, and a hub height at elevation z_h=90 m. It reaches a rated electrical power output of 5 MW at a rated wind speed of approximately 11.4 m s⁻¹. Its rotor blades utilize the Delft University (DU) and National Advisory Committee for Aeronautics (NACA) series airfoil profiles optimized for aerodynamic efficiency, structural integrity, and minimal fatigue loads, making the NREL-5MW turbine an essential tool for evaluating wind turbine performance, control strategies, structural design, and off-shore platform dynamics.

The dataset employs fixed but row-dependent rotor angular velocities determined through an initialization procedure. Initialization begins with all turbines operating at tip-speed ratio TSR = 7.5 (near-optimal for NREL-5MW turbines). In this initialization simulation (i.e., first part of Phase 2), the angular velocity Ω for each turbine is then computed dynamically using

$$\begin{array}{}\text{(7)}& \mathrm{\Omega }=\text{TSR}\times {\displaystyle \frac{\mathrm{1.087}{U}_{\mathrm{d}}}{(\mathrm{1}-a)R}},\end{array}$$

where U_d is the disk-averaged velocity; the numerator incorporates an empirical 8.7 % correction factor for LES filter-scale effects (ϵ=16 m), validated through single-turbine laminar inflow tests; and the induction factor a derives from rotor geometry (blade number N_b=3, radius R=63 m, and chord c=3–4 m) and local inflow angle ϕ via

$$\begin{array}{}\text{(8)}& a={\displaystyle \frac{\mathrm{1}}{\left(\mathrm{4}{\sin }^{\mathrm{2}}\mathit{\varphi }\right)/\left({\mathit{\sigma }}_{\mathrm{r}}{C}_n\right)+\mathrm{1}}},\end{array}$$

with a rotor solidity of ${\mathit{\sigma }}_{\mathrm{r}}=N_{\mathrm{b}}c/\left(\mathit{\pi }R\right)$ and a force coefficient of $C_n=C_{\mathrm{L}}\cos \mathit{\varphi }+C_{\mathrm{D}}\sin \mathit{\varphi }$. After approximately 40 min (minutes) of initialization simulation, the angular velocity Ω for each turbine is averaged within its respective row, which serves as the fixed operational values for the subsequent database simulations.

We also note that LESGO's ALM implementation includes detailed turbine operation control methods, such as pitching the blades (feathering) during region 3 operations. In the current simulation we chose to operate all turbines exclusively at optimal tip-speed ratio, “region 2” (also without including regions 1.5 and 2.5). This choice was made in order to avoid the need to store additional data relating to blade pitch (curtailment) and other complex turbine control actions. Since this practice deviates slightly from the reference NREL-5MW nameplate data, we refer to the turbine in our simulations as the NREL-5MW+ turbine. Specifically, the front turbines are allowed to rotate slightly faster than the maximum rotation rate of the original NREL-5MW reference turbine.

Figure 1 Schematic representation of the computational simulation domain (not to scale), showing (a) top view (x–y plane), (b) side view (x–z plane), and (c) front view (y–z plane). The precursor computational domain consists of the regions denoted as “P”, the precursor mapping region “P_M”, and the precursor spanwise shifting region “P_S”. The wind farm computational domain includes the wind farm region “W” and the fringe region “W_F” near the outlet. Both the precursor and wind farm computational domains include a Rayleigh damping region at the top (denoted as“R”). The turbine diameter D=126 m and hub height z_h=90 m are also marked.

3 Simulation parameters

We simulate turbulent flow through a 10×6 array of NREL-5MW+ turbines (with diameter D=126 m) in a $\mathrm{28.224}\times \mathrm{3.78}\times \mathrm{2}{\mathrm{km}}^{\mathrm{3}}$ domain, equally split between precursor and wind farm subdomains (each 112D=14.112 km long). Figure 1 displays the domain dimensions. The precursor domain includes the region denoted as P of length $\mathrm{5}{L}_x/\mathrm{8}$, the mapping region P_M of length $L_x/\mathrm{8}$, and the spanwise shifting region P_S of length $L_x/\mathrm{8}$. The wind farm domain features 14D of upstream buffer zone, a 63D turbine region, a 21D downstream wake recovery region (these three regions combined are denoted as W), and a 14D outflow fringe region (W_F). The turbines are spaced 7D (streamwise) and 5D (spanwise), with lateral boundaries 2.5D from the outermost turbines. Note that P_M, P_S and W_F regions have a length of $L_x/\mathrm{8}$ and that P_M extends from $\mathrm{5}{L}_x/\mathrm{8}$ to $\mathrm{3}{L}_x/\mathrm{4}$. Vertically, a 0.5 km Rayleigh damping sponge layer (denoted as R) is located between 1.5 and 2 km (see Fig. 1). We adopt θ₀=263.5 K as the reference potential temperature, consistent with the value chosen in studies by Gadde and Stevens (2021) and our prior simulations of stable boundary layer (SBL) and CNBL flows reported in Narasimhan et al. (2024a). This reference temperature was inspired by observations from the Beaufort Sea Arctic Stratus Experiment (BASE) and simulations by Kosović and Curry (2000). While the value of θ₀ is relatively low, it serves primarily as a relative additive reference that does not significantly affect the simulated flow dynamics or the physical interpretation of the results. For example, if we used 273K, it would change the implied thermal expansion coefficient in our Boussinesq approximation only by about 3 %.

The turbulent flow is driven by a constant geostrophic wind speed $G=\mathrm{15}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ at ${\mathit{\alpha }}_{\mathrm{g}}\approx -\mathrm{22.5}$° to the x direction, with the angle controlled by a PI controller (K_P=10, K_I=0.5) to align hub-height mean wind velocity with the x axis in the conventionally neutral boundary layer (Sescu and Meneveau, 2014; Narasimhan et al., 2022). The surface has roughness length of z₀=0.1 m and a reference potential temperature of θ₀=263.5 K. Initial conditions are set at $U_{\mathrm{g}}=\mathrm{15}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ (streamwise) and $V_{\mathrm{g}}=\mathrm{0}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ (spanwise), perturbed by random noise, while potential temperature decreases from 265 K at the surface with a 1 K km⁻¹ lapse rate, including random perturbations below 1 km.

The numerical simulation is conducted in three consecutive phases to ensure proper flow development and statistical convergence.

• Phase 1: coarse-resolution ADM spin-up. A 10 h simulation using the ADM is performed to establish a quasi-stationary atmospheric boundary layer and wind farm wake field. This phase leverages the computational efficiency of ADM, which approximates turbine forces without resolving actuator line-level aerodynamics.

• Phase 2: fine-resolution ALM convergence. A 1 h simulation using the actuator line model at finer spatial resolution transitions the flow from ADM-averaged to ALM-resolved turbine representation. Besides the turbine model update, two additional changes are introduced in this phase: (i) the time-stepping scheme is switched from a constant Courant–Friedrichs–Lewy (CFL) number of 0.0625 to a fixed time step of Δt=0.025 s. This adjustment has negligible impact on the results because, under these simulation conditions, CFL=0.0625 corresponds to Δt≈0.03 s. The slightly more restrictive Δt=0.025 s maintains numerical stability while preserving solution accuracy. (ii) The rotor control changes from a fixed tip-speed ratio (TSR = 7.5) to fixed rotor angular velocities that vary across turbine rows, as tabulated in Table 1. This adjustment has a negligible impact on the results because the prescribed angular velocities closely match the values achieved under TSR=7.5 conditions (see the calculation method in Sect. 2.3), ensuring nearly identical rotor dynamics.

• Phase 3: fine-resolution simulation for database construction. A final 1 h simulation is carried out to collect high-fidelity flow and turbine data. Flow field variables are recorded every 20 LES time steps (i.e., every 0.5 s) on a filtered and subsampled spatial grid (every other grid point in the x–y plane), while wind turbine data – both integral and blade-resolved – are stored at every LES time step (0.025 s). Note that we purposefully operate the NREL-5MW+ turbine in “region 2” during the simulation time, in order to avoid having to choose and document additional controller actions. As a result, during some times, some of the turbines operate “above rated conditions” but maintain self-consistent aerodynamic behavior of the blades and air flow.

Table 1 Rotor speed for each row of turbines.

The three phases of the simulation are illustrated through the time history of the boundary layer height z_i=h_ABL and the geostrophic wind angle shown in Fig. 2.

Table 2 Three consecutive phases and computational domain parameters

Figure 2 Time history of boundary layer height z_i=h_ABL and geostrophic wind angle α, indicating the three simulation phases (Phase 1: coarse-resolution ADM spin-up; Phase 2: fine-resolution ALM convergence; and Phase 3: fine-resolution simulation for database construction).

4 JHTDB-wind database construction

The LES data from the final 1 h sampling period are systematically ingested into the database and organized into two primary data types: (i) flow field data, consisting of 4D space–time fields captured across both simulation domains (precursor and wind farm domains), providing complete spatiotemporal information about the atmospheric flow, and (ii) turbine data, which are further subdivided into two subtypes. The first subtype is turbine-level operational data, comprising time histories of turbine power and thrust. The second subtype is blade-level data, which include time histories of aerodynamic quantities sampled at each discrete actuator point along each blade.

4.1 Flow field data

4.1.1 Domain of the dataset

As described in Sect. 3, the LES is conducted in the domain of dimensions $(\mathrm{2}\times \mathrm{14.112})\times \mathrm{3.78}\times \mathrm{2}{\mathrm{km}}^{\mathrm{3}}$ (see Table 2). When compiling the database, we exclude numerically imposed auxiliary regions: specifically, the final $L_x/\mathrm{4}$ of the precursor domain (which includes the precursor spanwise shifting region P_S) and the final $L_x/\mathrm{8}$ of the wind farm domain (i.e., the wind farm fringe region W_F), as visualized in Fig. 1. These regions serve purely numerical functions (periodicity enforcement and inflow recycling, respectively) without contributing to physical flow dynamics of interest. The resulting database domain has the extents of $(\mathrm{10.584}+\mathrm{12.348})\times \mathrm{3.78}\times \mathrm{2}{\mathrm{km}}^{\mathrm{3}}$, as shown in Fig. 3. The top 0.5 km sponge region is kept in the database for simplicity of data management and possible interest.

Figure 3 Schematic representation of the database domain (not to scale). This is the physical domain available in the database, merging the precursor domain (P+P_M) up to the end of the mapping region at $\mathrm{3}/\mathrm{4}{L}_x$ with the wind farm domain (W) and excluding the fringe region (W_F). A total of 60 turbines are shown, with only a subset labeled for clarity. The domain dimensions are $(\mathrm{10.584}+\mathrm{12.348})\times \mathrm{3.78}\times \mathrm{2}{\mathrm{km}}^{\mathrm{3}}$.

4.1.2 Spatial resolution of the dataset

To minimize storage, we applied spectral filtering on x–y planes for flow field data by truncating Fourier modes above ${\mathit{\kappa }}_{\text{max}}/\mathrm{2}$, where ${\mathit{\kappa }}_{\text{max}}=\mathit{\pi }/{\mathrm{\Delta }}_{\text{LES}}$ is the LES cutoff wavenumber. The filtered fields were then subsampled at every alternate grid point in the x and y directions, maintaining the original vertical (z) resolution. This approach reduces the dataset size by 75% while maintaining fidelity in capturing the dynamically significant larger-scale flow structures and turbine wake interactions. Thus, the flow field data have a grid size of $(\mathrm{576}+\mathrm{672})\times \mathrm{192}\times \mathrm{400}$.

4.1.3 Temporal resolution of the dataset

Field data are stored at intervals of 0.5 s (i.e., every 20 LES steps of 0.025 s), ensuring that fluid parcels advected at the maximum geostrophic speed (15 m s⁻¹) travel less than the horizontal grid spacing (Δx≈9.19 m) between snapshots. Although rotor blade tips move across several vertical grid spacings during this interval, the corresponding rotor force field is smooth (Gaussian filtered at scale $\mathit{ϵ}=\mathrm{16}\mathrm{m}>\mathrm{2}\sqrt[\mathrm{3}]{\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z}$), ensuring that the storage frequency of 0.5 s remains appropriate. Over the 1 h simulation period (i.e., 3600 s), the simulation advances through $\mathrm{3600}/\mathrm{0.025}=\mathrm{144}\mathrm{000}$ LES time steps, with flow fields stored at $\mathrm{144}\mathrm{000}/\mathrm{20}=\mathrm{7200}$ consecutive snapshots.

4.1.4 Final structure of the dataset

Consequently, the final stored data dimensions are $n_x\times n_y\times n_z\times n_{\mathrm{t}}=\mathrm{1248}\times \mathrm{192}\times \mathrm{400}\times \mathrm{7200}$. At each stored time step, six spatial fields are recorded: the three velocity components $u(x,y,z,t)$, $v(x,y,z,t)$, and $w(x,y,z,t)$; the (kinematic) pressure field $p(x,y,z,t)/\mathit{\rho }=p^{*}(x,y,z,t)-u_k{u}_k/\mathrm{2}$ (the SGS stress trace is not available and is anyhow negligible); the potential temperature field relative to the reference temperature ${\mathit{\theta }}^{\prime }(x,y,z,t)=\mathit{\theta }(x,y,z,t)-{\mathit{\theta }}_{\mathrm{0}}$; and the subgrid-scale eddy viscosity ${\mathit{\nu }}_{\text{SGS}}(x,y,z,t)$. In addition, the three components of the turbine force field, $f_x(x,y,z,t)$, $f_y(x,y,z,t)$, and $f_z(x,y,z,t)$, are also stored. Unlike the other flow field variables, these force components are stored only from the ground up to 200 m in the vertical direction. However, they are retained at the original spatial resolution (i.e., not filtered in the x–y planes). The detailed information of these stored field variables can be found in Table 3. It also needs to be mentioned that the concurrent precursor method ensures smooth transitions in velocity, potential temperature, and eddy viscosity fields between precursor and wind farm subdomains, by construction. However, due to the non-local nature of the pressure solution (solved separately in each domain via Poisson equations) and the velocity-only coupling between domains, the stored pressure field exhibits a minor discontinuity at the interface. This artifact does not affect the resolved turbulence dynamics or turbine wake interactions but needs to be taken into account if computing pressure gradients across the boundary separating the precursor and wind farm domains.

Table 3 Summary of flow field data variables.

These 4D field variables are stored using Zarr format (Miles et al., 2023). In Zarr-based storage, data are organized into chunks, the smallest units retrieved during a query. To ensure efficient data access, chunk sizes must be large enough to support common operations, such as differentiations and interpolations, which typically require access to a 3D neighborhood around the query point, while remaining small enough to avoid excessive memory usage. Based on extensive testing and prior experience with other JHTDB datasets, a chunk size of 64³ grid points provides optimal retrieval speeds and performance for typical data access modalities. We chose a similar chunk size but shaped according to $\mathrm{52}\times \mathrm{64}\times \mathrm{80}$ so that an integer multiple of the chunk size in each direction fits into the stored domain size. The total amount of data stored is about 15 TB. These flow field data can be queried using getData(...) calls from analysis programs such as Python, MATLAB, Fortran, or C, in the same manner as with other turbulence datasets available through JHTDB.

4.2 Wind turbine data

4.2.1 Turbine-level data

The turbine-level data are integral quantities characteristic of each turbine operation, which are derived from the actuator line modeling. This dataset includes high-fidelity time histories of power output, thrust force, and rotor angular velocity, sampled at Δt=0.025 s for all 60 turbines, as summarized in Table 4. In the present dataset, the angular velocity is held constant in time, but for other datasets (e.g., Xiao et al., 2025), this is not generally the case. For each variable, the dataset consists of 144 000 rows and 2 columns, where the first column represents time and the second column contains the corresponding values of the recorded variable. The turbine data are stored in files using Parquet format, which facilitates efficient access and querying from various programming languages. Turbine-level data can be accessed using the getTurbineData(...) function call from analysis environments such as Python or MATLAB.

Table 4 Summary of turbine-level data variables. Each dataset is a 2D matrix of size n_t×2, where n_t is the number of time steps. Columns 1 and 2 represent time and measured values, respectively.

Table 4 summarizes the turbine-level data variables. Note that, unlike the field data which are stored in kinematic (density-independent) units, the force and power data require a specified air density. The value used in the simulations to compute these forces is ρ_air=1.23 kg m⁻³.

4.2.2 Blade-level data

In addition to the integral quantities characteristic of each turbine's operation, more detailed information is captured along each turbine blade to enable blade-resolved aerodynamic analysis. This fine-grained dataset allows users to investigate the local aerodynamic behavior of blades under unsteady flow conditions, which is critical for understanding load distributions, fatigue effects, and control optimization strategies. The turbine blade-level dataset includes high-fidelity time histories sampled at 0.025 s for all 180 blades in the wind farm (i.e., 60 turbines × 3 blades each), with aerodynamic and geometric quantities sampled at 100 discrete actuator line points along the blade span. As summarized in Table 5, a total of 19 variables are sampled and stored, with each variable written to a separate file. For each variable, the dataset has dimensions of 144 000×3 rows and 103 columns. Each time step includes three rows corresponding to the three blades of a turbine, resulting in a total of 144 000×3 rows. Vertically, the first column represents time in seconds, the second column specifies the turbine number, and the third column denotes the blade number (blades can be identified by the time histories of the individual ALM point positions). The remaining 100 columns contain the values of the selected variables at each of the 100 actuator points from the blade root to tip. Similarly to turbine-level data, blade-level data are stored as Parquet files, allowing efficient access across multiple programming environments. Blade-level data can be accessed using the getBladeData(...) function call from analysis environments such as Python or MATLAB.

Table 5 Summary of blade-level data variables. Each dataset is a 2D matrix of size $(n_{\mathrm{t}}\times \mathrm{3})\times \mathrm{103}$. Here, (n_t×3) represents the total number of blade-wise samples, formed by concatenating the time series data from each of the three blades of a turbine. Columns 1–3 represent time and turbine number, and columns 4–103 store aerodynamic measurements at n_ℓ=100 discrete locations along each blade.

Figure 4 Contour plots of instantaneous flow field variables in part of the precursor domain (here between x=0 m and x=10 381.875 m), at time t=1800.75 s. (a) The streamwise velocity u, (b) the vertical velocity w, (c) the pressure p, and (d) the potential temperature deviation θ^′.

5 Web-accessible virtual sensor data access methods and examples

5.1 Flow field data

A defining feature of the JHTDB database system (Li et al., 2008) is its low entry barrier for data usage, enabling users to efficiently explore large-scale simulation datasets through web services and virtual sensor methodology. The JHTDB-wind system adopts the same approach, allowing access to wind farm data using these established tools. Users can develop analysis scripts or notebooks in familiar programming languages such as Python and MATLAB (as well Fortran and C) to run them remotely on their own machines or on SciServer, a cloud service dedicated to running code close to the data. Within these analysis environments, users specify space–time arrays by defining spatial locations (e.g., along a line, across a surface, within a subvolume, or scattered arbitrarily) and corresponding time instances; i.e., users specify the positions of virtual sensor arrays. These space–time arrays are then passed to the pre-defined function, getData(...), which returns interpolated values of the selected variables at defined coordinates. This framework enables targeted, on-demand data access without the need to download large volumes of raw simulation output.

Figure 5 Contour plots of instantaneous flow field variables in part of the wind farm domain (here between x=10 584 m and x=21 921.375 m), at time t=1800.75 s. (a) The streamwise velocity u, (b) the vertical velocity w, (c) the pressure p, and (d) the potential temperature deviation θ^′. The short black lines represent the location of wind turbines.

Figure 6 Python code snippet used to obtain the data to generate Fig. 5.

Figure 7 Contour plot of instantaneous streamwise velocity u in the entire database domain, ranging from x=0 to x=22 913.625 m, at time t=2505 s. It is noted that, although the total length of the database domain is $\mathrm{10}\mathrm{584}+\mathrm{12}\mathrm{348}=\mathrm{22}\mathrm{932}$ m, the data resolution in the x direction is 18.375 m and the grid points are located at cell centers. Consequently, the last data point is located at $\mathrm{22}\mathrm{932}-\mathrm{18.375}=\mathrm{22}\mathrm{913.625}$ m. The short black lines represent the location of wind turbines.

Figure 8 Instantaneous contours of turbine streamwise (i.e., x-component) force (as projected onto the LES grid using Gaussian smoothing as part of the ALM method) in the y–z planes (a) in the first row (i.e., Row 1, x=12 348 m) and between the relevant vertical range $z\in [\mathrm{2.5},\mathrm{200}]\text{m}$ and (c) in the second-to-last row (i.e., Row 9, x=19 404 m). Panels (b) and (d) show the x-direction pressure gradient distributions on the same planes, coincident with the turbines.

Figure 9 Python code snippet used to obtain the data to generate Fig. 8.

Figure 10 Vertical profiles of horizontal- and time-averaged (a) velocities $\langle u\left(z\right){\rangle }_{x,y,t}$, $\langle v\left(z\right){\rangle }_{x,y,t}$ and velocity magnitude $\mathcal{V}(z{)}_{x,y,t}={\left[\langle u(z{)}_{x,y,t}{\rangle }^{\mathrm{2}}+\langle v(z{)}_{x,y,t}{\rangle }^{\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}}$, (b, bottom axis) subgrid-scale eddy viscosity $\langle {\mathit{\nu }}_{\text{SGS}}\left(z\right){\rangle }_{x,y,t}$ used in the LES as a result of the Lagrangian scale-dependent dynamic model, and (b, top axis) potential temperature deviation $\langle {\mathit{\theta }}^{\prime }\left(z\right){\rangle }_{x,y,t}$ (i.e., the deviations from a reference temperature θ₀=263.5 K).

Figure 11 Python code snippet used to obtain the data to generate vertical profiles of $\langle u\left(z\right){\rangle }_{x,y,t}$: for the 100 heights z between z=2.5 m and z=1997.5 m, we query data on a regular mesh (not necessarily coinciding with stored grid points). For statistical convergence, we average over four times covering the entire hour (t=900 s, 1800 s, 2700 s, 3598.5 s).

Figures 4 and 5 display contour plots of flow field variables at the turbine hub height ($z=z_{\mathrm{h}}=\mathrm{90}\text{m}$) for the precursor and wind farm domains, respectively. Figure 6 presents Python code snippets that demonstrate how to query the JHTDB-wind database to extract snapshots of velocity, pressure, and potential temperature fields at a specific time, approximately in the middle of the stored 1 h dataset, namely at t=1800.75 s. As a first step, an array “points” is populated with spatial coordinates that define a 2D plane: in this case, an equally spaced grid of 950×200 points in the x and y directions at a constant height $z=z_{\mathrm{h}}=\mathrm{90}\text{m}$. These query points typically do not coincide with the actual simulation grid points, and users are not required to know the grid layout to access the data. The JHTDB-wind interface provides interpolated field values based on a user-specified interpolation method. Supported options include no interpolation (it returns the value at the nearest grid point); Lagrange polynomials of order 4, 6, or 8; and several spline interpolation methods (Li et al., 2008; Graham et al., 2016). In this example, we use eighth-order Lagrange polynomial interpolation in space. Similarly, if the requested time does not coincide with a stored time step, temporal interpolation is applied using the third-order Piecewise Cubic Hermite Interpolating Polynomial (PCHIP) method (Li et al., 2008). This user-friendly data access model eliminates the need for downloading and parsing simulation files. Instead, the Python application programming interface (API) returns arrays with the queried field variables, which can then be visualized directly within a Jupyter notebook (or MATLAB code). This approach was used to generate Figs. 4 and 5. It is important to note that the full 1 h dataset (comprising 14 400 time steps) is available for analysis, allowing users to query any time between t=0 and t=3600 s. For example, Fig. 7 shows a hub-height snapshot over the entire domain at time t=2505 s.

Similar queries can be made for the values, spatial gradients, and Hessians (second-order derivatives) of all variables listed in Table 3. For example, Fig. 8a and b show the turbine streamwise force field f_x and the x-direction gradient of the pressure field ($\partial p/\partial x$), respectively, on a y–z plane intersecting Row 1 (turbines #1–#6) at x=12 348 m (1764 m downstream of the wind farm domain), at time t=1000.013 s. Figure 8c and d present similar results on a plane intersecting Row 9 (turbines #49–#54) at x=19 404 m (8820 m downstream of the wind farm domain) at another time t=2000.67 s. These plots were generated using the Python code shown in Fig. 9. In these examples, the queried times are intentionally chosen not to coincide with the stored simulation time steps, demonstrating the temporal interpolation capabilities of JHTDB-wind.

Next, we provide examples of computed mean vertical profiles of fundamental flow quantities within the precursor domain, which features standard conventionally neutral atmospheric conditions. Figure 10 shows vertical profiles of horizontal- and time-averaged mean velocities, subgrid-scale eddy viscosity, and deviations in potential temperature, all obtained by averaging in the horizontal directions and over time. The data used to produce these profiles are retrieved using the virtual sensor framework, and an example code snippet demonstrating this process is shown in Fig. 11.

Figure 12 Time evolution of power from turbines during the 10 min time interval, i.e., $t\in [\mathrm{1000.33},\mathrm{1600.33}]\text{s}$. Panel (a) shows the total power from the entire wind farm, while panel (b) shows the power for the turbines in Row 1 (i.e., turbines #1–#6) and in Row 9 (i.e., turbines #49–#54).

Figure 13 Python code snippet illustrating the use of the function getTurbineData(...) as part of a loop over all turbines in the wind farm and subsequent summation to evaluate time series of total power used to generate Fig. 12a.

5.2 Wind turbine data

Wind turbine data, including both the turbine-level and blade-level data, are considerably smaller than the 4D flow field data, and one possibility would have been to allow users to download these data directly as files. However, such an approach would require users to identify specific files, understand naming conventions, and handle formatting, posing a barrier to seamless integration with flow field queries. To maintain consistency and usability across the platform, we adopt a similar virtual sensor data access paradigm used for the flow field data. Two dedicated query functions are developed: getTurbineData(...) for turbine-level quantities and getBladeData(...) for blade-resolved data. For getTurbineData (...), users specify the turbine number (ranging from 1 to 60) and desired time instances. For getBladeData (...), both turbine number and blade number (ranging from 1 to 3) need to be specified, along with an array of actuator point indices (ranging from 1 to 100) and times (ranging from 1 to 3600 s) at which the data are requested. Linear interpolation in time is supported to provide values between stored simulation time steps.

As an example, Fig. 12 presents the time series of total wind farm power output (panel a) and of rows 1 and 9 of six turbines (panel b). The code snippet specifying the getTurbineData(...) call is shown in Fig. 13. Similar calls can be made to extract any of the turbine-specific variables listed in Table 4.

Next, we illustrate the use of getBladeData(...) in Fig. 14, which shows (a) the time evolution of the lift and drag coefficients and (b) the lift coefficient as a function of blade angle. The blade angle is computed as $\mathit{\zeta }\left(t\right)=\arctan \left[z\right(t)-z_h)/\left(x\right(t)-x_T)]$ over a 60 s period. The results shown are for a particular turbine and blade (Turbine #28 in the central portion of the wind farm and Blade #3 – the latter being an arbitrary choice, of course). The Python code snippet shown in Fig. 15 illustrates the use of getBladeData(...), with the queried data plotted directly as a time series within the same script. Using a similar approach, variable data can be extracted along turbine blades and further processed to compute higher-order statistics. Figure 16 shows axial force, tangential force, and drag and lift coefficients for an upstream turbine (Blade #1 of Turbine #1) and a downstream turbine (Blade #1 of Turbine #60) at a specific time of t=1500 s. Any of the variables listed in Table 5 can be similarly queried (also in MATLAB).

Figure 14 (a) Time evolution of lift and drag coefficients on an ALM point 80 % along the span of Blade #3 for Turbine # 28. (b) Polar plot of lift coefficient for that point as a function of blade angle along its rotation. For this turbine, the rotational speed is fixed at $\mathrm{\Omega }=\mathrm{1.09}\mathrm{rad}{\mathrm{s}}^{-\mathrm{1}}$ (as obtained from getTurbineData(...)), corresponding to approximately 10.5 revolutions during a 60 s period.

Figure 15 Python code snippet used to obtain the data to generate Fig. 14.

Figure 16 Distributions of ALM quantities along the turbine blade at a specific time (t=1500 s for Blade #1 of Turbine #1, blue lines; Blade # 1 of Turbine #60, orange lines): (a) axial component of the local force (on each Δℓ=0.615 m segment) F_axi, (b) tangential component of the local force (on each Δℓ=0.615 m segment) F_tan, (c) lift coefficient C_l, (d) drag coefficient C_d.

6 Conclusions

In this paper, we have introduced JHTDB-wind, hosting datasets from high-fidelity LES simulations of wind farms. We extend the standard “virtual sensors” data access methods (Li et al., 2008; Yu et al., 2012; Graham et al., 2016) that have been successfully used for democratizing access to more fundamental turbulence datasets. Besides velocity, pressure, potential temperature, and SGS eddy-viscosity fields, JHTDB-wind adds 4D space–time data on aerodynamic turbine force distributions as seen by the flow and the time series of turbine- and actuator-line-specific aerodynamic data along each of the turbine blades, modeled using ALM. We explain the simulation details, provide background on the numerical method and flow parameters, and provide detailed examples and explanations of the user-friendly data access methodologies. It is hoped that these data will provide useful insights about the complex fluid dynamic processes occurring in wind farms.

We realize that in generating a dataset for a representative conventionally neutral boundary layer case, with a relatively large wind farm with 60 turbines, many other choices could have been made (flow parameters, turbine model and control scheme, usage of a particular LES numerical code, numerical resolution, and so on). We anticipate that different members of the community would have made different choices, and we look forward to conversations about how to further improve such datasets. We believe, however, that the case selected is representative of CNBL wind farm dynamics that have been studied by many others before, with a well-tested numerical code. Hence, the authors hope that the data can be of some use and interest to researchers in wind energy.

As a final note, we have additionally prepared a second dataset for JHTDB-wind featuring an eight-turbine wind farm over a full diurnal cycle, capturing both strongly stable and unstable atmospheric boundary layer regimes at different times of the day and night (Xiao et al., 2025).