This is a timely article that rigorously quantifies values of shear stress prediction using a range of methods.  This topic is important as shear stress is essential to the accurate prediction of sediment flux and sand dune evolution in a range of aeolian environments. The article provides a good starting point for further exploration of the rapid calculation of shear stress over dunes, particularly in the context of landform evolution numerical modelling tools.  

There are however some details that require scrutiny before publication:

Some key details missing in the CFD simulations sections (3.2.1.) that require additional scrutiny before publication.

Section 3.2.1.  – what were the discretisation schemes used in the CFD?

Section 3.2.1. – What are the computational boundary dimensions (Length and Height). Please make clear whether in this section that all calculations are 2D.

Section 3.2.1. - Where we the ‘dunes’ placed in the context of the upwind and downwind boundaries?

Section 3.2.1 – What z0 value was used in equations 7 to 11? How does this compare to the smallest cell size? (z0 must be below 50% of the bottom cell height as wall function cannot extend above this. See Figure 3 in https://doi.org/10.1016/j.atmosenv.2006.08.019)

Line 226 – What is meant by a uniform base discretization?

Line 228 – Please quantify what a ‘sufficient drop’ is.

Figures 2, 3, 5, 6, 8, 10 and 11 lack text within the legend. It was therefore difficult to confidently assess the validity of these results as well and the descriptions of the data contained within the results section and discussion. This was particularly pertinent to Figure 10 and the discussion in Section 4.3.

Discussion – The authors responsibly highlight that the results only relate two-dimensional structures. Some brief discussion should also be made that they are also relatively simple idealised structures as well.

Using CFD to test the validity of linear approximations of the hydrodynamic equations outside their intended parameter range is certainly a necessary, and sometimes laborious, undertaking. However, there are several important technical and conceptual issues with this work that should be resolved before eventual publication.

1. As shown in Charru et al, Annu. Rev. Fluid Mech (2013) (Fig.3a), it is known that the linear approximation obtained by Jackson and Hunt (1975), and Hunt et al. (1988) (Eqs. 1-3 in the manuscript), overpredicts the value of the constant A (Eq. 5a) for L/z0 < 10^8 compared to the full solution of the turbulent boundary layer in the limit of small amplitudes (H/L << 1). For the value L/z0 = 25000 considered in the manuscript, the KSH model predicts A around 5, whereas the full solution would give A around 4. The authors should add this to their description of the KSH approximation and the discussion of model accuracy.

2. There is a mistake in Eqs.5b and 5c, in the actual approximation by Kroy et al. (2002), the term 2*ln(pi/2) is neglected. This term comes from the term 2*ln(kL) in Eq.2, since for a periodic modulation of wavelength lambda we have k = 2*pi/lambda, and by definition L = lambda/4.  

2.1. For clarity, in Eq.2 the term ln(k) should be replaced by ln(kL) since all lengths are rescaled by L in this model.

3. The value L/z0 = 25000 used to evaluate the prediction from KSH seems arbitrary. If L = 50m and z0 = 10^-3 m, as seems to be the case for the CFD, then L/z0 = 50000. Furthermore, given the known overprediction of this approximation, the authors should add another prediction using z0=10^-4 m, which is closer to the actual roughness of sand. 

4. There is a major confusion with the dune aspect ratio. In Fig.1, the aspect ratio H/L defines L as the dune's toe-crest length, whereas in all models L is defined as the half-length at half-height (the length between a point at half height and the crest, if I understood correctly). This means that there is no obvious way to compare the different simulated cases (H/L) with the relevant values in Fig.1. Here I have several suggestions:

4.1 Please use different symbols for the different lengths, for example keeping L for the half-length at half-height but using L_base for the toe-crest length. You could also add a diagram illustrating the different terms to avoid confusion. 

4.2. Given the toe-crest length is around twice the half-length at half-height, you could add another axis in Fig.1 to convert the standard aspect ratio plotted (using toe-crest length) to the new one used for the simulations and all other figures. This is crucial because the simulated range in H/L = [0.01, 0.5] only represents up to H/L_base = 0.25 in Fig.1, thus missing many relevant cases. 

4.3. You should consider extending the simulated ratio H/L (using half-length at half-height) all the way up to the avalanche angle (H/L = 1.3 and H/L_base = 0.65) to include the whole relevant range of foredunes' slopes. Also, please add the avalanche slope in Fig.1 for reference.

5. There is a problem with using a Bump as a dune profile to compare to CFD. Because a Bump, as defined in Eq.6d, is not a smooth profile (there is a discontinuity in the slope and/or curvature at the base) is will experience flow separation at the lee-side and flow stagnation at the toe. Therefore, the predicted shear perturbation from KSH is not defined at the base (tau' < -1) and has to be replaced by tau' = -1, which means a vanishing bed shear stress (tau = 0). Therefore, there is no meaning in using the values at these extreme locations to evaluate model accuracy (Fig.3e,f). Furthermore, the CFD predictions for the Bump (Fig.3c,d) don't really show the stagnation and flow separation expected at steep slopes, which suggests an issue with the simulations. The authors should either avoid the Bump case (and restrict only to smooth profiles) or properly discuss the complex physics involved.

6. The symbolic regression result shown in Eq.19 is physically wrong, since a constant term in the relation between the Fourier transforms of the bed shear stress perturbation and the surface implies that the shear perturbation is proportional to the surface itself which doesn't make sense (the flow field is essentially scale-invariant in the limit of large L/z0). 

6.1. There is a similar problem with Eq.20, which is quadratic on the wave number k. There are strong physical reasons behind the form of Eq.5a. If the authors want to improve the KSH model, the simplest way is to fit the parameters A and B to the CFD results on smooth profiles. 

7. The way to measure the differences between the models, using MSE (Eq.18) is not particularly useful because it amplifies minor horizontal shifts in the simulations. For example, even a perfect prediction shifted horizontally by 1m would lead to a large MSE. You should also include the difference in the maximum (and minimum) shear stress perturbations between the models. 

8. The examples in Fig.12 for the predicted transport rate seem completely arbitrary. The linear models have been successfully used to simulate unvegetated dunes in many different situations mainly because the aspect ratio of dunes (using toe-crest length) is always less than ~0.15. Steeper dunes have to be vegetated, but then there is essentially no sand transport. The authors should either use actual dunes shapes or at least the aspect ratios of actual mobile (unvegetated) dunes to compare sand transport rates.

8.1. The transport equation (Eq.21) is wrong for aeolian sand transport. Please use the one in Martin and Kok, Science Advances (2017), which is consistent with several wind tunnel measurements and transport simulations.

8.2. The data in Fig.12f seems inconsistent with Fig.3b. In Fig.3b, both CDF and KSH are quite similar for H/L = 0.4, but they are quite different in Fig.12f for H/L = 0.5.   

9. In lines 100-105. It is not true that the KSH model cannot be applied to scarps. In fact, a better version of this model (given by Weng et al. 1991) is already implemented in the Coastal Dune Model which certainly simulates dune recovery after erosion and scarp formation (Duran and Moore, Nature Clim. Change, 2015). In this case, similarly as for the Bump, the shear stress perturbation reaches < -1, and the only physically meaningful condition is tau = 0 (which follows from flow stagnation).

10. In addition to the technical issues mentioned above, I see two important conceptual problems. The manuscript introduces the main limitation of the linear models (e.g. KSH) using the steep slopes of the foredunes or vegetated dunes at the beach. However, they never compare the models using an actual foredune profile. Furthermore, the relative difference of the maximum shear stress perturbation between the KSH and the CFD for a gaussian hill seems to improve for larger H/L (Fig.3a,b) instead of getting worse, which would undercut the main purpose of the manuscript. There is also a consistency issue here because steep foredunes (where model's predictions would be worst) are covered by vegetation (that is why they can be so steep in the first place) and therefore there is no active sand transport on the crest and lee side. On the other hand, active dunes would have milder slopes and thus the model would be more suitable. These things should be properly discussed to help the reader understand the scope and meaning of the actual problems with flow simulations.

11. The second conceptual problem is related to the general assumption behind the manuscript, that somehow there should be a simple shortcut to improve the calculation of the bed shear stress under arbitrary topographic conditions. These models (in particular the 3D version in Weng et al. 1991) are already very powerful since they allowed for the first time physical simulations of 3D dune formation and interactions, formation of dune fields (Duran et al., ESPL, 2010), the formation of parabolic dunes (Duran and Herrmann, PRL, 2006), coastal dunes, dunes on Mars (Parteli and Herrmann, PRL, 2007), etc. This implementation of the Weng et al. (1991) approximation is very fast to compute and works (with the proper extension for the flow-separation region) for arbitrary topographies. Of course, it is not perfect and it is important to evaluate its accuracy, but it should be clearly stated that it represents a actual physical approximation, not a phenomenological fitting, and cannot be easily replaced by a data-based approach that requires many 3D CFD simulations to provide data. I would suggest to focus on using CFD to improve the parameters in KSH as a way to move forward in the field without having to start from scratch with data-based approaches after +40 years of physical insight and meaningful applications.