We thank the reviewer for their excellent comments. Our response, denoted as AUT, is shown below. Modified text of the paper is shown in italic between quotes, while the reviewer’s comments are denoted as REV.

REV: I have read the manuscript with great interest. The manuscript and figures have been 

carefully prepared. The study introduces a new simulation code to the community and new 

simulation strategies that can help study the effect of gravity waves and wind farm blockage. The 

results for the wind farm in section 5 are not directly compared to related literature references. 

Given that the manuscript already discusses so many aspects, that is logical [the corresponding 

section in the introduction needs improvement; see below]. Statements on the benefit of the new 

control methods should be made more carefully [as now these statements are very bold] as 

indicated below. The ideas behind the presented controls are very interesting and the manuscript 

clearly outlines why this needs to be considered. Overall, this is an excellent and very important 

study for the community and certainly qualifies for publication in Wind Energy Science. Below I 

prepared a list of questions/suggestions/recommendations. I hope this helps the authors to improve 

the manuscript. 

REV: comment 1, line 15: "The proposed methodology overcomes these issues by ensuring that the 

same hub wind and thermal stratification are used in successor simulations, regardless of the 

adopted code or precursor run time." --> This should be formulated more carefully; see also below. 

Also, with the proposed methodology, there can be differences between codes.

AUT: rephrased to "The proposed methodology overcomes these issues by ensuring that inflow 

conditions produced from different codes feature the same hub wind and thermal stratification, 

regardless of the adopted precursor run time". Although the shape of the velocity profile, the wall shear stress, the minimum heat flux or the geostrophic wind can differ slightly between codes, hub wind and potential temperature profile are exactly the same, as they are controlled variables. We rephrased the sentence so that it specifies that these variables are the same in the inflow condition, while they can differ slightly in the successor due to the interaction between fringe region and gravity waves, as mentioned later in the paper. We also highlight that such methodology is more focused on the wind farm rather than the CNBL comparability. In fact, ensuring the same hub wind yields very similar wind turbine power, at least at the first row, while the same temperature profile yields the same gravity waves physics, which has been demonstrated to influence the wind farm operation.

REV: comment 2, line 30: references in this section must be checked for accuracy and relevance. 

Some of these studies consider wind farm blockage, not just turbine-induced blockage, as 

suggested. Gribben and Hawkes does not seem to be a peer-reviewed publication.

AUT: we removed Segalini and Dahlberg 2019, as it is an experimental campaign that studies wind farm blockage as a function of the wind farm layout. This study may not be relevant for the present paper. Moreover, we moved those studies that introduce linear superposition of individual effects by rephrasing the sentence slightly. Regarding Gribben and Hawkes, the reviewer is right to say that it is not a peer reviewed article, but the Rankine half body is now an industry standard in firms such as UL Renewables and RWE, so we believe the study is relevant. We hope that the reviewer will agree on this aspect.

REV: comment 3, line 33: no reference is given.

AUT: if the reviewer refers to the superposition of individual effects, now the reference has been added. For the following sentence: “However, recent studies suggest that this could underestimate --- if not totally misrepresent --- wind farm-level blockage, which is heavily influenced by the mutual interaction between the wind farm and the density-stratified atmospheric boundary layer (ABL)” references are given at the end of the sentence.

REV: comment 4, line 35: same as line 30. Make sure the correct references are cited in the correct place.

AUT: the references Smith, 2010; Wu and Porté-Agel, 2017; Allaerts and Meyers, 2017, 2018, 2019; Centurelli et al., 2021 all point to the conclusion outlined above, i.e. that wind farm-level blockage could be underestimated by applying superposition of individual effects alone. All these studies suggest that wind farm – ABL interaction has a leading role in the blockage produced by large wind farms.



REV:comment 5, line 40: Anderson, 2009; Openwind; are not scientific publications. Within the context of this paper, you may not need so many references here.

AUT: these references have been removed.  

REV: comment 6, line 69: the selected references are not the common SOFWA references, in particular Matiz-Chicacausa and Lopez, 2018 and Zhang et al., 2022.

AUT: these two references have been removed and replaced with two more cited references (Churchfield 2012a and 2012b).

REV: comment 7, line 74: better to replace the reference to one of the new grand challenges papers.

AUT: rephrased using “Scientific challenges to characterizing the wind resource in the marine atmospheric boundary layer”, Shaw et al. 2022, Wind Energy Science.    

REV: comment 8, line 82/83:  the results indicate some artificial effects are still observed. The statement "avoids artificial wave reflections" is too strong. Effects in the top of the domain by design are limited, but some effects are still visible in the remainder of the domain. 

AUT:  substituted "avoids" with "limits", which is indeed more correct. Regarding some other 

effects pointed out by the reviewer, we provide some more explanation below:

1. If the reviewer refers to the wave images from the periodic spanwise sides of the domain, we note that these are not reflections but interactions with waves that are physical – in the sense that they are a physical result of the actual simulation setup – as they are produced by the wind farm images mirrored w.r.t these boundaries. To remove these waves completely, a larger domain would be sufficient.
2. Conversely, wave reflections are introduced by the fringe region and by the Rayleigh damping region. These are not physical as energy cannot be radiated back to the wind farm in reality. Such spurious effects will be always present to a certain extent, as the Rayleigh damping and fringe region are proportional controllers, i.e. an error has to exist for the damping to be actually different from zero. In practice, there is no way to avoid reflections and some wave energy will always be reflected inside the domain. Moreover, the fringe region itself triggers gravity waves as it displaces the BL back to its original state, and also this effect is unavoidable. To make sure that these do not interact with the physical ones we use the advection damping region, which we found in later studies to be crucial, especially when a strong capping layer is present. By killing the horizontal advection of vertical velocity, these spurious perturbations are not transported, thus they cannot interact.               

REV: comment 9, line 98: does this mean TOSCA only works with "openMPI" or would it also work with other MPI libraries? I assume the latter.

AUT: yes, although TOSCA has only been tested with the OpenMPI implementation of the MPI standard.    

REV: comment 10, line 107: the references "(Loken et al., 2010; Ponce et al., 2019)" are not about the Niagara system and should be removed.

AUT: we refer to the following web page https://docs.scinet.utoronto.ca/index.php/Acknowledging_SciNet, where these are stated as the publication to be used to acknowledge Niagara and the SciNet datacenter.    

REV: comment 11, in section 2.4: the manuscript compares with Martínez Tossas et al. (2012). 

However, an updated version of this work is provided at https://doi.org/10.1002/we.1747. Please 

update the reference accordingly. Also, please check whether the results have been updated 

because some numbers have been changed compared to the 2012 version.

AUT: we thank the reviewer for having brought this aspect to our attention. Velocity plots seem exactly the same, while power vs grid spacing plots have been changed. In particular, turbine powers from OpenFOAM are higher now and they compare better with TOSCA. Also, the reference has been changed to Martinez Tossas 2015.

