The authors thoroughly investigate the effect of blade flexibility on loads and wake behaviour using LES simulations with a flexible actuator line coupled to a structural solver. Their results and conclusions present a significant contribution in their field and they are presented with great quality. Overall it is very nice to read with maybe a few wording issues here and there. 

There are few comments that the authors should address in the attached document before acceptance, mainly concerning the actuator line implementation, numerical sensitivity and analysis of the flowfield.  

“Investigation of blade flexibility effects over the loads and wake of a 15 MW wind turbine using a flexible actuator line method” investigates the effects of blade flexibility on the forces and wakes of the IEA 15-MW wind turbine, using a LES simulation with the actuator line method (ALM) coupled with a blade structural dynamics solver. Three models for the turbine are considered: rigid undeformed blades, rigid deformed blades and flexible blades. For the case with uniform inflow, with and without turbulence, the results are compared to the results of a BEM (blade-element momentum theory) model and a free vortex wake method. The turbine is also simulated in a neutral atmospheric boundary layer.

The detailed analysis of this work includes the analysis of the distribution of forces and displacements using different methods and flexibility models, in time and frequency domains, in addition to characteristics of the wake. It is a clear contribution to the field of research and within the scope of the journal. The manuscript is well-written, with clear objectives. The results are original, with associated discussions of high scientific relevance. There are a few remarks that should be addressed before publication, please see below.

1. In lines 93 and 465: reference to “advanced” ALM. However, this concept of “advanced ALM” is ill-defined. For example, previous works (Churchfield et al., 2017) have defined it by a collection of developments, however, not all of these developments seem to be implemented in the present work. Also, “advanced” may be interpreted as a common adjective, because it is not capitalized. In that sense, the concept of “advanced ALM” is questionable in the present context. For these reasons, I suggest that the term “advanced” be avoided.

2. In line 107: The mentioned corrections are not simple tip corrections, a more modern term is “smearing corrections”, because they correct for the smearing of the forces, not only at the tip. If you prefer to be more specific to this class of smearing corrections based on vortices, one usual term is “vortex-based smearing corrections”.

3. In line 111: The use of “much” in "much decreasing the need for a tip correction" might be misleading. It decreases the need for a correction, but not by orders of magnitude. Despite having a lower error than a 3D Gaussian, Caprace et al. (2019) showed that there is still a relevant error if using the 2D Gaussian without any correction.

4. In line 181: A wake length of 10 rotations without a far wake model might be too short. Please show that this choice of parameters does not have an effect on the distribution of forces (results of table 1 and figure 4).

5. In figure 4: The results from OLAF show a small peak near the tip, which looks unusual. Please explain it and show that it is not dependent on the choice of parameters.

6. In line 202: Please define “lateral forces”.

7. In line 205: Please cite (Caprace et al., 2019) when mentioning the error of the mollification of the ALM. That work clearly showed an error in the induced velocity for the choice of smearing used in the present work (2D Gaussian).

8. It is relevant to cite (Mikkelsen, 2004) in the introduction and in the discussion of section 3.1 (lines 202 to 206). This work showed, in 2004, that an uncorrected 2D Gaussian had better results than an uncorrected 3D Gaussian. Also, the same work showed that the ALM with a 2D Gaussian regularization without tip correction predicted higher forces than methods (ALM and ADM) with Prandt’s (or Glauert’s) tip correction. (Mikkelsen, R.F., 2004. Actuator disc methods applied to wind turbines, PhD Thesis, DTU, Chapter 8).

9. In lines 230 to 235: It is mentioned that the blade natural frequencies are calculated using the undeformed configuration. However, the flexible blades oscillate around the mean deformed configuration. For a non-linear solver such as BeamDyn, are the differences in the natural frequencies negligible?

10. In line 269: It is not clear that the value of 20 rotations for the discarded time is sufficient to reach a statistically stationary flow (it seems to be much lower than one flow through time). Please show that this choice does not affect the statistics.

11. In line 394: Analogously to the previous remark. Please show that this choice of “last 150 rotations” does not affect the statistics. How many flow through times were discarded?