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. 

The specific comments are as follows.

1. Abstract

It would be appropriate to add one or two sentences to highlight the key findings from this work.

2. Introduction

The novelty of the paper shoud be better highlighted. Application of Cholesky decomposition approach to generate a 4D wind field has already been reported in some literature. 

3. 4D wind field generation

In Section 2.5, it would be appropraite to indicate the maximum relative difference between the simulated values and the theoretical values when doing the comparision. 

3. Lidar simulation with integration of 4D wind fields

It would be appropriate to add references for the equations that are not derived by authors.

4. Other comments

It would be appropriate to present the Nomenclature and Abbreviations before the Introduction section. 

The structure of the paper could be improved. It would be appropriate to have a separate section to present Results and Discussion. 

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 computationally-efficient 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.

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. 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. 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.

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.

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".
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.
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?
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.
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.
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?
7. Pg. 5, ln. 117: "and a FFT factor". Please explain what the "FFT factor" is.
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.
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. 
10. Pg. 6, ln. 133: "gamma_{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. And it would also be helpful to explain what you mean by "magnitude coherence" (vs. "magnitude squared coherence").
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.
12. Eq. 10: Typo in the equation: the lower right term should have "x" instead of "yz".
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).
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.
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?
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?
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.
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?
20. Table 1: How are these parameters chosen? For example, are the turbulence parameters the default values recommended by the IEC standard?
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?
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.
23. Eq. 21: Is the "r" in "v_losP(r,t)" supposed to be "s+r_0"?
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?
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.
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.
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?
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?
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)?
30. Pg. 17, ln. 402: Please define "U bar" (unless if it was previously defined)
31.  Pg. 18, ln. 407: Please discuss what you mean by the "critical wavenumber"
32. Pg. 18, ln. 420: "between both coherence curves". The plots are showing the auto-spectra, not the coherence, right?
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 Delta s_k <= 1 /(2*k_max) to be more conservative?
34. Fig. 6 and Fig. 10: Please explain what discretization and total number of points are used to calculate the "theoretical" curves.
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.
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.  
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?  
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.
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?
40. Eq. A2: Please cite where this formula comes from (e.g., IEC standard)
41. Pg. 25, ln. 559: "where d_{yz} is a matrix…": More accurately, d_{yz} is a specific element of a matrix.
42. Pg. 26, ln. 563: "where sigma_i is the variance… sigma_iso is the…": These are described as variance here, but since they are not squared, they are only the std. dev. values.
43. Pg. 26, ln. 563: "l is the scale parameter": How is this parameter selected? Is there a default value you use?
44. Table B1: In the Simley and Pao and Kristensen models, can you describe what the parameters "sigma" and "L_u" represent?