REV: comment 12, line 378:  the manuscript compares with Martínez Tossas et al. (2012). However, an updated version of this work is provided at https://doi.org/10.1002/we.1747. Please update the reference accordingly. Also, please check whether the results have been updated because some numbers have been changed compared to the 2012 version.

AUT: the sentence has been rephrased according to the reviewer comment to "This could be beneficial for example when making comparisons between different codes, [...]".

REV: comment 13, figure 1: concerning the introduced inflow generation method. The only aspect discussed is the spanwise periodicity, for which it is claimed that it can be removed completely by the Munters 2016 method. However, if the domain is as small as indicated in Figure 1 there could be additional effects as the large-scale flow features in the boundary layer are not fully resolved (see, for example, the work by Javier Jiménez and coworkers). So there would be limitations in size.

AUT: yes, there are limitations in size for CNBLs, but our argument is that, for large wind farms, the limitations on the domain imposed by the gravity waves and by the farm itself are much bigger than the ones that are required to resolve a CNBL. When using the concurrent precursor technique, precursor and successor should have the same spanwise and vertical size, so performing the spin-up on a precursor as big as the successor would be a waste of resources if a smaller domain – which still resolves the CNBL – could be used. Our method decouples the two domains, allowing to perform the spin-up with the domain size that is suitable to resolve the CNBL, and the successor phase on a different domain that suits the requirements of the wind farm – atmospheric interactions.

REV: comment 14, line 441-444: : the formulation, as is, suggests atmospheric gravity waves are not real flow physics, which is not what the authors want to convey.

AUT: the spin-up phase and, in general, the precursor phase, do not contain gravity waves as there is no perturbation to the CNBL (there are small gravity waves, generated by the turbulent structures that hit the inversion layer, but they are negligible). For this reason the requirement on the domain size is much more computationally feasible at this stage. In any case, the paragraph has been rephrased to make this concept more clear, namely "Nevertheless, the proposed hybrid method [...] is very convenient as it allows to reduce the overall computational cost of the ABL spin-up phase,where wind farm-induced gravity waves are not yet present. In fact, this initial phase is run on a domain whose size is dictated by the current flow physics, rather than on quantities that will only become relevant at later simulation stages".         

REV: comment 15, section 3.1: compare to https://doi.org/10.1002/we.1747 instead of the 2012 version.

AUT: the reference and data have been modified according to the new version.     

REV: comment 16, line 462: what is the azimuthal direction here?

AUT: this is the discretization along the direction that is tangential to the disk perimeter. This is an angular discretization, as its actual size changes with radius, that is why sometimes is referred to as the azimuthal discretization.

REV: comment 17, Table 1: double-check the numbers in the table.

AUT: updated to Martinez-Tossas 2015 as suggested by the reviewer.       

REV: comment 18, line 482: "while TOSCA’s ALM predicts a slightly higher power". --> The difference is nearly 15%, certainly more than "slightly higher". One would expect a closer agreement since this concerns a uniform inflow case. Perhaps also compare to more recent simulation studies of that case performed with different codes, see e.g. https://doi.org/10.1063/1.5004710, https://doi.org/10.1002/we.2714.

AUT: both data from TOSCA and OpenFOAM have been changed. The latter has been replaced with data from Martinez-Tossas 2015 paper, which show a higher power for the ALM. TOSCA data have been initially calculated without excluding the initial transitory by mistake. Now data have been averaged so that the initial transient of the flow and turbine controller is avoided, as stated in the paper. Data are indeed closer (6% difference w.r.t. the average between the four powers for the ALM).



REV: comment 19, line 519: the reference is incorrect, as Chamorro and Porté-Agel (2010) only consider the single turbine case. 

AUT: reference has been changed to Chamorro and Porté-Agel (2011) throughout the paper, this was a typo. 



REV: comment 20, line 586: what exactly is the model parameter c (perhaps I missed the definition)?

AUT: this parameter is a very important one within the Rampanelli and Zardi model (defined after Eq. 13 in their paper), as it essentially controls how sharp the capping inversion layer is. As one can imagine, very different results can be obtained by solely changing c. For this reason, we believe it is important to specify the choice of such parameter when detailing the initial condition of the simulation. In our study, we set it to 0.33 as Rampanelli and Zardi suggest 1/3 as a suitable value for CNBLs.

REV: comment 21, figure 6d: is the S2 case supposed to look like the case PT? This is unclear from the text. The temperature profile for case PT looks like it is intended.

AUT: we rephrased the paragraph to make it a bit clearer: “This results in a higher inversion height for a given set of ABL parameters, and can be observed by comparing cases G and S2, which feature the same wind-angle controller, but which differ in the obtained profile of potential temperature. This leads to an increased surface temperature predicted by TOSCA and a different wind veer profile between the two codes. Although we note that such differences are accentuated by the fact that statistics from SP-Wind are collected at an earlier time, i.e. the CNBL has grown by a lower extent, case S2 seems to be more aligned to case PT, where the average potential temperature profile is kept constant by the controller. The difference in mixing between the two codes also affects the average hub-height velocity, which differs by 0.33 m/s between case G and S2".

This means that case S2 should be exactly the same as case G if SP-Wind and TOSCA would be the same, as these two cases feature the same controller. Conversely, case P and PT feature the pressure controller with geostrophic damping, which is slightly different. One aspect of notice, also mentioned in the text, is that statistics from S2 have been collected at an earlier time, hence the CNBL had less time to mix (lower ground temperature, lower inversion height etc).  The fact that case S2 and PT agree simply means that SP-Wind predicts less mixing than TOSCA, as in SP-Wind the potential temperature profile only moves slightly w.r.t to case PT, where it is forced to be constant in the spatial average, equal to its initial condition. This even more highlights the need to use a temperature controller when running precursor simulations as, by the time that turbulence has developed and the CNBL has formed, the potential temperature profile has shifted from its initial shape.

To sum up, approaches G and P are almost exactly equivalent if applied to the same code base, while the use of the temperature controller allows to compare different code bases. If one doesn’t want to control temperature but desires to quantify the amount of mixing between different codes, the comparison has to be made by starting and ending the averaging at the same time. If this doesn’t happen, potential temperature evolution could be confused with more or less mixing in the code.

We hope this clarifies the reviewer comment/doubt.

REV: comment 22, figure 6: one can see that the different controls have some effect on the flow. 

AUT: the effect is mainly due to the temperature controller, which affects the heating history (cases P or G vs case PT). Conversely, the different velocity controllers (P or G) can be used interchangeably. Controller G is preferred when the geostrophic wind is to be controlled, while controller P + geostrophic damping is better suited to control the wind inside the boundary layer.

REV: comment 23, around line 620 : keeping the temperature profile fixed is also an idealization. I agree there are reasons to do so, but keeping the temperature profile fixed also changes boundary layer dynamics as the way buoyancy is represented in simulation changes a bit



AUT: we do not believe that fixing the temperature profile is more neither less physical than not doing so. It is only more reproducible as it rules out the dependency of its evolution on the adopted code, numerical discretization and turbulence closure both in the momentum as well as potential temperature equations. This applies to Boussinesq based solvers; we didn’t investigate the effect of fixing potential temperature e.g. to anelastic formulations.

