This paper examines the ocean radiant heating parameterization and its impact on prognostic schemes simulating the diurnal cycle of SST.  Typically, in such prognostic schemes (e.g., Zeng & Beljaars (2005)) the diurnal warm layer depth is fixed at 3m.  This assumes negligible diurnal warming and minimal radiative heating below this depth; however, this paper notes that this is not necessarily the case and suggests that the warm layer depth in the prognostic scheme should instead vary spatially and temporally depending on the e-folding depth derived using chlorophyll data. This seems to me to be a worthy improvement and an interesting study.

The default parameterization for solar radiative heating in such schemes uses a 9-band decomposition of the solar radiation with invariant coefficients. However, it is well known that water properties can influence solar penetration; hence, the development of parameterizations with improved skill that make use of ocean color data. For example, Morel & Antoine (1994), Ohlmann & Siegel (2000), and Lee et al. (2005).  These types of schemes are not mentioned in the manuscript, but I think provide helpful background.  How does the new scheme compare to these and why not implement one or more of the above schemes for comparison? 

This study uses code from established and widely-accessible models.

Specific comments and technical corrections:

Line 19-20: Clarify this sentence – as the new scheme also truncates solar transmission! The difference seems to be between a fixed warm layer depth vs a spatially and seasonally adjusted warm layer depth.

Line 23: “… improves the diurnal signal …” compared to what? When you refer to the “old one” which old one as there are 2!? It seems to me that modradnemo is clearly better than wrfskin, but similar (or even worse) than nemoskwrite.  What is the evidence/argument from the results that modradnemo is better than nemoskwrite?

Abstract: “… we develop a new scheme …” “Our “chlorophyll-interactive” method”.  Can you clarify, what is the actual new contribution here? Is it to try out (implement) an already available method within NEMO (as I understand from lines 270-276) within a prognostic scheme for skin SST?

Line 67: suggest rewording

Line 77: should be 2001 not 2000

Line 79: “… by most of existing …” --> “… by most existing …”

Line 100-101: suggest rewording final sentence: “We redirect the documentation …”?

Line 125: I don’t think “resumes” (and “resuming” in the table caption) are the correct or best word choice, perhaps “outlines” (and “outlining”).

Line 149: remove “at day” (or possibly replace with “today”)!?

Line 191: Do you consider sea surface albedo? Is R_s net; i.e., after albedo reflection?

Line 205: reword “(already known from decades at those times)”

Line 215: “neglecting the effect of solar radiation”.  Saunders (1967) did in fact recognize solar radiation within the skin layer, see their equation (6).  It’s probably least confusing to just replace Q in your equation (6) with Q + f_s*R_s where f_s is given by equation ? (currently in supplemental materials). 

Lines 215-216: This sentence can be improved. 

Note that the limiting case of low wind was noticed before Artale; e.g., Fairall et al. (1996) had already corrected for this.  I think it is helpful to provide all the equations used in the code (e.g., include the Artale expression for lambda), perhaps in the supplemental materials. 

I encourage the authors to thoroughly check the equations to make sure they are stated correctly and implemented in the code(s) as intended.  For example, I believe While et al. (2017) has a typo in their thermal skin layer description.  They also use the Artale et al. (2002) scheme, but their notation (equation 1) introduces a constant beta, but the value they state (beta=864) does not match the value that would be derived from the Artale et al. (2002) (equations A1 and 3) or Tu & Tsuang (2005) (equations 1 and 6) descriptions!

Line 220: No need for the line break

Line 250: Why does a1+a2+a3 not equal 1?  Why is the scheme designed like that? Is it designed differently to equation (13), where the a_i’s sum to 1 (as expected)?

Line 259-261 (also Table 1 caption): Can you be more careful here.  Paulson & Simpson (1981) propose/use the 9-band parameterization and take the coefficients from Defant (1961) based on clear water data.  Soloviev and Schlussel (1996) provide a variety of values for b_1 (your notation) based on different water types (from Jerlov, 1976).  However, without knowing what the Jerlov water type is this is no use. Instead, what has often been done is to take the mean value of b_1 for Jerlov types I, IA, IB, II, and III.  Hence, b_1 becomes 0.1488.

Note, that Gentemann et al. (2009) actually keeps the original Paulson & Simpson value for b_1 (i.e., pure-water), but includes solar angle in the parameterization (see their equation (14)) - is that a helpful modification and is it implemented today!? If not then the reference to Gentemann et al. (2009) is probably not helpful.

Line 263: Is it right to call this ratio a “coefficient”?

