I think the work has significance for the scientific community. As described by the authors, similar development has been done in the past, but the present work removes some of the assumptions behind those, to get higher accuracy/realism in the simulations. This is well explained in the text, which is also a summary of the VG modelling by meaning of panel code and highlights the next steps to further improve

The paper is nicely written and easy to read. First, a CFD RANS simulation is made over a flat plate with VGs to investigate the boundary layer parameters for zero pressure gradient. Based on this some functions are made that can with some success model the spanwise variation of the shape factor, H, see Figure 6a and an excellent agreement of the spanwise averaged value, Figure 6b. This is only shown at a distance of 10 heights from the VGs and it could be nice to also see if this iis equally valid at some other distances. It is a bit strange to apply a compressible solver even though this is capable with a preconditioner to model an incompressible flow. And, also a grid of 6 million cells is not impressive and maybe just at the limit even for RANS.

There is something wrong with the units in the extra terms in Eqs. 6 and 7, e.g. the extra term including the induced pressure from the VGs has dimension N/m^2 and the other terms are dimensionless in the equation. The last term in Eq 7 is dimensionless, but not the second last.

It is written that pi,VG is the induced pressure from the presence of the VGs, but it is not clear how this is computed. A lot more details on the coupling and the closure terms  are needed. And an important reference to Ramos-Garcia et al is missing, that also made some improvements to XFOIL and where results are quite good for clean airfoils, “A strong viscous-inviscid interaction model for rotating airfoils, Wind Energy 2014 (17)”

The first test case of a thin airfoil with a very low angle of attach gives reasonable results for the BL parameters, but the skin friction Cf I Figure 11c is quite off, also without VGs, but that is not contributed to the extra VG terms but the general IBL.

Then it was chosen to compare results against two old measurements, one made in the Velux tunnel and another using the TUDelft WT. The Velux tunnel data are old and of not so high quality and made at very high inflow turbulence and the TUDelft WT much better. None of the results are very good and overshoots Cl,max, but it is known that computing thick airfoils is difficult and wind tunnel measurements are also not easy since the flow is 3-D. But it is not so convincing that RFOILVG is much better than the older versions.   

The authors mention that there still is work to do with respect to include higher pressure gradients and I recommend that the work is not yet ready for publication in WindEnergyScience.

Overview

This paper proposes an integral boundary layer approach to modeling the effect of vortex generators (VGs) for use in reduced-order airfoil aerodynamic analysis tools such as XFOIL / RFOIL. The model is based on flat plate, zero pressure gradient (ZPG) boundary layers. The effect of VGs is modeled through shape factor and viscous dissipation. The model is implemented in RFOIL and is tested for polar predictions for a variety of airfoils. Vortex generator models in XFOIL/RFOIL already exist, however, they model the effect of VGs as an additional source of turbulence in the boundary layer. Such models are highly empirical, employing several tunable coefficients that require training and work well for cases that are used for training, but not as well as for unseen cases. The key contribution of the paper is that the new model is analytical (does not require tuning) and it accounts for changes in flow/momentum in the boundary layer and not just model the VGs as a turbulence source.

The idea pursued in the paper is good but I see a fundamental issue -- at high Reynolds numbers (relevant for wind turbines), stall is caused by turbulent boundary layer separation that progresses (with increasing $\alpha$) from the airfoil trailing edge to the leading edge. Boundary layer separation is determined almost exclusively by the adverse pressure gradient in such cases and hence models built on a ZPG flow assumption have inherent limitations. This can be seen in the results presented in the manuscript. I have the following suggestions/comments; the first and the main comment requires significant work and modification to the manuscript.

Questions/Comments

• A systematic study can be performed with the numerical setup discussed in Section 3 to investigate the effect of adverse pressure gradients. The upper boundary in the simulation can be modeled as a wall and its geometry modified to prescribe a pressure gradient distribution on the bottom flat plate with VGS. CFD simulations can be performed for varying magnitudes of adverse pressure gradient (APG) and the proposed model can be a) tested against the data, and b) potentially improved using the data. Given the critical importance of APG, I feel its inclusion is critical. Even if it is not included, a comparison against CFD for a systematically performed analysis of the effect of APG would add tremendous value.
• Comparing Figs. 12 and 13 in the manuscript shows a significantly increased discrepancy between the $C_l-\alpha$ curves with the VG models (not only the proposed model, but also the existing RFOILVG model). This increased discrepancy is seen in every case. What causes this? Surface $C_p$ distributions can be plotted to compare if the loss of lift is occurring due to a higher trailing edge flow separation, or something else.
• The authors claim in some places that the stall position and post-stall characteristics are predicted better with the proposed model. Looking at it from a physical perspective, the VGs should be ineffective (or at least its modeling would be highly inaccurate) post moment stall as the flow over the airfoil is fully separated. What is the rationale that the ``improved'' model is because of accuracy.
• A passing remark is made about improved robustness of the proposed model in the conclusion. This should be elaborated on and put in the main paper.