REV: comment 24, line 623: "Results from different codes would thereby be more aligned with each other" ==> This is not shown and may depend on the quantity of interest.

AUT: rephrased to “This would lead to a better agreement on wind farm power predictions and on 

the actual inflow conditions used in successor simulations, allowing to compare precursor 

simulations that did not run for the same amount of time.”

REF: comment 25, figure 7: one can see some oscillation also remains towards the end of the simulation, which is why statements that all oscillations are completely removed should be adjusted.

AUT: referring to Fig. 7, if the reviewer refers to oscillations in the G case, this does not feature geostrophic damping for the reason mentioned in the paper (i.e. the geostrophic wind is known a-priori and it is possible to provide a consistent initial condition), hence oscillations are present, especially within the boundary layer where the initial condition cannot be consistent (these are small as they are limited by turbulence which acts as a damping term). Conversely, case P does not show any noticeable oscillations after the damping time mentioned in the paper.

REV: comment 26, line 627: "differs slightly" --> please define.

AUT: Rephrased to “differs by 0.25%”. See Tab. 5, geostrophic wind (G) incases G and P is 12.00 m/s (prescribed) and 11.97 (calculated) m/s respectively. 

REV: comment 27, line 656: "Although these streaks do not alter the simulation results" --> probably depends on what quantities are considered. The periodicity of the flow in a spanwise direction is visible in Figure 8

AUT: rephrased to “If averages are gathered for a sufficient amount of time, these streaks do not alter the simulation results from the wind farm performance point of view. Nevertheless, as statistics convergence can become extremely slow, this issue can be alleviated by using the so-called shifted periodic boundary condition [...]”. We believe that the periodicity in the streaks does not alter any average results if averages are taken for a sufficient amount of time. We also carried out our simulations further forward in time to ensure such statement. Another way to remove the streaks is to increase the spanwise size of the fringe region, so that they do not remain locked in position between the inlet and outlet of the concurrent precursor. As mentioned in the paper, we also think that such phenomenon is enhanced by the choice of equivalent roughness height, which is higher than other studies.

REV: comment 28, figure 9: can you comment on the origin of the oscillations shown in Figure 9b?

AUT: these gravity wave patterns are free of spurious (i.e. non-physical) waves. The waves that seem reflected are instead waves interactions produced by the fact that the side boundaries are periodic. Hence, to the eyes of our LES set-up, it is as if many domains were staked laterally, and wave trains that are exiting from one side are re-entering from the other side. However, it is important to highlight that these are wave transmissions, not reflections, across the periodic boundary, not a numerical artifact such as gravity waves reflecting due to a boundary condition (which would happen for example with a zero-gradient BC). Moreover, there isn’t any attempt in current state of the art wind farm LES in avoiding such lateral interactions, as they do not affect wind farm operation if they happen far and downstream with respect to the wind farm. That is why we chose a domain which is fairly large in the spanwise direction, allowing waves to be close to the outlet when they reach the domain sides. By doing so, their images will re-enter the domain already close to the outlet and far downstream the wind farm   

REV: comment 29, figure 10: what do the symbols in Figure 10b represent?

AUT: they are sample probes, for this reason they are finer where we expect the blockage to have a high gradient, transitioning from global to local induction. We removed the symbols as they are misleading.

REV: comment 30, line 676: as different cases are considered, these numbers cannot be compared one-to-one. 

AUT: We removed the reference to Wu and Porté-Agel, which were provided for reference. Another aspect to notice is that their results might be affected by high wave reflectivity, as the vertical domain size, set to 2.4 km in their simulations, was insufficient to vertically resolve the gravity waves, characterized by an estimated vertical wavelength of 10 km and 4.8 km for the weak and strong lapse rates respectively.      

REV: comment 31, line 679:  "Overall, these wave patterns agree well with linear theory (Allaerts and Meyers, 2019)" --> This is not shown in the manuscript. This statement should either be removed or made more specific. 

AUT: this was a qualitative comparison that can be assessed e.g. from Allaerts & Meyers 2019, but we removed the whole sentence for simplicity, as the case shown in Allaerts & Meyers 2019 is run with different CNBL parameters and could be confusing.     

REV: comment 32, figure A1: the low-level jet looks quite sharp. Is there any reason for that?

AUT: it may be a consequence of the chosen CNBL input parameters. In this case, the vertical resolution has been refined to 4.6 m within the inversion layer to capture the Ellison scale.



REV: comment 33, line 757: "Brunt-Vaisala frequency" add accents.

AUT: accents added.

We would like to thank the reviewer for the time dedicated to revising the paper. We proceed with answering and clarifying, where possible, the proposed comments.

Our response, denoted as AUT, is shown below. Modified text of the paper is shown between quotes, while the reviewer’s comments are denoted as REV.

REV: The manuscript provides a detailed overview of a new LES model for wind energy applications based on finite volume discretization. The manuscript also includes presentation of model tests for simulation of conventionally neutral boundary layer (CNBL) as well as simulations of an isolated wind turbine rotor and a wind plant. The simulation of the CNBL is compared with previous studies that utilized three different LES models. The model performance in reproducing CNBL is comparable to that reported in previous studies. The simulation of an isolated wind turbine rotor is comparable to that reported by Stevens et al. (2018).

AUT: We would like to clarify that the isolated wind turbine is not compared against Stevens et al. (2018), but rather against Martinez Tossas 2015 (it was compared against the 2012 preprint in the manuscript, but that has been already changed and the data updated, according to comments from the first reviewer). Furthermore, there are two more studies in the paper, not mentioned by the reviewer, which are in our opinion relevant. First, we compare an infinitely wide wind farm without thermal stratification against numerical results from Stevens et al. (2018) and experimental results from Chamorro and Porté-Agel (2011), showing that TOSCA agrees well with both previous numerical models and against wind tunnel measurements. Secondly, we simulate a finite wind farm under CNBL (neutral stratification within the BL and stable aloft, with an inversion layer), quantifying blockage for this specific scenario, providing valuable details on the setup of this kind of simulations, that are useful in order to avoid contaminating the results by unphysical wave reflections. In particular, we provide considerations on the domain size, damping regions and propose a novel precursor methodology that allows to resolve gravity waves at a lower computational cost. We also note that—to the best of our knowledge—prior to our study, resolving gravity waves in an LES of a wind farm has only been achieved by the group of Meyers, which uses the non-open source SP-Wind code. For this reason, we believe that showcasing this capability of TOSCA while at the same time providing the source code used to achieve it, are two aspects that strengthen the novelty and importance of the publication.

REV: The authors present several flow field controllers that enable specifying stationary flow conditions. They introduced a novel hybrid off-line precursor/concurrent-precursor method for generation of turbulent inflow for non-periodic streamwise flow conditions. In addition to the LES model itself this is the novel component of the current study.

