This is a well-written paper that showcases novel lattice-Boltzmann techniques for atmospheric boundary layer simulations. With clarity, it accomplishes what it sets out to do, which is to demonstrate the accuracy and speed of the implemented approach in the VirtualFluids solver. The discussion is thorough, addressing the theory and formulation, placing these in the greater research context, and connecting these to the usage of computational resources.

There are no clear mistakes or omissions in the progression of the paper, and this "important step for the lattice Boltzmann method" should proceed to publication.

As a minor note, the last sentence at the beginning of section 2.6 is missing a word, as it should read "layer has to be implemented." Also, enlarging the figures that plot curves (Figures 3, 4, and 5) would help with readability and distinguishing between the results from each solver.

Review of the manuscript “Large Eddy Simulation of Thermally Stratified Atmospheric Boundary Layers with a Lattice Boltzmann Method” by Henry Korb, Henrik Asmuth, Martin Schonherr, Martin Geier, and Stefan Ivanell, submitted for publication in Wind Energy Science.

In the manuscript “Large Eddy Simulation of Thermally Stratified Atmospheric Boundary Layers with a Lattice Boltzmann Method” the authors present development and testing of a large-eddy simulation (LES) model for simulation of flows in an atmospheric boundary layer based on the Lattice Boltzmann method. The model is evaluated by simulating a conventionally neutral boundary layer (CNBL) and a stably stratified atmospheric boundary layer and comparing the results to previous studies by Berg et al. (2020), Kuhn et al. (2025), and Gadde et al. (2021).

General Remarks

The manuscript is written well, and it clearly conveys the characteristics of the Lattice Boltzmann model, VirtualFluids, and results of the study. Important new development presented in the manuscript is introduction of lower boundary condition for atmospheric boundary layers based on the Monin-Obukhov similarity theory. The numerical study is based on the design of previous studies. The results of the simulations of a conventionally neutral and a stably stratified atmospheric boundary layer (ABL) with the Lattice Boltzmann model are in general comparable to the results obtained with models that use finite volume/finite element discretization. The flow is, as expected, less well resolved in comparison to pseudo-spectral model.

The development of the Lattice Boltzmann model is motivated by the computational performance, however, limited assessment of model performance is presented. The simulations are executed on one GPU. Since the VirtualFluids code also includes parallel processing capability based on the Message Passing Interface (MPI), it would be important to demonstrate the scaling properties of the model as it is commonly done (e.g. Min et al. 2024, Int. J. HPC App.; Sauer and Muñoz-Esparza, 2020, JAMES). Furthermore, the definition of buoyant forcing in the manuscript is not correct. Considering that the simulations of the stably stratified ABL produced reasonable results it is likely that the Equation (5) where buoyancy force is defined is not correct and not the implementation in the code. There are also several statements that would need to be clarified or made more precise as outlined below under Specific Remarks.

Taking all the above into account the manuscript can be accepted for publication in the journal Wind Energy Science after comments and suggestions for major revisions are addressed.

Specific Remarks

Abstract – line 9 – The statement about real time execution needs to be qualified, since real time execution depends on the grid cell size and wind speed.

Lines 16-18 – The statement starting with “The example…” should be qualified. RANS models may not be able to accurately capture complexity of marine boundary layers, but large-eddy simulations models perform well (e.g., WES; Santoni et al. 2025, Phys. Rev. Fluids; Chatterjee et al. 2025, PRX Energy).

Lines 19-22 – The statement about gravity waves and blockage is irrelevant for the paper and should be omitted. Although there are several papers claiming that gravity waves excited by a wind farm can result in significant power production reduction recent study by Khan et al. (2025, WES) demonstrates importance of proper design of a numerical simulation including Raileigh damping. Also, LES of a finite size wind farms by Sanchez Gomez et al. (2023, WES) did not exhibit impactful gravity waves. This study was conducted using a compressible numerical weather prediction model in LES mode.

Line 36 – It is not clear what “populations” are referenced to.

Line 74 – Reference Stoll et al. (2020) is not the first one where Boussinesq approximation was used in LES.

Line 84 – Coriolis parameter depends on the speed of Earth’s rotation and latitude, so these should be used as parameters.

