The paper by Maurin ZOUZOUA et al. exhibits a demonstration study in which they propose a data-driven statistical model to better represent the turbulent heat flux exchanges in the atmospheric surface layer. The aim of this improvement serves as a potential alternative to the classical parameterization schemes that are inherently violated in weather and climate numerical models. The developed model uses approximately 10 years of near-surface sensible and latent heat fluxes over the supersite "Météopole" of the French national research infrastructure ACTRIS-FR, located in Toulouse, for training and testing.  

The authors selected thirteen variables to characterize the environmental forcing at an hourly time scale to estimate these near-surface fluxes (heat and latent fluxes). The 10-year observational data was split half-half between training and testing. The suggested model suffers from predicting negatively large fluxes, a signature of the stable regimes which current models suffer from too.   

Technically, the work is well executed. The title and the abstract are informative, succinct, and clear. The manuscript is well organized, logically and clearly conveyed. The authors took care of clearly explaining the motivation behind their work, depicting all the deployed measurement tools, and formulating the research topics to which the manuscript is answering. This is both instrumental for understanding the scope and content of the manuscript, making its reading easy and pleasant.

Although I think some of my comments below are significant and I would like to see the responses, there are number of points that must be clarified before accepting it for publication.

I recommend acceptance with minor revisions.  

General comment: Some of the grammar/language shall be revisited

Specific comments:

Line 25: Mention also the first important source of bias.

Line 74: In section 5

3.1:

Why not splitting the training datasets in each stability category (Stable, Neutral, Unstable) and different friction velocity ranges

Line 157: (ii) instead of (iii)

3.2:

Coarse time resolution (half-hour averages data may mask physical information for smaller time scales, at what time step the simulations were run?)

Lines 201-202:

The most conventional meteorological variables, such as T and RH at 2 m agl are also available. The lowest level as mentioned is within 20 m (pressure level akin to sigma coordinate which dynamically changes with time), so half eta-level is still higher than 2m. Any interpolation here?

Although, you mentioned this here Lines 303-305:

Moreover, the meteorological variables are directly taken at the first half-eta level (M=1, around 8 magl) instead of diagnostic variables as much as possible.

But 2m and 8m are quite far especially under stable conditions… Maybe, you need to comment on this..

Line 202:

The data are stored at a temporal resolution of 3 hours …

How did you compare the instantaneous 3^rd hour snapshot with observations averaged in half-hour chunks?

Line 204:

Meanwhile, the surface data, mostly provided by ORCHIDEE, consist of time-centred mean over a 3 hours window…

Specify what surface data…

Lines 209-213:

Could you show in Fig. 3 the grid-layout that shows the real observation geographic coordinate and the 2 nearest grid cells considered in the analysis.

In this context, given the locality of space (dx=dy=20 Km) and the point observations collected from a tower at specific coordinates, have you done sensitivity analysis on the effect of spatial averaging? Also, any comments on the local dynamics forcings in the considered two grid cells compared to the observations? You considered the surface type aspect as a main criterion in selecting these 2 grid cells, but what about the topography effects and the local environmental dynamics (wind speed and directionality, temperature, …) when comparing outputs in these 2 grid cells to the observation’s dataset at that specific grid location?   

3.3.1:

What’s the criteria for selecting the training dataset?

Is it season independent?

Have you done any sensitivity analysis on choosing different dataset i.e., dataset convergence in terms of variability in the final weights and biases of MLP?

3.3.2:

Lines 298-300: Define these variable acronyms especially:

Eventually, 4 trigonometric temporal coordinates are added for the description of seasonal (dx, dy) and diurnal (hx, hy) cycles.  

General Comments

The paper by Zouzoua et. al shows a promising demonstration of a method for using MLPs to better understand the sources of uncertainties and errors in surface layer schemes. The work clearly outlines the method, and demonstrates it’s applicability at a representative site. The quality of the work is generally high, especially in setting up the background and motivation, however some additional polish, revision for clarity, as well as addressing a few critical questions, could improve the manuscript, especially in the presentation of the results from the demonstration site. The method has strong implications, and I greatly look forward to seeing this applied to a larger set of sites. 

Specific Comments

The introduction, background, and methods are clear and well written, setting up for a promising manuscript. They detail both the significance of the work, potential applications, and provide sufficient information for others to apply a similar method to alternative flux tower sites. One exception to this may be figure 1. While a schematic figure is a good addition to the work, I found it hard to interpret. Color coding was used but not described, and it was a challenge to understand the processes shown in the figure without frequent back and forth with the text. The schematic would be stronger if revised to show a bit more detail so that it could stand independently from the text. A schematic diagram has very significant potential and I appreciate the choice to include one, however the current iteration adds little to the text.

The results section could use some additional context and discussion, as well as revision for clarity. The connection between the background/motivation/introduction and the results gets lost at times. One small change that could assist this is a more specific naming of the cases. In section 5 in particular, I found myself frequently confusing the two MLP based fluxes and struggled at times to immediately understand the comparison being made. Perhaps assigning abbreviated case names (one for each of the four (or five if you count the different grid cells): estimated fluxes, observed fluxes, fluxes in the same environment, simulated fluxes) as well as text reminding what exactly they represent within the section could improve this clarity. The authors could use the naming already present in the figures in the text, for example, to have strong consistency. Language throughout sections 4 and 5 connecting back to the goals and motivation in the beginning of the paper would also help promote cohesion and make it easier to interpret.

Finally, there are a few differences between the simulations and the tower observations (and the MLPs based on them) that hinder comparison. While the authors do not necessarily avoid talking about them, the discussion on them is scattered throughout and could be enhanced with a more detailed and focused discussion. In particular answering: 

• How does the mismatch of temporal resolution (30 min vs 3 hour) affect the results, particularly since we would not expect fluxes to be stationary over 3 hours (especially during the mornings/evenings)?
• How effectively can we compare between 20 km grid cells with different (and heterogeneous) land cover and a tower with ~4m agl flux readings which is likely only reading a small area from a grassland?
• What is lost by neglecting the soil/surface temperature? Those should have a strong correlation with the sensible heat flux in particular.

Technical Corrections

20-21:... key drivers of atmospheric boundary layer (ABL) processes, such as …

25-26: Rephrase “... source of biases in simulations with the numerical models”

107-109: I would rephrase for clarity “... we propose an evaluation approach dedicated to a full numerical simulation…” 

395: I think it would be good to elaborate here, as the use of coarse grid cells with a different land use type, versus the real grassland with the flux tower site creates a mismatch in the comparison that deserves further explanation/attention as it could explain some of the results.

524: remove “well”: “... simulated heat fluxes are not very well sensitive”

Figure 1: Color coding should be described in the figure description. Also, cyan color on white is a challenge to read. More overall comments on figure 1 in section above.