AUT: We are glad that the reviewer found this aspect of the paper useful. Nonetheless, we believe that in addition to the methodological advances, our LES results showing gravity waves triggered by a finite wind farm are also novel both due to the scarcity of similar studies in literature and because the code used to produce them is made available to other researchers.

REV: The presentation of the new open-source LES model and the off-line precursor/concurrent-precursor methods are focused on model development and that brings two questions:

First, the focus on the manuscript is on model development and not a well formulated scientific question, suggesting that this manuscript is better suited for a model-focused journal, like Geoscientific Model Development, rather than Wind Energy Science which expects exploration of a scientific question.

Second, the motivation for the development of a new model is to provide an open-source platform, and that is a noble goal. However, considering that there are numerous other open-source models for wind energy applications like SOWFA, PALM, and Exawind/AMRwind (as mentioned in the manuscript), it is not clear what is clear distinguishing characteristic of the TOSCA model other than it is a finite-volume model. It is stated that SOWFA is not “sufficiently parallel-efficient,” but no evidence is given for that statement. No comparison in performance between TOSCA and other models is provided.

Third, the authors ignore a significant body of work focused on creating realistic simulations of wind turbine and wind plant flows including comparisons with observations (e.g., Mirocha et al. 2014, Aitken et al. 2014, Marjanovic et al. 2017, Arthur et al. 2020, Sanchez Gomez et al. 2023). Rather than include comparisons with observations, the authors instead opt for creating different flow controls to achieve stationary conditions, without inertial oscillations or wind veer, essentially creating idealized conditions that are not commonly observed in the atmosphere. Better justification of new idealized scenarios would be required.

Finally, there are a number of imprecise or inaccurate statements made that are listed below, under Specific Remarks.

Taking all the above into account, Wind Energy Science journal may not be appropriate journal for publication of the work primarily focused on model development and idealized simulations. Perhaps Geoscientific Model Development journal would be more appropriate, but only if the numerous issues regarding inaccurate statements about the state of the science addressed and modeling choices better justified.

AUT: The reviewer’s opinion that the paper is primarily focused on model development is fair. Nonetheless, we note that Wind Energy Science has published numerous studies that primarily report on model development, such as Martinez-Tossas et al. 2021; Bradstock and Schlez 2020; Blondel, F. and Cathelain 2020; to name a few.

The research question that motivated the development of TOSCA is detailed in the introduction. Specifically, as wind farms are getting larger and closer to each other, there is a need of better understanding blockage effects and cluster wakes. This is a complicated task which cannot be accomplished using only observation data. As a consequence, the need for high-fidelity models is always present, as these are run in a more controlled environment and data are more easily available at the desired locations and in the desired time window. However, the available high-fidelity models are not uniformly capable of capturing the physical phenomena involved in wind farm blockage, slowing scientific advance. We argue that our paper adds to the scientific inquiry of wind farm blockage by (i) contributing an open-source tool that is specifically tailored for investigating these phenomena, and (ii) suggesting a consistent method for initializing wind farm simulations.

Regarding other open-source codes, we agree with the reviewer that at least three more platforms are available. Regarding SOWFA, we mention that this is not well-suited for large wind farm simulations as it is based on OpenFOAM, a general-purpose code with limited performance at scale. This is due to its intense I/O, which results in the generation of a massive number of files when running on thousands of cores (one folder per processor with several field files). We are also aware that OpenFOAM features a “collated” I/O option, but this definitely is too slow for large size problems, as it features multiple MPI_AllReduce calls for every field writing. We do not want to include a scaling performance comparison in our paper, as a fair comparison would require the same HPC platform and case. Besides, the NREL itself acknowledged such limitations as it stopped developing SOWFA and it is now concentrating on the ExaWind project. While ExaWind may become a leading tool in the future, its development is still at an initial phase. It does not feature an IBM method and, as well as PALM, lacks a concurrent precursor method, which we found to be essential if one wants to at the same time prescribe a realistic turbulent inflow and damp gravity wave reflections, as mentioned in the paper. TOSCA is at the moment the only finite-volume open-source code that allows both inlet/outlet (off-line) precursor field mapping and concurrent-precursor technologies. Making it to our knowledge the preferred open-source tool when the objective is to capture gravity waves around a large wind farm.

Nevertheless, to make this concept clearer, we rephrased to “In addition, although SOWFA has been used in several research studies in the last decade […], it is not sufficiently parallel-efficient when running with thousands of cores, as the number of produced files increases drastically with processor count. While some of these shortcomings have been addressed and solved in the NREL Exawind project […] the latter is not yet at a production stage. Moreover, all SOWFA, PALM and Exawind platforms do not feature the concurrent precursor technique, making TOSCA the only finite-volume open-source code to possess such capability.”

Regarding the third comment of the referee, we would like to point out that our CNBL simulations do feature wind veer, as the Coriolis force is considered in the governing equations (not in the infinite wind farm and isolated wind turbine case, but those cases are matching other validation data which didn’t model wind veer). Secondly, the velocity and temperature controllers are models to represent existing phenomena such as large-scale horizontal pressure gradients and subsidence respectively. Our choice to simulate a non-evolving flow (mainly in terms of BL height) is explained in the paper, and we believe that this is not more idealized than running a precursor for big amount of time without temperature control, as in reality CNBLs are only observed for short lapses of time, especially in transition between night and day. Without temperature controller, by the time turbulence has developed and reached a statistical steady state, the BL height has increased and the inlet conditions for the wind farm in terms of potential temperature structure are not the ones initially prescribed in the precursor, and should be recomputed.

We agree with the reviewer that the most physically-realistic simulation platform would be the one that is coupled with mesoscale models (such as WRF), but this is out of scope for the two objectives that have driven the initial development of TOSCA, i.e. blockage and cluster wakes. Nevertheless, it is not excluded that we will add such feature in the future, when interested in simulating wind farms in real-life conditions also from the mesoscale point of view. On the same note, we would also like to highlight that the numerical implementation that is necessary to drive LES simulation with realistic mesoscale information is essentially the same that we propose, with the only difference that the average reference state changes in time (Allaerts at al. 2020, Allaerts at al. 2023). Despite this, the message that we want to convey is that it may be useful, when studying the ABL response to the wind farm in terms of gravity waves and wakes, to impose a certain wind speed and a very precise average thermal stratification (imposing a given roughness height is trivial), as these are the parameters upon which such response mainly depends. Having certain parameters such as the hub-height wind speed, BL height and inversion strength varying between codes, given the same initial condition, introduces uncertainties which we argue could be eliminated by adopting the proposed controlling methodology. Finally, within the paper, we mainly talk about realistic inflow conditions, or realistic turbulent inflow, which refer to prescribing an inflow that is realistically turbulent as it results from a prior simulation.

Regarding the reviewer’s conclusive comment, we believe that our CNBL simulations are as idealized as the majority of CNBL precursors available in literature. They feature a realistic velocity and temperature profile, as shown in the paper, with the addition that the proposed controlling techniques allow to have direct control on the variables of interest. We also believe that trying to study blockage effects or wake effects in varying mesoscale pressure gradients, BL height or thermal stratification would provide an inconclusive picture of the underling physics, as too many variables would be changing at the same time.

