Authors’ response: 4D wind field generation for the aeroelastic simulation of wind turbines with lidars

First of all, we would like to all the reviewers for their time taken to read our manuscript and their constructive comments. We have considered all the comments in detail and revised our paper accordingly. We believe that these comments have helped us to further improve the quality of our paper. Please find below our responses to reviewer comments. The reviewer comments are repeated in black text, our responses are given in blue text, and if necessary, the corresponding revisions are provided in red text. This extends the 3D method to 4D to enable the modelling of wind evolution along the wind direction. A two-step Cholesky decomposition approach is proposed for the factorization of the coherence matrices in the wind field generation. Overall, the paper is well written.

This paper extends the 3D method to 4D to enable the modelling of wind evolution along the wind direction. A two-step Cholesky decomposition approach is proposed for the factorization of the coherence matrices in the wind field generation. Overall, the paper is well written.
We would like to thank the referee for the interest in this research and the positive feedback on the manuscript.

Specific comments
1. Abstract It would be appropriate to add one or two sentences to highlight the key findings from this work.
Thanks for your comment. All the key findings have been properly summarized in the abstract.

Introduction
The novelty of the paper should be better highlighted. Application of Cholesky decomposition approach to generate a 4D wind field has already been reported in some literature.
Thanks for your comment. One of the novel findings of this work is the application of a two-step Cholesky decomposition, which is applied differently from the commonly used one.
3. 4D wind field generation In Section 2.5, it would be appropriate to indicate the maximum relative difference between the simulated values and the theoretical values when doing the comparison.
Thanks for your suggestion. In fact, the simulated time series follows Gaussian process and each Fourier coefficient in the frequency domain is Gaussian. When we simulate a Gaussian process with a large number of samples, the ensemble mean of the samples will be close to the theoretical mean value (please see: Law of large numbers). In our application, if we simulate the wind field with a larger number of samples, the estimated power spectra, cross-spectra, and coherence will be closer to the theoretical value. Also, the estimated power spectra, cross-spectra, and coherence depend on the method of spectra estimation, e.g. what type of window we use, whether the overlapping is applied, and the Fourier transform length etc. Overall, we think it will not be representative to show the errors between simulated values and the theoretical values, because one can reduce the errors by increasing the number of random samples.
4. Lidar simulation with integration of 4D wind fields It would be appropriate to add references for the equations that are not derived by authors.
Thanks for your comment. The authors have went through the paper and added references for the equations not derived by us. Specifically, these equations are Equation (1), (28) and (42).

