GENERAL COMMENTS:

This paper reports a very witty use of automatic calibration algorithms to highlight an unphysical flaw in surface heat flux parameterisation in the widely used Delft3D hydrodynamic model. Such issue is often not evident in common uses of the model, in which temperature calibration is performed over the full water depth and calibrated turbulence parameters compensate the error. The issue was thus tracked down optimising surface temperature only and blocking turbulence parameter calibration, taking benefit of the optimal solutions found through unbiased automatic calibration. A correction in the shortwave radiation absorption parameterisation is then formulated and included into the model, making it able to reproduce the micro-stratification often observed in the surface-most part of the epilimnion of open waters.

The research is well developed and mostly efficiently conveyed in the manuscript. The topic is relevant to the journal and the results are useful to the scientific community. I’m highlighting below some prompts for further discussion of results and of their implications. I’m also suggesting a careful final proofreading of the paper aimed at optimising general clarity.

Given all the above considerations, I think this paper can be published after minor revisions.

SPECIFIC AND TECHNICAL COMMENTS:

L32: “allowing for extending” is not proper English

L33: the meaning of “surrogate models” should be disclosed here

L55: “in occurrence with” is not proper English

L77-78: Delft3D-FLOW solves the shallow-water equations when used in 2D mode, the RANS equations with the Boussinesq assumption when used in 3D mode. This should be better specified.

L78-80: if I recall correctly, the vertical momentum equation is still present, yet it is simplified to a hydrostatic equilibrium condition

L92: the c_gamma = 1.7 factor actually dates back to the pioneering work of Poole & Atkins (Poole, H.H., Atkins, W.R.G., 1929. Photo-electric measurements of submarine illumination through-out the year. J. Mar. Biol. Ass. U.K. 16, 297-324)

L95: “in the near” is not proper English

L102: add “in”

L132-134: define “phi” at first occurrence

L135-136: although I get the meaning of this sentence, it should be clarified

L139: remove “’s”

L140: use “i” as generic subscript for added clarity

L163-164: improve the definition of “T*”

L187: replace “(14)” with “(12)”

L194: ok, but continuing your reasoning, also Eq. 15 drops to zero if D_s is very large. I get the meaning, but the reasoning and the relevant outcomes of this passage should be better explained

L228: specify that D_s with the hat is the one used in the model and the one without hat is the one from field data

L230: explain why you are reporting only the first values in Table 1

L243-244: what do you mean by “speed of convergence”?

L282-284: ok, but surface temperature measurements are always somehow elusive, especially under no-wind or weak wind conditions for which surface micro-stratification builds up, as there is no consolidated methodology on “how deep should I plunge down the thermistor to make a surface temperature measurement?”. This could be discussed.

L300-302: ok, but there are other reasons for which the Secchi depth parameter could still be calibrated, such as to lump intrusions by tributaries, baroclinic mixing effects not simulated by the model due to the Boussinesq assumption, imperfect turbulence modelling, imperfect atmospheric boundary forcing. This could be discussed.

L319-326: I think that this is the most relevant result of this work, the fact that by improving the model you are able to reproduce the typically observed micro-stratification in the top-most part of the surface mixed layer (Figure 4b). You should try to make this achieving more evident, starting from the abstract.

L334-336: if I got it right, this is close to what the GOTM (General Ocean Turbulence Model) 1D model already does, following the Paulson and Simpson parameterisation (Paulson, C. A., Simpson, J. J., 1977. Irradiance measurements in the upper ocean, J. Phys. Oceanogr., 7, 952-956.), in which also the near-infrared shortwave radiation has an exponential decay, yet over a much shallower layer than the ordinary shortwave radiation. This could be discussed.

L338-342: it should be pointed out here that automatic calibration is (mostly) unbiased, as opposed to user-sensitive manual calibration, and as such can better highlight objective flaws in models

This is an overall well-written and interesting manuscript demonstrating how optimization toolboxes can help uncover limitations in simplified models/parameterizations within complex numerical models, with an application to a widely used 3D hydrodynamic model. The authors show that in the present version of this hydrodynamic model, the parameters of the shortwave radiation module often must be tuned to unphysical values to achieve the best response in lake temperature. They convincingly explain that this is due to a simplification in the shortwave radiation penetration model and provide a physically based improvement to the code.