REV: Abstract, line 9 – It is stated that “…prescribe a desired average hub-height wind speed while avoiding inertial oscillations above the atmospheric boundary layer (ABL) caused by the Coriolis force, a known problem in wind farm LES studies.” However, inertial oscillations are not necessarily a problem: they are a component of a realistic atmospheric boundary layer flow in near-neutral and stably stratified conditions (e.g., Baas et al. 2011, QJRMS). Spurious inertial oscillations that are part of model-spin-up are a problem, but realistic inertial oscillations should be included in simulations that seek to match the real world.

AUT: We agree with the reviewer that inertial oscillations of the atmospheric boundary layer are physically realistic and important. Nevertheless, when one desires to impose a given hub-height wind speed instead of a geostrophic wind, the inertial oscillations that are triggered are spurious, as they arise because of a mismatch between the initial condition and the equilibrium solution (which is initially unknown as explained in the paper). As a consequence, these oscillations associated with the initial spin-up of the model should be removed, as they would not be part of the solution if, for example, one forced the BL using what is referred to as the geostrophic controller in the paper. Moreover, again in the context of gravity waves, their vertical wavelength depends on the geostrophic wind magnitude. As a consequence, when trying to characterize their impact it is evident that having a non-stationary free atmosphere wind (an effect that is in the equations but it is spurious in this case) would probably lead to wrong conclusions if for example averages are not performed within one oscillation period (which is rather long and unfeasible computationally for a whole wind farm simulation).

REV: Abstract, line 14 – It is stated that “The proposed methodology overcomes these issues by ensuring that the same hub wind and thermal stratification are used in successor simulations, regardless of the adopted code or precursor run time.” Based on this statement it would seem that the goal of this work is to make all codes produce the same results regardless of the numerics. Observations do not seem to be relevant. It is not clear that this is a scientifically valid goal.

AUT: Our goal is to propose a method that would facilitate comparisons across codes by ensuring that wind farms are simulated under consistent initial conditions in terms of thermal stratification, wind speed, and direction. We believe that this will help systematize computational investigations of large-scale wind farms. Moreover, our statement does not convey that such approach should be always used (for example it is not of interest when trying to couple LES and mesoscale models) but we deem it useful when trying to simulate the wind farm behavior under a specified thermal stratification and wind speed and direction, which many researchers do.

REV: Abstract, line 17 – It is stated that “periodic wind farm are also presented to assess TOSCA’s ability to model large-scale wind farm flows accurately and with high parallel efficiency.” Here, it is not clear how is accuracy or parallel efficiency assessed – no comparison to observations or assessment of computational performance is presented.

AUT: We asses the accuracy of the model by comparing our results against existing and well validated numerical codes. Moreover, the periodic wind farm validation case is also compared against wind tunnel measurements. Finally, in the appendix we included TOSCA’s strong scaling test results that indeed prove that the code can efficiently run at scale (576M cells with 32.000 cores).

REV: Line 27 – Lundquist et al. 2019 (Nature Energy) presented simulations of wind plant losses due to wind plant wakes.

AUT: Citation added, thanks.



REV: Line 35 – It should be pointed out that Smith 2010 presented a 2D linearized analysis and that Wu and Porte-Agel 2017, Allaerts and Meyers 2017, 2018, 2019 used incompressible codes with Boussinesq approximation, and Allaerts and Meyers use 25 km vertical domain size, which is not appropriate for incompressible simulations (c.f., Lilly 1996 Atmos. Res., Wood and Bushby 2016, JFM). By citing papers based only on incompressible simulations with an inappropriate domain, the authors present an inaccurate picture of the state of the science.

AUT: The references mentioned by the reviewer and provided in the paper are indeed pointing to the very conclusion that we try to convey, i.e. that thermal stratification has an important effect on wind farm blockage. From the date the manuscript was submitted, more papers appeared in this regard: Smith 2023: The wind farm pressure field, and Lanzilao and Meyers 2023: A parametric large-eddy simulation study of wind-farm blockage and gravity waves in conventionally neutral boundary layers.

Regarding the results from Smith, we highlight that they employ linear gravity wave theory, which also features the Boussinesq approximation. Several books (Nappo 2012, Lin 2007) or articles (Allaerts and Meyers 2019, Smith 2002, Smith 2006, Smith 2010) use linear theory to model free atmosphere gravity waves, where the latter are produced by vertical disturbances in the flow generated by terrain features or wind farms.

In in the context of the present work, one should pay attention that the vertical length scale of interest for the validity of the Boussinesq approximation is not the vertical domain size, but rather the scale height of the atmosphere, typically around 10 km. Regarding the references pointed out by the reviewer, Lilly 1996 demonstrates that the Boussinesq approximation does not strictly conserve mass, a general property of such approximation. Conversely, Wood and Bushby 2016 perform an analysis that addresses the predictive difference in the onset of convective instability between Boussinesq and compressible models. As such instabilities eventually produce a turbulent convective flow in the entire domain (they also cite the Rayleigh-Benard convection) in their case the whole domain length is of interest. In our study, we model a simple CNBL without moist, deep or shallow convection or cloud formation and precipitation. Here, stability is neutral up to the inversion height, where it is capped by a stable inversion layer. This prevents the BL from reaching the top of the domain, and it is followed by a stable lapse rate aloft, where the flow is laminar, referred to as the free atmosphere. Atmospheric gravity waves would develop in the free atmosphere, inducing disturbances in the flow variables that can be studied using the Boussinesq approximation, as done in almost every LES study published so far with e.g. the ETH code, SP-Wind, SOWFA, AMRWind, PALM and TOSCA. In particular, in order for the latter to hold, the following hypotheses should be verified (Lin 2007): (1) the vertical dimension of the fluid motion is much less than any scale height and (2) the motion-induced fluctuations in density and pressure do not exceed, in order of magnitude, the total static variations of these quantities. These hypotheses are verified, as wind farm gravity waves only slightly perturb the pre-existing equilibrium, inducing extremely small motions. These waves can become unstable and break into clear air turbulence, but this is likely to happen only with variable lapse rates, when the local Richardson number becomes less than its critical value (Klemp and Lilly 1978). We use a constant linear lapse rate aloft.

Conversely, we believe that the vertical size of the domain is instead dictated by the vertical wavelength of gravity waves. In fact, it has been shown by several studies (Klemp and Lilly 1978, Allaerts and Meyers 2017, 2018 and Lanzilao and Meyers 2022) that this should be split between a physical domain region and a damping region, which damps out gravity waves before they reach the top boundary and reflect, contaminating the solution. This is achieved by a damping region that is greater than the dominant wavelength, calculated as 2piG/N, where G is the geostrophic wind and N the Brunt-Vaisala frequency (again assuming linear lapse rate). As this is usually of the order of kilometers (around 12 in our study), and it increases with decreasing lapse rate, it appears clear how vertical domain sizes of 25/30 km are easily achieved in literature (at least two wavelengths, one for the damping region and one for the physical region). If this is not satisfied, results may be contaminated by gravity waves reflections produced at the top boundary, or simply unable to resolve gravity waves due to the insufficient space. For example, in Wu and Porte-Agel 2017, the vertical domain size is set to 2.4 km. This is insufficient to vertically resolve the gravity waves, characterized by estimated vertical wavelengths of 10 km and 4.8 km for the weak and strong lapse rate that they analyzed respectively.