Other comments
It would be appropriate to present the Nomenclature and Abbreviations before the Introduction section.
Thanks for your comment. We think it would be easier for readers to find variable names when we present the Nomenclature in the same section where the variables are used.
The structure of the paper could be improved. It would be appropriate to have a separate section to present Results and Discussion.
Thanks for your suggestion. Our manuscript discusses two main topics: the 4D wind field generation approach (an extension of the Veers method) and its application in combination with lidar simulations. Considering that some readers might be interested in only one of the two topics, we've therefore chosen this parallel structure and presented the results of the two topics in the corresponding section. This is an interesting and overall well-written manuscript describing a 4D wind field generator that introduces wind evolution to standard 3D wind fields via longitudinal coherence for realistically simulating preview-based lidar-assisted controllers. The most valuable contributions of the paper are: 1) providing a computationallyefficient method for computing 4D wind fields using the two-step Cholesky decomposition method (given some assumptions) that significantly reduces the computational effort when computing the Cholesky decomposition for large matrices; 2) describing how existing 3D wind fields can be used to create the evolving 4D wind fields, which makes it more computationally efficient and easier to integrate the 4D wind fields into existing aeroelastic simulation workflows; and 3) providing an open-source tool for computing the 4D wind fields. This is a relevant topic because of the need for easy-to-use, realistic wind field simulations to assist with the design and certification of wind turbines with lidar-assisted controls.
We would like to thank the referee for the interest in this research and the positive feedback on the manuscript.
While Section 2, describing the 4D wind field generation technique is very strong, parts of Section 3 on lidar simulations in 4D wind fields feel disconnected from the rest of the paper and could be improved. Specifically, Sections 3.3 and 3.4 look at the impact of range weighting function discretization and interpolation methods on the accuracy of the lidar measurement auto-spectrum, but do not consider the impact on the coherence between the lidar measurements and the rotor effective wind speed, which is just as important for lidar control applications. It would be good to include the impact on measurement coherence in Sections 3.3 and 3.4 and consider some real-world examples involving rotor effective wind speeds and lidar scan patterns.
Very grateful to your suggestion. We'll add the contents requested by the reviewer. More specifically, we'll first include the formulas for the coherence between the rotor-effective wind speed and the rotor-effective wind speed estimated by lidar (hereafter referred to as lidar measurement coherence) in Section 3.2. Then, in Section 3.3 and 3.4, we'll add the analysis of the impact on the lidar measurement coherence using a scanning configuration from a typical commercial lidar as an example. The lidar measurement coherence is compared using different range weighting function discretizations and interpolation methods.
Further, Section 3.3 discusses the "critical wavenumber" or "maximum relevant wavenumber" as a way to determine an acceptable spatial resolution of the range weighting function, but real-world examples of critical wavenumbers are missing, so it isn't clear what spatial resolution would ultimately be acceptable.
Thanks for your suggestion. We'll add examples to make the results more intuitive.
There also don't seem to be any conclusions drawn about the number of points that should be used to approximate the range weighting function.
The conclusion is that the number of points is not relevant but the step width (see Pg. 18 ln. 419). Therefore, we suggest to consider step width rather than the number of points when applying a discrete range weighting function. We'll highlight this conclusion in the corresponding text.
Another general comment is that since once of the major contributions of the paper is a computationally-efficient 4D wind field generation method, comparing the computational time of the method to the computational time of a single 3D wind field or a 4D wind field without the two-step Cholesky decomposition method would be very interesting to most readers.
Thanks for your suggestion. If you mean the comparison between the 4D wind field generation with and without the two-step Cholesky decomposition approach, corresponding to Eq.(14) and Eq. (7), respectively, we did consider to do it at the beginning. We didn't end up doing it because the algorithm that we've implemented in our tool evoTurb is mainly based on Eq.(17). This implementation improves the computational efficiency mainly due to the possibility of using pre-generated 3D wind fields, which is made possible by the two-step Cholesky decomposition approach. Even if we don't consider the use of pre-generated wind fields, evoTurb first calls TurbSim (Fortran based code) or Mann turbulence generator (C++ based code) to generate 3D wind fields and then finishes the rest of the steps in Matlab or Python. When it comes to comparing computational time, it is very important to keep the tools in the same code environment. In this sense, it is not justified to compare the computational time of evoTurb with an implementation of the direct generation of 4D wind fields (Eq. (7)) in Matlab or Python. Nevertheless, we can add a short discussion about the computational time of applying Cholesky decomposition to a huge matrix compared to two small matrices on a theoretical level, if it is interesting to readers.