Line 266: “… depth at which transmission drops by 1/e from its surface value …” don’t you mean “to 1/e” and not “by 1/e”!?, i.e., the depth it reaches when it decreases *to* ~37% of the surface value. (see also line 285 and line 303).

So the e-folding depth is 3m in the Soloviev (1982) scheme, but it is 1.6m in the Soloviev and Schlussel (1996) scheme – is that correct?  And it looks like using 3m in the S&S (1996) scheme takes you to 25% of the surface value (if I’m correctly reading and understanding your Figure 2a). 

The issue you are raising is not really whether d=3m is a good choice or not, but that the 9-band parameterization with fixed coefficients is not sufficiently accurate within the top 3 metres.  This is already well known.  Hence, we have solar transmission parameterizations whose coefficients/parameters change based on ocean color/chlorophyll data.  Key references include Morel & Antoine (1994), Ohlmann & Siegel (2000), and Lee et al. (2005).  These parameterizations have all previously been used and assessed in diurnal warming modelling studies.

Line 272-276: Please explicitly state this parameterization with formula and table of values as necessary.  Determining this e-folding depth seems to be your key update to the prognostic model, but you say very little about it.  Have other studies verified/assessed its accuracy? Can we have the full details, use the supplemental materials if necessary. 

If I understand correctly, you are using this approach to estimate the e-folding depth and then using this value of d in equation (10), basically no longer requiring the 9-band parameterization (i.e., do not use the right-hand side of equation (13)). 

You have shown that equation (13) is not very good at predicting the e-folding depth, i.e., the fixed transmission coefficients typically used are not very good at predicting solar absorption within the diurnal layer.  For the Mediterranean Sea at least, this leads in an underestimation of solar penetration (as you note on line 271 and we see in Figure 2), making the simulated diurnal layer too warm on average. (And I’m trying to link this understanding to your results/figures – so if the simulated diurnal layer is too thin and too warm it would limit the development of large diurnal cycles in SST, hence why the magnitude of the diurnal cycle is increased with the new scheme – am I on the right track? It would be helpful to spell out some interpretation/explanation in the results/summary sections to help the reader understand what is going on).

How does this approach compare to the schemes that use non-invariant transmission coefficients; e.g., the Ohlmann & Siegel (2000) and Lee et al. (2005) parameterizations?  Why not just use a state-of-the-art parameterization to begin with?

Section 3.3.1:  Can you specify in the text the temporal frequency in which you determine the e-folding depths, and how this average is obtained (is it the mean of daily values over the period 2019-2021!?).  Is the intention to use these seasonal e-folding depths for all future uses of the new scheme; i.e., you only need to compute this dataset once, a prior? Although obviously would need to be done again for a different region of the ocean.

Line 280: This is the first mention of Takaya et al., (2010) and assumes readers are already aware of what those refinements are.  Please briefly explain here or earlier the changes/improvements T10 made to ZB05. 

Also, just to confirm WRF uses ZB05 without T10 and A10 (your case 2: wrfskin), but otherwise is identical to the NEMO implementation of ZB05 (your case 3: nemoskwrite). Correct?

Section 3.4:  Table 1 is helpful and needed as the text description is not always so clear at first reading.

Line 299 (also Table 1): As mentioned earlier, I think the Gentemann et al. (2009) reference is potentially misleading, as I don’t believe the scheme is what G09 state/use in their paper.

Line 305: “… (see section 3) below.”  Do you mean “(see section 3.3.1)” which is actually above!

Line 314-315: suggest rewording “… simulations outputs with data … assess methods performances …”

Equation (14): Do you mean to start the index at i=0? If so, wouldn’t the divisor be N_seas + 1?  Also, I think it’s better/correct to write max SST – min SST (i.e., don’t pull SST(h_i) outside the bracket as you currently have).

Line 331: “… for every of the ...” --> “… for all of the …”

Line 342: “over estimation of the wrfskin … except winter” Isn’t the modrdnemo even more over estimated than wrfskin!?

Line 344: Your wording should be more careful in the first line of this paragraph.

Line 350-352: modradnemo looks similar, but slightly worst, to me compared to nemoskwrite (except for May)!?  Maybe a table of values would also complement the figure.

Section 4.4: Can you try and explain why you see these changes, why are the different methods influencing the results?  I’m trying to make sense of these results (figures 4 and 5) based on what I see in Figure 2. Please can you help the reader in this regard.

Line 359: “can be also” --> “can also be”

Supplemental Materials:

The first term 0.065 modifies Fairall et al. (1996)’s value of 0.137 (The Fairall value of 0.137 is obtained by summing a_5 to a_9 from your Table 2).  Wick et al. (2005) suggested using 0.067 (0.137-0.07) as an ad-hoc adjustment to provide results that were comparable to the better Ohlmann & Siegel (2000) parameterization.  However, I am unsure why the value of 0.065 is stated in Zeng & Beljaars (2005) (also in Zhang el al. (2021) so perhaps has become established) and repeated here in your supplemental materials.  Which value is used in the actual code?

This number can make a difference.  For example, suppose delta was 1mm, then f_s (equation (17)) in Fairall et al. (1996) would give 0.1009; i.e., 10% of solar radiation absorbed (or using the reduced value in Wick et al. (2005) it would be 0.0309, or using the value from Zeng & Beljaars (2005) 0.0289).  I suppose I don’t really get why this approximation is even needed in the first place.  Why not just plug delta into the full Paulson & Simpson (1981) 9-band parameterization?  You have the complete expression on line 202. (This is in fact is what Paulson&Simpson (1981) did see their equation (6)).  Wouldn’t a chlorophyll-based scheme help here as well?

The delta equation in your supplemental materials, although identical to equation (6) in Zeng & Beljaars (2005), is not the same as the referenced Fairall el al. (1996) equation (14) [note that this is an equation for lambda not delta!], I believe there must be a typo in Zeng & Beljaars (2005).  The correct Fairall et al (1996) equation is reproduced in Tu and Tsuang (2005) (see their equation (4)) and Zheng et al. (2021) (see their equation 4).  Please clarify.  And perhaps more importantly is the correct expression used in the actual code(s)! 

In this paper Authors propose an updated scheme to simulate heat propagation at the ocean-atmosphere interface. The importance of the task is related to heat budget estimates, assimilation of satellite data.

The study proposes an interesting and important upgrade, but in my opinion there are points that could be improved:

• Description of the model/new parameterization. It would be useful to add a diagram illustrating the layers, the cold skin layer and the warm layer. In the description of the model (e.g. Eq 5 and 10) it would be useful to understand what the prognostic variables are and how this parameterization is related to the OGCM internal variables. (e.g. Nemo potential temperature). Nemo vertical discretization in the upper layer could be also shown in the diagram.

• Authors should provide a better description for the choice in the variability of the terms in the new approach. Considering the formula for the solar transmission equation, the exponents b_i should also be space/time dependent because light attenuation should be modulated by chlorophyll (and other optical constituents). This is suggested in Table 1 (R-G-B + chl e-folding) but this should be reported with a clear formula. This variability of attenuation coefficients b_i  is directly connected to the change in d. In particular authors imply that transmission is not constant due the chlorophyll variability and then  they assume “constant transmission throughout the basin, but with a spatially and temporally varying e folding depth and defines our new prognostic scheme for skin SST warm layer calculation”. I think this is a central part of the work and should be better illustrated and explained.

Below I report some technical points.

P2 Line 64 “Simplified approaches” rather then “Simplified models”

P4 Line 114 “regions where the percentage of model data is higher than 50% have been masked out both in 115 CMEMS MED DOISST and our experiments”

Not clear “higher than 50%” higher in respect to what?

P6 Lines 183 and 185 please put units of measure of quantities considered.

P6 Line 191 “assuming this constant” which constant, maybe the const in eq 2? 

Please clarify.

P7 Lines 212,214 please put units of measure of the symbols.

P7 Line 230 Eq 8 is not clear. Please write better the argument of \phi function.

P8 Line 241 “Assuming a temperature of dependence”

Please explain better.

P8 Lines 245,246  I would specify if the equations are coupled e.g  “In ZB05 scheme (Zeng and Beljaars, 2005), eqs. (5, 10) are  the coupled equations for the cool skin and warm layer respectively.” As already reported above, it would be useful to make clear which prognostic variables are taken into account.

P8 Line 246 “ within this layer” which layer? Maybe “these layers”?

P12 Line 337 “Looking at the mean profile averaged over all grid points in the given area, the agreement is better for all simulations during summertime months, both for the eastern and the western region (see figs.  8c, 9c), showing in particular that the modradnemo simulation outperforms the nemoskwrite one”

I don’t see where modradnemo simulation outperforms the nemoskwrite one, they seem equivalent in terms of skill, could authors be more quantitative? Some tabulated skill metric could be useful.

P12 line 376 “On the other hand, in the western Mediterranean all simulations tend to overestimate the signal, with our modified scheme doing a better job. ”

 Could authors be more quantitative, adding some statistics to support their statements? 

Figure 2 panel (a) the plot could be improved. I suggest putting the z axis vertically (with negative ticklabels for depth as used in equations).

Figure 10. Region 0 should be removed, it’s not a sea region.