REV: Line 39 –Again, this literature review is inaccurate. The statement “wind energy predictions are made using low-cost but fast, often analytical, reduced-order models (Katic et al., 1987; Ainslie, 1988; Larsen, 1988; Anderson, 2009; Openwind; Bastankhah and Port.-Agel, 2014; Niayifar and Port.-Agel, 2016),” is incorrect. The Ainslie model predicts wind flow in wake of a turbine, Katic is for cluster efficiency, and Larsen is wake model, Anderson model is wake model, while Openwind is a commercial tool for wind plant design (and its name is deceptive, because it is not open), so these papers do not refer to wind energy prediction (i.e. wind energy resource assessment of forecasting). The sentence should be rewritten and properly focus either on wind turbine wake prediction or different references should be provided.

AUT: Rephrased to “In industry, annual energy captures are made using low-cost but fast, often analytical, reduced-order wake models” and updated references (removed Openwind and changed 1987 Katic and Jensen to 1983 Jensen, which more clearly refers to the top-hat wake model).

REV: Line 41 – The statement “While they [the models listed above] have been used effectively for hundreds of wind energy projects, the majority of these models currently struggle in accurately reproducing wind farm blockage and farm-farm wake interactions” is not supported by any references. Are these models really used for hundreds of wind energy projects? Perhaps Openwind, but others are research and not design tools. Vague and/or inaccurate comments should be avoided.

AUT: Wake models such as the Jensen or Bastankhah wake models are indeed industry standard in large companies such as RWE or Orsted. Other models are also very popular such as the EV model by Ainslie, later reformulated analytically by Larsen with some approximations (preferred in Openwind, which lately started to propose TurbOPark as its core model). We do not believe that such statement is inaccurate information. Nevertheless, we added the reference from Nygaard 2022, which clearly shows that the park model (i.e. the Jensen model) underestimates wake effects and additional physical phenomena have to be accounted for, such as a variable wake expansion coefficient to better capture the wind farm wake evolution. Regarding blockage effects, it is easy to show that those models completely ignore it, as they do not predict any velocity deficit upstream the wind turbine. Moreover, the work by Centurelli et al (2021) clearly shows that current analytical blockage parametrizations do not agree with LES, but instead underpredict blockage in some cases.

REV: Line 66 – The statement “do not implement the concurrent-precursor method, making it difficult to properly capture gravity waves effects and inlet/outlet reflections” is confusing, it is not clear what is meant with “properly capture… inlet/outlet reflections.” Such reflections should be avoided not captured. Furthermore, such problems are not encountered when using compressible codes like the Weather Research and Forecasting (WRF) model (c.f. Mirocha et al. 2014 JRSE, Aitken et al. 2014 JRSE, Vanderwende et al. 2016 JAMES, Marjanovic et al. 2017 JRSE, Arthur et al. 2020 , Dubrowa et al. 2020 Wind En., Kale et al. 2022 Renew. En., Sanchez Gomez et al. 2023 WES, Bui et al. 2023 EGUsphere).

AUT:  That was a typo, thanks for pointing it out. Rephrased to “making it difficult to properly capture gravity waves effects at the same time avoiding inlet/outlet reflections”. Regarding the references pointed out by the reviewer, it is not clear whether gravity waves could even be captured in those studies. Mirocha et al 2014 and Aitken et al. 2014 could not observe any gravity waves as the physical domain height (1 km) was too small to actually resolve them, so they would be damped immediately by the Rayleigh damping layer (0.7 km). In the paper of Vanderwende et al. 2016 it is difficult to say if gravity waves are present in the domain, as no large-scale pictures or data is shown. Marjanovic et al. 2017 is not relevant for the present study, as they simulated a precipitating event over the Balkans and their domain was extremely large (whole Europe and part of Asia). Arthur et al. 2020 and Bui et al. 2023 might have picked conditions where gravity waves were not strongly excited, upwind propagating or simply died out before reaching the domain boundaries, it is difficult to say as the authors do not mention anything about it. We prefer to avoid speculating in this sense. Kale et al. 2022 simulated isolated wind turbines, so the wind farm – induced boundary layer displacement, which is the triggering mechanism for gravity waves, is minimal and gravity waves are probably absent. In Sanchez Gomez et al. 2023, as for Wu and Porte-Agel 2017, the vertical domain size is not sufficient to resolve gravity waves. Moreover, the mentioned Froude number is a meaningful parameter for surface waves only, which only propagate within the inversion layer. Internal waves always propagate upstream and upwards, but in some cases they can reach the top Rayleigh damping before the inlet, so that they are damped before reflecting (Gabriele Centurelli on PALM, personal communication). Given the very small vertical domain size, this seems to be the case for Sanchez Gomez et al. 2023.

As a general remark, it is difficult to talk about the behavior of gravity waves in WRF without exact knowledge of how the domain nesting technique is implemented, on how the data is transferred from one domain to the next, or given the extremely high variability of the flow properties in such simulations featuring dynamic evolution of the background flow. We believe that WRF would face the same issues in a simulation where the background stratification is steady and the simulation is carried out for a sufficient amount of time. In fact, as pointed out by Klemp and Lilly 1978, wave reflectivity increases over time for a given steady background state. It may be that energy accumulation due to wave reflection might represent a smaller issue for those WRF simulations where the background flow is evolving. These issues are worth investigating, but we feel fall outside the scope of the present paper.

REV: Line 78 – It is stated that “TOSCA is designed to enable LES simulations of large finite wind farms, operating in realistic atmospheric conditions,” however, for simulations of wind plants in realistic atmospheric conditions it is necessary to couple mesoscale and microscale simulations (c.f. Haupt et al. 2023, BAMS). This requirement should be acknowledged.

AUT: Rephrased to “TOSCA is designed to enable LES simulations of large finite wind farms under realistic turbulence inflow and thermal stratification”. In the reminder of the paper, we always refer to realistic turbulent inflow conditions. See above for our response regarding the mesoscale to microscale coupling.

REV: Line 79 – It is not clear what is meant by “mid-fidelity actuator model.”

AUT: Rephrased to “Wind turbines can be modeled using the actuator line […] and the actuator disk […] models”.

REV: Line 86 - In simulations with the compressible WRF model by Sanchez Gomes et al. 2023, WES, there are not gravity waves.

AUT: See above. In particular, too small domain combined with Rayleigh damping layer could completely eliminate their presence.

REV: Line 88 – Horizontal grid stretching can be detrimental to resolving well boundary layer turbulence.

AUT: We agree with the reviewer regarding this aspect. The favorable agreement with the presented validation studies suggests that the current level of grid stretching is not too detrimental.