Specific comments
1. Pg. 2, ln. 25: "Veer's method": Here and throughout the paper, since the author's name is "Veers", not "Veer", this should be referred to as "the Veers method".
Thanks for pointing this typo out. This has been corrected throughout the manuscript.
2. Pg. 2, ln. 25: "simulates stationary and homogenous multidimensional random processes. . . ": The Veers method does not require the turbulence to be homogenous in that it allows different auto-spectra and different TI values at different grid locations.
Veers (1988) mentions that the 3D wind field simulation method is based on Shi-nozuka and Jan (1972), and the method presented by the latter is for homogeneous process. If this is disputable, we'll remove the word "homogeneous". The references mentioned above: 3. Pg. 2, ln. 51: "Wind evolution refers to time-dependent variation of turbulence structure": Can you be more specific about what you mean by "turbulence structure" here?
It means eddy structure here. This has been added to the corresponding text.
4. Pg. 3, ln. 72: "both methods only can generate unfrozen turbulence on two different vertical planes. . . " This statement is a little misleading, because the extension of the Veers method in Laks et al. could fundamentally be applied to multiple planes. So the method is not necessarily limited to only two planes. However, it is true that the authors did not attempt to simulate more than two planes.
In that sentence what we actually mean is that the formulas are not intended to be directly applied to multiple planes. We'll improve the formulation to avoid misunderstanding. 5. Pg. 3, ln. 84: A reference to the evoTurb Github repo would be appropriate when the software is introduced. There doesn't seem to be any reference until the "Code availability" statement at the end of the paper.
Thanks for your suggestion. Now we've added the Github link to the first mention of evoTurb. 6. Fig. 1: The axes are a little confusing here. The "positive" x direction is toward the right, but the arrow indicating the u wind speed component direction points in the opposite direction. Can you clarify the spatial and wind component coordinate systems?
Thanks for the careful review. Indeed the coordinate in the figure is misleading, the turbulence should usually be defined to propagate towards the positive x direction. We have updated the figure to match this point. And the spatial Cartesian coordinate (xyz) is defined based on the right hand rule. 7. Pg. 5, ln. 117: "and a FFT factor". Please explain what the "FFT factor" is.
Thanks for your comment on this point. For brevity, we didn't give the equation of A u explicitly in the manuscript because it is considered as a constant in all subsequent derivations. However, if this might cause confusion, we'll remove the word "FFT factor" throughout the manuscript and add the equation of A u to the corresponding part.
8. Pg. 6, ln. 130: "A u is the amplitude composed of the auto-spectrum of the u component". To obtain the desired auto-spectrum using the Veers method, the amplitude A u should be proportional to the square root of the auto-spectrum.
Thanks for your comment. Now, this sentence has been rephrased and the equation of A u is given explicitly: A u is the two-sided Fourier coefficient obtained from the auto-spectrum of the u component S u at this frequency with ∆f the band width of the frequency in Hz. A factor of 1/2 needs to be multiplied in the square root if S u is a one-sided power spectral density. For a specific frequency component, A u is a constant.
9. Pg. 6, ln. 133: "Cholesky decomposition (Press et al., 1992)": Here and throughout the paper, terms don't need to include the citation every time they are mentioned. The citation to Press et al., 1992 probably only needs to be included when the Cholesky decomposition is introduced. This comment applies to other terms in the paper as well, for example Veers method and Kronecker product.
Thanks for your comment. We've removed citations for terms except for the first mention.
10. Pg. 6, ln. 133: "γ u,yz,i,j ": It would be helpful to define the basic coherence formula when this is introduced for readers not as familiar with coherence.
Thanks for your suggestion. Now, a general definition for the spatial coherence has been included in Appendix A as follows: Before giving the formulas of both turbulence models, a general definition of the spatial coherence between the wind components i and j (i, j = u, v, w) at two spatially separated locations k and l is given as follows: with f the frequency in Hz, S i,k (f ) and S j,l (f ) the respective auto-spectra and S ij,k,l (f ) the cross-spectrum.
And it would also be helpful to explain what you mean by "magnitude coherence" (vs. "magnitude squared coherence").
An explanation for the magnitude squared coherence has been added to Appendix B.
11. Pg. 7, ln. 178: "which will cause the Cholesky decomposition to be very slow.": I would suggest describing how the computational time of the Cholesky decomposition scales with the size of the matrix. This would be very interesting and strengthen the motivation for using the two-step method. I believe the computation time is roughly proportional to the cube of the number of grid points, but am not sure.
Thanks for your suggestion. We'll add a short discussion about it.
12. Eq. 10: Typo in the equation: the lower right term should have "x" instead of "yz".
Thanks for pointing this typo out. Now this has been corrected.
13. Pg. 8, ln. 193: "provides a very useful property". Is this property described by Press et al., 1992 or did you derive it? If it is from Press et al., then this would be one case where the reference should still be included. And if it is derived by Press et al., then I don't think it is necessary to provide the proof in the appendix. The appendix is probably only necessary if it was derived by you (the authors).
This property was derived by the authors. We've modified the text to make it clearer.
14. Pg. 10, ln. 233: "because the Kaimal model only considers the spatial coherence of the u component." Is this a realistic model for lidar simulation? For example, with no correlation between v or w components at different locations, volume averaging along the lidar beam could unrealistically average out the contributions of the v and w components to the line-of-sight velocity. If they were correlated, then line-of-sight errors would be larger (and probably more realistic). Some discussion on this point would be nice.
Thanks for the nice suggestions. We have added some discussions about the unreality of neglecting the correlation between v or w components at different locations in Pg. 10, ln. 235. The added content is: However, in reality, the atmospheric air flow is assumed incompressible for normal wind turbine applications. It implies that either the v or w component is spatial correlated due to the continuity of the incompressible fluid as discussed by Mann (1994)." And in Pg. 10, ln. 372, we've added some discussion on the impact of this unrealistic assumption on lidar simulations. The content is: As discussed in Section 2.4, it is unrealistic to ignore the spatial correlations of the v or w component at different locations from the physical point of view. The volume averaging of LOS speeds contributed by the uncorrelated v or w component can be unrealistically low because they are averaged out. In the case that the lidar beam is misaligned from the longitudinal direction significantly, further study is necessary to quantify the errors caused by neglecting the spatial coherence of the v or w component.
15. Pg. 10, ln. 245: Are the upstream planes time shifted to account for the longitudinal separation before they are saved as binary files in step 5?
Temporal shifts will not be applied when the upstream planes are saved. In-stead, the positions in x of different 3D fields (y,z,t) will be saved to make it easier to implement in the current OpenFAST lidar simulator (https://github.com/ sowentoDavidSchlipf/openfast/tree/f/lidarsim).
16. Pg. 10, ln. 248: "If the same wind fields already exist": Do you mean wind fields with the same random seed?
17. Pg. 10, ln. 253: "Export the 4D wind field as binary files": Are they exported in the same format as the original 3D wind field files?
Yes, the format of the wind field are the same. The identical scaling and offset can be applied to convert the binary data to decimal data (please see https: //github.com/OpenFAST/openfast/tree/main/modules/turbsim/src). But we actually separate the rotor plane 3D wind field (".wnd" extension) from the upstream 3D wind fields (".evo" extension) as two files. If lidar is not simulated or the Frozen theory is applied, only the ".wnd extension" file will be read in. If the lidar is simulated with turbulence evolution, then both ".wnd" file and ".evo" file will be read in.
18. Pg. 11, ln. 255: "additionally contains the spatial coherence of different v and w components": Does something like "spatial coherence of the v and w components" make more sense here? I.e., I'm not sure what "different" refers to here.
Thanks for the nice suggestion. After carefully reading, we thought the previous content here is not very clear. Now we have updated the text. Though the sentence is longer, we think the text is much clearer now.
19. Pg. 11, ln. 255: "and the coherence between the u and w components": But aren't you already accounting for this by correlating the w components (Eq. 19)? And couldn't that be extended to the v component as well to more realistically add wind evolution to the Mann turbulence fields?
Thanks for the nice suggestion. The previous text is not very clear, here we mean the coherence between u and w component at the same position. Actually, the u and w components are both spatial correlated, which means the u at position 1 is correlated with u at position 2 and the w at position 1 is correlated with w at position 2. In Mann model the correlation between u and w at position 1 or 2 is also modelled. Based on the spectral tensor of Mann model: Φ ij , it is possible to obtain the coherence between any two velocity component i, j in the wind field. In principle, it is possible to extend the evolution coherence to v component, however, as we discussed in the same paragraph and in the conclusion, we have not find information from existing literature to support us to reasonably model the longitudinal coherence of v component. So, this remains to be further investigated in the future. 20. Thanks for the nice suggestion. We have added the source of these parameters now.
21. Pg. 12, ln. 286: "in Fig. 4d-f that after applying the longitudinal coherence. . . " Can you explain whether the time shift between the different planes is applied in this example?
We are grateful for the suggestion. No, the temporal shifts due to the longitudinal separations are not shown in Fig.4d-f in order to make it easier to observe the effect of wind evolution. 22. Pg. 12, ln. 289: "retains more eddy structure of the original wind field as shown in Fig.4c." This is hard to see. Figs. 4d-f all look pretty similar to Fig. 4a/d. Fig. 4d-f are supposed to look similar to Fig. 4a/d because they have been constrained to it, which means evoTurb can do its job properly. To explain the constraining process more clearly, the following text has been added to the paragraph: More specifically: Fig.4a and d are identical since the wind field at x = 0 m is inherently regarded as the reference wind field in the constraining process; the wind fields in Fig.4e and f are generated by constraining the wind fields in Fig.4b and c to a with the coherence at ∆x = 50 m and ∆x = 100 m at the same time. Because the smaller the longitudinal separation, the higher the coherence, Fig.4e is more similar to Fig.4a compared to Fig.4f, whereas Fig.4f retains more eddy structures of the original wind field in Fig.4c, e.g. the strong eddies at z of 100 m to 150 m in the first 10 s. 23. Eq. 21: Is the "r" in "v losP (r, t)" supposed to be "s + r 0 "? Substituting s = r − r 0 or r = s + r 0 in Eq.(21) is equivalent.
24. Pg. 15, ln. 328: "only the wind fluctuations are considered in the lidar simulations": What do you mean by this? E.g., the mean wind speeds are not included?
Thanks for the suggestion. We have modified the text. In fact, the mean value should be considered to simulate the LOS speed correctly. 25. Eq. 22: How is the normal vector of the beam direction defined? Does this equation need to be multiplied by -1 to make the math work out? For example, if the beam direction vector is pointing away from the lidar, but the line of sight velocity is positive when flowing toward the lidar, then the "-1" term may be needed.
Thanks for the suggestion. We have modified the text in the paragraph and added an equation for calculating the unit vectors. Also, a figure is provided to show how we define the coordinate system. We calculate the unit vectors using the laser beam azimuth and elevation angles. The azimuth angle is defined as the angle between the negative x axis and the line that the laser is projected to xy plane. And the elevation angle is the angle between the laser beam and the xy plane. With this definition of coordinate system, the "-1" is not required. If the wind blows towards the lidar, the positive LOS will be obtained by the lidar. 26. Section 3.2: I think it would be worth discussing how the cross-spectra used in the derivations in this section are determined from the known coherence and auto-spectra for readers not familiar with this.
Thanks for the suggestion. We have added equations to connect between the coherence and the cross-spectrum in the paragraph. Now it is more clear how we can obtain cross-spectrum by knowing coherence and auto-spectrum.
27. Eq. 34: In reality, the wind speed would be estimated from multiple beams. Are you considering that in this analysis? Or is this intended to model a single beam?
Thanks for the suggestion. After adopting the reviewer's comment, we thought it is more interesting to investigate the rotor effective wind speed estimated by a lidar which is obtained by the signals from multiple beams. Please see the corresponding section for the updates.
28. Eq. 36: Similarly, the wind speed at the rotor plane is usually modeled as the average velocity over the rotor plane. Are you modeling this here? Or just investigating a single point?
Thanks for the suggestion. After adopting the reviewer's comment, we thought it is more interesting to investigate the rotor effective wind speed which is often defined as the mean longitudinal wind components over the rotor swept area. Please see the corresponding section for the updates. 29. Section 3.3: Is there a conclusion that can be discussed about the number of points in the range weighting function that is acceptable (3, 5, or 7)?
The conclusion is that the number of points is not relevant but the step width (see Pg. 18 ln. 419). Therefore, we suggest to consider step width rather than the number of points when applying a discrete range weighting function. We'll highlight the conclusion in the corresponding text.
30. Pg. 17, ln. 402: Please define "U bar" (unless if it was previously defined) Thanks for your comment. Now the definition of "U bar" is given where it is first mentioned. 31. Pg. 18,ln. 407: Please discuss what you mean by the "critical wavenumber" The critical wavenumber of the additional coherence γ s , denoted as k c , indicates the peak of γ s (see k 1 to k 4 in Fig. 7e and f). This explanation has been included in the corresponding part. 32. Pg. 18, ln. 420: "between both coherence curves". The plots are showing the auto-spectra, not the coherence, right?
Thanks for your comment. The "both coherence curves" means "curves of the magnitude of the FFTs (see Fig.7e)", which is equivalent to coherence. We've modified the expression here to avoid misunderstandings.
33. Eq. 40: This rule would still allow the first higher-frequency "peak" in the auto-spectrum to occur at k max . Would it make more sense to set the rule as something like ∆s k <= 1/(2 * k max ) to be more conservative?
Thanks for your suggestion. We've modified it to ∆s k < 1/k max . The suggestion we gave is considered as a minimum. One can either set a more conservative k max or add a factor to Eq.40 (as the reviewer did) to avoid "peaks" in the auto-spectrum. 34. Fig. 6 and Fig. 10: Please explain what discretization and total number of points are used to calculate the "theoretical" curves.
In Fig.6 the "theoretical" curves were calculated with the points from -60 m to 60 m with spacing of 0.1 m. In Fig.10 the "theoretical" curves were calculated with the points at -15 m, 0 m, and 15 m. 35. Pg. 22,ln. 472: "requires a massive amount of computational effort for 4D wind fields". But the two-step Cholesky method should reduce the computational effort significantly. Is the computational time still a barrier? As previously mentioned, this could be clarified by discussing the computational effort required for the Cholesky decomposition as the matrix size grows.
Thanks for the suggestion. We have modified the text and explain what is the number of 3D wind fields needed for the 4D wind field in the full-grid approach. The two-step Cholesky method indeed reduces the computational effort compared with the one-step approach. However, when we simulate a 4D wind field with many x positions, the computation time increases by the number of x positions. In the semifrozen approach, we reduce the number of x positions to shorten the computational time without violating the simulated lidar spectral properties. In the text, now we say we want to further reduce the computation time by introducing the semi-frozen method.
36. Pg. 22, ln. 479: "and thus wind evolution is negligible for such a short instant." I don't quite follow this argument. Are you saying that because the wind speeds along the probe volume are measured simultaneously, the wind evolution between those locations can be ignored? But the wind speeds on opposite ends of the probe volume are still separated by 30 meters, so shouldn't the longitudinal coherence still be considered? Also, the spectra from many pulses within the 0.25 second measurement period will be averaged to estimate the velocity, so isn't that a more relevant time period to consider than the 10 − 7 second pulse time? Lastly, if wind evolution within the probe volume should be ignored, then why would you consider simulating the full grid at all? I think it is fine to investigate the assumption of no wind evolution within the probe volume for this comparison, but more justification should be added if you are claiming that "Taylor's (1938) hypothesis is considered valid within the probe volumes." Similarly, if this is claimed, then more discussion should be included about how longitudinal coherence should be modeled when wind speeds are sampled at the same time at different longitudinal locations.
Thanks for the nice suggestion. We have modified this part. We think it is necessary to consider the "full grid" when the evolution is strong and the lidar probe volume is small, but for commonly used wind lidars in control application and most typical wind evolution conditions, the semi-frozen approach should give good approximations.
37. Pg. 22, ln. 481: "coherence between lidar-estimated u component in the upstream wind and the u component on the rotor plane". More details would be helpful here. Are you averaging the lidar measurements at the two range gates, and are you taking the rotor disk average of the u component at the rotor plane?
Thanks for the suggestion. After adopting the reviewer's comment, we thought it is more interesting to investigate the coherence between lidar-estimated rotoreffective wind speed and the actual rotor-effective wind speed in the rotor swept area. We will average the lidar-estimated u component and the u component within the rotor area in the future version.
38. Pg. 23, ln. 496: "LAC is a preview control based on the upstream. . . ": I would suggest "LAC is a type of preview control", or "LAC is a preview control strategy" here.
Thanks for your suggestion. The sentence has been modified accordingly.
39. Pg. 24, ln. 525: "coherence between the interpolated point and the corresponding neighboring points". Doesn't the formula only depend on the coherence between the neighboring points?
Thanks for the careful reading. The reviewer is correct, the weighting factors only need to be multiplied with the coherence between the neighboring points. 40. Eq. A2: Please cite where this formula comes from (e.g., IEC standard) The citation of IEC 61400-1:2019 is given in ln. 568 for the whole Appendix A. 41. Pg. 25, ln. 559: "where d yz is a matrix. . . ": More accurately, d yz is a specific element of a matrix.
Thanks for your comment. This has been modified accordingly.
42. Pg. 26, ln. 563: "where σ i is the variance. . . σ iso is the. . . ": These are described as variance here, but since they are not squared, they are only the std. dev. values.
Thanks for pointing this typo out. These should be σ 2 i and σ 2 iso in the text.
43. Pg. 26, ln. 563: "l is the scale parameter": How is this parameter selected? Is there a default value you use?
Thanks for pointing this out. After consideration, the authors now use the origi-nal form of the Mann (1994) model equations. Previously the equations are based on the IEC 61400-1:2019 where some of the equations have been made non-dimensional using the turbulence length scale parameter l. The reason to use the original form is that the original equations are easier to understand and track in a lot of existing literature. Now we have avoided introducing the parameter l.
44. Table B1: In the Simley and Pao and Kristensen models, can you describe what the parameters "sigma" and "L u " represent?
Thanks for your comment. Now an explanation of the variables has been added to the table caption.