In general, I find this an interesting study worth publishing. However, some context is missing, and the wording should be adjusted in some places. The authors demonstrate for one, no doubt interesting and relevant case, that their framework can work well, providing a simple way to improve a model with a few lines of code. However, it is my understanding that other widely used 3D numerical hydrodynamic models already have a better implementation of the shortwave radiation penetration, e.g., those based on Jerlov water types, as implemented in the MITgcm and ROMS (see Paulson & Simpson, 1977; https://doi.org/10.1175/1520-0485(1977)007%3C0952:IMITUO%3E2.0.CO;2). Thus, the improvements the authors suggest are no revolution as such.

I am not saying the model improvements are irrelevant and I agree that they should certainly be implemented into all newer versions of Delft3D. However, the title promises a bit more (“improve the physics behind complex numerical models”), in the Abstract it says that “similar outcomes in other applications” can be expected, and in the Conclusions it says that “users and developers can leverage these toolboxes to detect code-related issues” and that “the insightful analysis of the results of the optimization may offer similar outcomes in other applications”.

Yet no further examples are given anywhere in the text. I am wondering if the authors can think of other areas where optimization toolboxes could help to improve lake/ocean models.

Another example of a rather general statement that might need to be clarified: in line 28 it says:

“Therefore, some processes that are crucial in lakes, but not for the original environment, may not be correctly reproduced.”

I agree that the purpose for which a model has been developed is something to keep in mind and can be important. However, in the case of penetrative shortwave radiation, some of the mentioned models already have a more advanced parameterization than Delft3D even though they have been developed for the ocean (e.g., MITgcm and ROMS; see comment above). Later in line 100, it says that “many numerical models” already employ a shortwave radiation model that includes beta. So, all in all, it sounds as if Delft3D is the outlier here. This does not imply that the results of this study are unimportant, but I think some context is missing.

Other minor comments:

- line 33: Briefly explain surrogate models?

- line 43: “Many examples in THE literature”

- line 43: “Many examples”: could you name a few?

- line 50:  You refer to Schwind et al., (2022) but do not say what their study was about or why we should care about it. Give key results?

- line 53: “has been long recognized”:  some references?

- line 61: In large lakes, D_s can also vary significantly in space (see, e.g., https://doi.org/10.1016/j.jglr.2020.03.013). Maybe this is worth mentioning? Also, “can significantly vary” is very general. Could you provide some typical values for D_s throughout the season?

- line 111-113: It is not clear here why non-penetrative terms should usually represent a sink/be negative. You get back to this later when introducing the concept of thermal equilibrium. --> Maybe mention here that this will be explained later or give a brief explanation already?

- line 125-126: Can you name other lake models where the beta term is absent so that colleagues using them become aware of these limitations? As far as I can tell, MITgcm and ROMS (two other widely-used models, although admittedly mostly in large lakes or the ocean) already include beta (see comment above).

- line 139: “water density” sounds more natural than “water’s density”

- lines 150-151: The wording is a bit unclear. Isn’t the “overall” lake temperature, as in the total heat content, always determined by the balance between the surface fluxes unless we assume geothermal heat fluxes or other sources/sinks, e.g., river inflows? Doesn’t equilibrium (surface) temperature rather mean that the net heat flux is zero?

- line 175: “second term on the right-hand side of equation (13)”: unclear which line of equation (13) you refer to. Introduce (13.1) and (13.2)?

- line 187: “expressed by equation (14)”: do you mean by equation (12)?

- line 212: typo in “correctly”

- line 213-214: It is unclear where and when exactly these temperature measurements were taken. What does “surface temperature” mean exactly? I assume measurements were taken at some depth near the surface. Could the depth of this “surface temperature” have an impact on the optimization? Maybe summarize the measurement locations and periods in a short section or table?

- line 219: You optimize based on full-depth or surface temperatures. As shortwave radiation penetration decays exponentially, I assume the upper few cm or m are most strongly impacted by the improved parameterization. Would it make sense to define a cost function that focusses on these upper layers? For example, by taking one e-folding length scale as the range of the cost function. Or do you expect the results to be insensitive to this?

- line 241: If you give the computation time, also give the number of nodes/cores.

- lines 255-259: Can you comment on when/why it is valid to neglect advection and diffusion and the role of the wind speed? I assume low wind means less wind-induced turbulent mixing and advection, so vertical temperature gradients due to wrong surface forcing can build up more easily, driving stronger “artificial” convection in the model to compensate for that.

- line 276: “fully consistent with literature values”: Could you please also give these literature values here and not only in the caption of Figure 3?

- line 287: “ON three selected days”

- line 287: Why did you select these days, and does that choice matter for the results? For example, do you expect different results for calm vs. windy conditions (see comment above on wind-induced advection and turbulent diffusion)?

- lines 329-336 (see also lines 290-291): Does this imply that the level of improvement one can expect from the improved representation of shortwave penetration directly depends on the vertical grid spacing near the surface? If so, this seems worth mentioning/discussing further.

General minor comments:

- Not sure if this is on purpose or just a formatting error: When giving the units, there is a full stop/period instead of a space between units, e.g., W.m^-2, instead of W m^-2

- Write “in Appendix A/B”, instead of just “in A/B”?

- Consistent AE, BE spelling (modelling and modeling are both used)

- There are a few minor English mistakes, as pointed out by the other reviewer. I trust that the authors will correct them in the revised version.