REV: Line 92 – Previously it was mentioned that CNBL will be simulated and not stably stratified ABL.

AUT: The CNBL is an idealized type of boundary layer characterized by neutral stratification below a strong inversion layer and a stable lapse rate aloft. This can be imposed as an initial condition (Sec. 4 and 5) or reached starting from an initial condition with fully stable stratification and uniform non turbulent wind. As turbulence develops, the temperature profile will depict a region where the mixing generates a neutral stratification, an inversion layer that caps such region and the initially prescribed stable lapse rate aloft (Sec. 3.2). At line 92 we do not mention that the ABL is stably stratified (which would require a heat flux), but rather that it is thermally stratified, conveying that a dynamic evolution is present. This is explained in detail in Sec 3.2.

REV: Line 98 – The sentence starting with “TOSCA is written…” reiterates what was already mentioned in the Introduction in line 75.

AUT: Eliminated.

REV: Line 146 – 155 - This is standard textbook curvilinear coordinate transformation that can be put in an appendix or a reference could be provided.

AUT: Moved to appendix.

REV: Line 230 – Here, ‘k’ should be under the filter overbar.

AUT: Corrected for better readability also in the following lines.

REV: Line 231 - Filtered value of a product is not equal to product of filtered values regardless of discretization, but only approximately equal.

AUT: If the mesh is uniform, the filtering operation on the face area vectors has no effect. Hence the equality holds exactly. If it is stretched, this is an approximation as stated.

REV: Equation 19 –line 237 – This is a standard dynamic model derivation and as such can be put in an appendix or a reference can be provided.

AUT: Moved to appendix.

REV: Line 264 – It is not clear what is meant by “uniform actuator disk (UAD) models.” This is not  standard terminology.

AUT: Rephrased to: “To represent wind turbines, different models have been implemented, namely the actuator line (AL), actuator disk (AD), and uniform actuator disk (UAD) models. Following Tossas Leonardi (2015), the first two models require detailed blade information (i.e. airfoils, twist and chord), while the UAD only requires turbine thrust coefficient and general rotor information such as diameter and hub height.



REV: Line 265 – Reference to actuator disk and actuator line models are provided, while they share elements it would be better to describe and reference them separately.

AUT: Referring to our actuator disk and actuator line models, their implementation is identical. The only difference consists in the fact that points physically rotate in the AL model, while they are stationary in the AD model. Secondly, in the point force calculation the AL uses c*dr as local area element (where c is the local chord and dr the line segment), while the AD uses c*dr*Nblades/Ntheta, where Nblades is the number of blades and Ntheta is the number of points in the azimuthal (tangential direction). This approach is exactly as in SOWFA. In the UAD, the force at each actuator point is only axial and it is evaluated by using the conventional turbine thrust formula, scaled by dA/Adisk, where dA is the actuator element area and Adisk is the turbine disk area. This is explained in the paper, and plenty of references are given. Our implementation follows what is done in previous literature studies, so we believe that describing the models in detail is not necessary.

REV: Line 267 - Not all ADs account for blade rotation.

AUT: The reviewer is correct. We rephrased the first paragraph so that it is clearer we are talking about the actuator models implemented in TOSCA. The paragraph now reads “To represent wind turbines, different models have been implemented. In TOSCA, they are referred to as the actuator line (AL), actuator disk (AD), and uniform actuator disk (UAD) models. Following […], the first two models require detailed blade information (i.e. airfoils, twist and chord), while the UAD only requires turbine thrust coefficient and general rotor information such as diameter and hub height […]. The idea behind actuator models is to represent the wind turbine as a distribution of points, each associated with a Lagrangian force. For the UAD model, the sum of forces from all points must be equal to the total wind turbine thrust, while the AL and AD models in TOSCA additionally include rotor torque, as they also model blade rotation.”

REV: Line 286 – It is not clear what is meant by “blade loading is uniform:” is this uniform thrust (and torque) along the blade?

AUT: As mentioned earlier in the paper, the UAD model is different from the AL and AD models in which the blade info is accounted for. In the UAD model, only thrust coefficient is considered, such that the thrust is uniform along the blade and torque is null. If one would plot the radial distribution of blade loading, it would be constant, as the sampled velocity is the same for all points. 

REV: Line 368 – References for NCAR-LES and Wire-LES should be provided.

AUT: References added. Thank you.

REV: Line 403 – A long review of precursor methodology is given but no reference is provided.

AUT: References added.

REF: Line 407 – The question is if the plant induced gravity waves are a numerical artifact or they are real.

AUT: We believe the observed gravity waves are consistent with what is published on the characteristics of gravity waves in literature.  Their amplitude is dependent on the disturbance and on stratification, and their scale (in terms of wavelength) is extremely large and also depends on stratification. All this has theoretical and mathematical support (Gill 1982, Nappo 2012, Lin 2007, Smith 1980, 2007, 2010, Allaerts and Meyers 2018, 2019, Lanzilao and Meyers 2022a,b,2023).

We included an addendum that clarifies such concept through a simple but insightful model.

REV: Line 434 – The simulations were conducted for 120000s, or more than 33 hours, it is very unrealistic for the atmosphere to be steady for 33 hours.

AUT: The spin up phase itself requires more time than what a CNBL would be observed in real life. The use of the proposed controllers allows to run such idealized flow (in terms of its temporal evolution) in truly statistically steady state conditions, while e.g. not fixing the average temperature profile would produce a physical, but unreal because too slow, evolution of the ABL height. This whole concept is explained in Sec 2.5.2.

REV: Line 442 – It is not clear what is meant by “real flow physics.”

AUT: Rephrased to “Nevertheless, the proposed hybrid method, sketched in Fig. […], is very convenient as it allows to reduce the overall computational cost of the ABL spin-up phase, where wind farm-induced gravity waves are not yet present. In fact, this initial phase is run on a domain whose size is dictated by the current flow physics, rather than on quantities that will only become relevant at later simulation stages.”

REV: Line 448 – It would be better to reference the journal publication instead of the conference preprint, i.e., Martinez et al. 2018.

AUT: Revised also according to the first reviewer comment.

REV: Line 458 – It would be appropriate to provide a cross-section of the grid so that one can see the shape of the grid.

AUT: We do not believe it would add much information. Details are provided to exactly reproduce the mesh for anyone that is interested.

REV: Line 460 – Smagorinsky coefficient for standard Smagorinsky model is between 0.18 and 0.23 (cf. Lilly 1966 NCAR Tech. Note).

AUT: Our objective is to compare with Martinez Tossas (2015), so we set up the simulation in such a way that the two are exactly comparable.

REV: Line 462 – It is not clear what is meant by 72 points in azimuthal direction.

AUT: This is the discretization along the direction that is tangential to the disk perimeter. This is an angular discretization, as its actual size changes with radius, that is why sometimes is referred to as the azimuthal direction, for example in the SOWFA source code.

REV: Line 478 – It is not clear why is nearest neighbor approach used, it should be easy to implement interpolation to the geometric center of the Gaussian function.