Equation 5 – The definition of buoyancy force for incompressible flow where buoyancy effects are introduced through the Boussinesq approximation is not correct. The numerator should be a difference between local temperature and reference temperature, not the horizontal average of the temperature. Simulations with buoyancy defined as in Equation 5 would not result in correct simulations of a stably stratified ABL.

Line 173 – Instead of “fluid” it should be “momentum.”

Equations 30, 31, 32, and 33 – Superscripts (1 and 2) for momentum and heat stability functions can be confused with powers, it would be better to use subscripts instead.

Line 209 – It is not clear why the constants from Beare et al. (2006) and Arya (2001) where combined, why were they not taken from a single study. This should be addressed.

Algorithm 1 – Superscript “t” over q_w should be omitted – this is not a tangential flux.

Algorithm 1 – Turbulent surface stress can also have a second component in cross-flow direction because of wind veering due to the Coriolis force.

Line 224 – It is not clear what is meant by "linkwise-manner."

Line 229 – Instead of “Coriolis parameter” it should be latitude.

Line 242 – Instead of Allaerts and Meyers (2017) more appropriate reference would be Khan et al. (2025, WES) since it is not likely that the damping was applied well in Allaerts and Meyers (2017).

Line 248 – Smagorinsky (1963) is not the correct reference – Smagorinsky used the strain rate magnitude as a numerical viscosity to make sure that his global simulations are stable, the actual development of a subgrid model for parameterization of inertial range turbulence in large-eddy simulations was presented by Lilly (1966, https://opensky.ucar.edu/islandora/object/manuscripts%3A861, and 1967, in the Proceedings of IBM Scientific Computing Symp. Environ. Sci., p. 195).

Line 271 – It is not clear what is meant by "free [atmosphere?] lapse rate." Where is this laps rate applied?

Line 274 – What is “original description” should be made more explicit.

Figure 1 – The difference between left and right panels should be explicitly addressed in the caption.

Lines 277-279 – It is not clear why is the NCAR LES described here – a reference would be sufficient.

Line 283 – Why were C and D grids used? Why not B and C? This must be addressed.

Figure 2 – Same as Figure 1.

Figure 4 – The total stress and TKE should be shown, resolved + subgrid, like in Berg et al. 2020, Fig. 9.

Line 296 – Based on Equation 4 the geostrophic wind and mean wind should be aligned above the boundary layer. Some plausible explanation should be offered for the reason why they are not aligned.

Line 307 – Considering that the model does not include an equation for subgrid kinetic energy the TKE defined here is only the resolved component.

Line 310 – The statement in parentheses is not correct, TKE is shown in Fig. 12 of Berg et al. 2020.

Figure 7 – Again, total turbulent stresses and fluxes should be presented. The ratio of resolved to subgrid stresses depends on gird size and numerical scheme used (i.e., effective resolution) and the SGS model used.

Figure 8 – Why are these not compared to results in Beare (2006)?

Line 390 – These are not "observed" rather "predicted."

Line 399 – It would be good to try to provide more information about computational performance.

Appendix C – The title should be “Convergence study.”

General Comments

This study focuses on the development of a novel method for conducting LES in stratified atmospheric boundary layers in homogeneous conditions. The manuscript is very well written and clearly organized. I appreciate the explanation of methods tested and expansion of the reasoning for selecting one approach over another in Appendix A. The manuscript documents the model methodology and compares the simulation of both a neutral and a stably stratified boundary layer to relevant state-of-the-science models.

I recommend the manuscript be published after minor revision / technical suggestions.

Minor Revisions

• In the explanation of applying MOST at the lower boundary, it was unclear to me how the mean and perturbation are calculated: is the mean calculated in space (entire domain) or in time (local to each cell)? I bring this up due to the application of spatial averaging when dealing with heterogeneous conditions (e.g., including wind turbines in the domain) as opposed to time-averaging which would work for heterogeneous conditions. (Note after finishing the full paper: is this maybe hinted at in line 428?)

Technical Suggestions

• Many of the plots use red and green lines for the simulations of interest which are difficult for colorblind readers. They use different line styles to compensate, but a different combination of colors could be helpful for figures such as Fig. 7.