AUT: Yes, this will be added in the future. With the revised data from Martinez Tossas (2015) data are more closer with TOSCA, so it does not seem that the interpolation makes a big difference.



REV: Line 487 – Why are sinusoidal perturbations added, why not random?

AUT: Random perturbations do not satisfy continuity and would be killed by the Poisson iteration. We follow the same approach used in SOWFA, where the introduced oscillations are sinusoidal and divergence free.

REV: Line 488 – These are surprisingly large wavelength perturbations- can the authors show why is that the best option?

AUT: We do not understand the comment of the reviewer in relation to line 488.

REV: Line 541 – The tower is just represented by drag - not actual actuator which is supposed to represent the rotor.

AUT: Yes this is the case. It is exactly the same approach, with the only difference that lift coefficient is zero. It is exactly the same as the rotor points that are located at the blade root for AL models, they have usually a lift coefficient of zero.

REV: Table 3 – It should be mentioned how S2 case differs from G case.

AUT: It is, see line 575: “In particular, we compare case S2 from Allaerts and Meyers (2017) against results obtained from TOSCA using both pressure and temperature controllers at the same time (case PT), and pressure and geostrophic controllers with no temperature forcing (case P and G, respectively).”

REV: Line 594 – 72000s is 20 hours. Again, is this realistic to expect stationarity for this long?

AUT: See previous comment and Sec 2.5.1 and Sec 2.5.2, this concept is explained.

REV: Line 612 – It is stated that “predicts lower geostrophic wind…” Isn't that just due to mass (i.e., more precisely volume) conservation?

AUT: This is in our opinion a consequence. Dissipation produces more mixing, and the inversion height raises. Then, due to mass (volume) conservation, geostrophic wind increases a bit.

REV: Line 614 – This is not due to “higher amount of dissipation” but due to mass (i.e., volume) conservation.

AUT: See previous comment.

REV: Line 639 – The streamwise spacing is smaller than used in any recently developed wind plant. The authors should explain why such unrealistic spacing is used.

AUT: We agree that the spacing is small, but not excessively in our opinion. Smaller spacing (down to 3.3D can be found in the Lillgrund wind farm for some wind directions)

REV: Line 641 – Is the 7.312 K a step change, or over some height?

AUT: Added “The temperature profile is initialized according to RampanelliZardi (2004)”. In Appendix A, where the whole set up is described.

REV: Line 646 – Low reflectivity may not be enough – reflected waves can be amplified.

AUT: We agree, but it depends when the reflectivity is calculated. Reflectivity should be always calculated on the latest available samples, as reflections usually build up during the simulation.

REV: Line 654 – The subgrid model plays a role in this, too, c.f. Mirocha et al. 2010 MWR, and Kirkil et al. 2012 MWR.

AUT: Yes, even though we did not observe such streakiness with lower values of equivalent roughness height.

REV: Line 671 – While the wind farm blockage may be amplified by gravity waves, it still exists in the absence of waves (Sanchez Gomez et al. 2023 WES)

AUT: We respectfully argue that reviewer’s comment is only partially true, that in the absence of gravity waves, wind farm blockage reduces to the superposition of the pressure increase created by individual turbines. If thermal stratification is modelled, the domain should be designed such that gravity waves can be captured, which is not the case in Sanchez Gomez et al. 2023 WES. In the absence of stratification, streamlines will be displaced more, as they do not encounter any resistance from buoyancy. As a result, induced pressure gradients – and consequently blockage - will be smaller (favorable at the farm exit and unfavorable at the wind farm entrance). This is because mass imbalance produced by velocity decrease is mainly compensated by streamline upward displacement. Pressure response is absent, as neutral stratification does not generate any density anomaly when streamlines are displaced. All this can be proved mathematically, as shown in the attached addendum. As a result, we believe that, without any stratification, blockage reduces to the superposition of individual turbine pressure increase, which has nothing to do with stratification. When stable stratification is present in the free atmosphere, large scale pressure gradients will be always produced to some extent, as they are a consequence of overtopping density anomalies generated by streamline displacement.

In conclusion, we do state that there is an effect that is not related to gravity waves in the sentence at line 649 “The gravity wave footprint inside the ABL can be clearly noticed in the pressure field, together with the small-scale pressure increase in front of each rotor”, but the dynamics is much more complex than simply asserting that blockage can be observed without waves.  

REV: Line 677 – What is the relevance of “an infinitely wide wind farm” for real wind plants?

AUT: We removed such comparison. It is misleading, as Wu and Porte Agel 2017 used a vertical domain size that didn’t allow to resolve gravity waves. On the other side, we believe that infinitely wide wind farms are useful, as they allow to simplify the physics while still providing insights, enabling model developments that could then be further extended to fully general cases.

REV: Line 693 – It is stated that “We also described a new methodology for simulating finite-size wind farm…” Strictly speaking, spanwise infinite wind plant is not finite-size.

AUT: Maybe the reviewer did not notice that a spanwise-finite wind farm is simulated in Sec 5, and the proposed methodology refers to this type of simulations (although is also applicable to infinitely wide wind farms). There are also Figures that show the simulated wind farm in relation to the domain (Fig. 9 and 10).

REV: Line 699 – The sentence “The concurrent-precursor domain is usually bigger than required, as its size is determined by the successor domain that runs concurrently” is not clear, why is the domain bigger than required?

AUT: We have added an explanation at line 443 (in the respective section) “Nevertheless, the proposed hybrid method, sketched in Fig. […], is very convenient as it allows to reduce the overall computational cost of the ABL spin-up phase, where wind farm-induced gravity waves are not yet present. In fact, this initial phase is run on a domain whose size is dictated by the current flow physics, rather than on quantities that will only become relevant at later simulation stages.”

REV: Line 702 – It is stated that “… TOSCA is able to capture wind farm gravity wave interactions and large scale blockage effects,” however, it is not clear that if the gravity waves are realistic – no comparison with observations is provided, so the term “capture” cannot be justified. This work just reproduced some of the previous numerical results that have not been verified in comparison to observations.

AUT: We agree that our study lacks comparisons to observations, as it is difficult to observe gravity waves in real life given their small amplitudes compared with the boundary-layer turbulence and their massive spatial scale, and we are not aware of any project that was dedicated, with suitable equipment, to these kind of observations. As such, we concede that the word “capture” may be better replaced with “simulate”.

Regarding the physicality of the simulated gravity waves, we argue that their existence is supported by the cited literature. If streamlines in a stably stratified flow like the free atmosphere are perturbed upwards by an object lying inside the BL, waves will be triggered. Their amplitude is dependent on the disturbance and on stratification, and their scale (in terms of wavelength) is extremely large and depends on stratification. All this has theoretical and mathematical support (Gill 1982, Nappo 2012, Lin 2007, Smith 1980, 2007, 2010, Allaerts and Meyers 2018, 2019, Lanzilao and Meyers 2022a,b,2023). Please refer to the addendum for the mathematical explanation.