I agree with the main premise & conclusion of this manuscript that Stoke's drift, very near surface currents via coupled modeling are critical to dispersion of surface material (such as oil spills, and microplastics). The manuscript is solid, well written, includes insightful test cases. The authors are informed scientists in this field.

The main deficiency of the manuscript is that it is modeling only and lacks real-data verification... which is very very hard, so understandable.

Also, some relevant work is not cited:

Curcic, M., S.S. Chen and T. M. Özgökmen, 2016: Hurricane-induced ocean surface transport and dispersion in the Gulf of Mexico. Geophys. Res. Lett., 43, 2773-2781.

Laxague, N.J.M., T.M. Özgökmen, B. K. Haus, G. Novelli, A. Shcherbina, P. Sutherland, C. Guigand, B. Lund, S. Mehta, M. Alday, and J. Molemaker, 2018: Observations of near-surface current shear help describe oceanic oil and plastic transport. GRL, 45(1), 245-249.

The first paper involves testing with a truly coupled ocean-wave-atmosphere model with using real data and the second involves current measurements going all the way to top 1 cm near the air-sea surface (which was very hard to pull off); could be interest to the authors. (I had to declare my identity since I am posting these papers.)

Review of “Non-negligible impact of Stokes drift and wave-driven Eulerian currents on simulated surface particle dispersal in the Mediterranean Sea” by S. Rühs et al.

This manuscript analyzes the impact of including ocean surface gravity waves when simulating ocean currents that are used for Lagrangian advection.  The conventional approach is to consider the effect of surface waves on Lagrangian particles through their phase-averaged Stokes drift.  Previous studies have therefore considered adding the Stokes drift to the output of ocean models that did not account for wave driven processes.  This study employs an ocean circulation model that is coupled to an ocean surface wave model, and thus considers the impacts of surface waves on the ocean circulation model physics, including through the impact of waves on vertical mixing, surface fluxes, and on the resolved scale currents/advection in the ocean model.  The results show that Lagrangian parcel simulations using the ocean circulation model with full wave coupling yields a different result from simulations using the ocean model without wave coupling.  This difference can not be accounted for by adding the Stokes drift to the output of the non-coupled model.  The study therefore suggests that the feedbacks of ocean waves to modify the background current should be properly accounted for when using output of numerical ocean circulation models to drive Lagrangian particle simulations.

I found this study to be well-written, thoughtful, and important for the ocean modeling community.  I have a few comments on the presentation that I think can better clarify the result and its place within the scope of similar literature.  

General Comments

1. A general comment for definitions of the decomposition in the equation at L50.  The Lagrangian current is routinely separated into an Eulerian component and a Stokes drift component that is attributed to the surface wave field.  In this study the Eulerian current is further separated into a “non-wave” component U_Enw, and a wave-driven component U_Ew.  The non-wave component is later defined from the non-coupled simulation and the wave component is defined from the residual of the coupled simulation minus the non-coupled simulation.  This is useful conceptually to decompose the Eulerian current and explain the results.  The approach here is pragmatic, but I do think it is important to make the choice of this definition clear early in the paper (e.g., as is conveyed later in Table 2).  There are non-linear terms when you back substitute U_Ew and U_Enw into the momentum equations, such that one could consider more refined ways to decompose the wave-driven part, not just from this residual approach.
2. It is important to clarify that for the Wagner et al. (2022, also disclosing that I am a coauthor on that work) work we conducted our experiments on a 25 km ocean-wave model and resolved wave-current interactions at that scale.  There are other differences between this work and our simulations (e.g., more wave physics impacts than just resolved-scale wave-current interactions are considered in this work), but I personally did not find it surprising or controversial that a different result may be found at 4km resolution with wave-current interactions resolved at smaller scales.   The results may in fact be compatible, perhaps yielding insight into the scales where the resolved scale impacts are important.

Specific Comments

L53:  I suggest using different language, "explicitly resolved" to me implies simulating the surface phases of waves directly by the ocean model.  But I think it is meant that they are sometimes coupled to spectral wave models.

L115:  I don’t think Craig and Banner (1994) is the best reference for Langmuir turbulence.  Perhaps McWilliams et al. (1997, doi:10.1017/S0022112096004375) and other more recent reviews (e.g., Belcher et al., 2012 doi:10.1029/2012GL052932, D’Asaro, 2014 doi:10.1146/annurev-marine-010213-135138)?

L148:  But it does neglect the feedbacks of U_Ew on U_Enw.  This is a benefit of the approach here, could that be tested here?

L152:  “controversy” seems like an overly strong word to me.

L166:  WaveWatch -> WAVEWATCH.

L169: despite -> except

L175:  What is the first cell thickness?  I can imagine the results could be sensitive if the first cell is particularly thin or coarse.

L200:  Cell horizontal interfaces or vertical interfaces or both?

L220:  Are there any citations for this?  Which specific “TKE” scheme is it?  A k-l type, a Mellor-Yamada, GLS, etc.?  Is the momentum flux directed only down the Eulerian vertical current shear or also down the Stokes gradient (e.g., Harcourt 2013, doi: 10.1175/JPO-D-12-0105.1)?  This detail could be important since down-Stokes mixing can be an additional source of “anti-Stokes” current.

L225:  Is there a citation for this?  Otherwise, it may be better to say the parameterized Langmuir turbulence is expected to be more realistic with the simulated, sea-state dependent Stokes drift than the wind speed based Stokes drift.

L3.1.2:  Some discussion of the wave model performance in this configuration would be helpful.  Are there obs comparisons in a previous study that can be cited?

L229:  Is there a spectral tail added for the Stokes drift computation?

Table 1:  It would be inconsistent to include some of these Stokes drift impacts in NEMO without others, so I suggest not splitting into three subcolumns when intermediate experiments aren’t attempted.  I didn’t catch the distinction between the modified TKE scheme [note typo in manuscript] and Langmuir turbulence parameterization, if more details are given as noted by comment at L220 it would help here.

L244:  Why not spin-up the coupled version from rest?  Is it possible that analysis in early 2019 is contaminated from the initial adjustment?

L252:  I don’t expect these results to be particularly sensitive to this detail, but I’m surprised that the atmospheric fields are only updated every six hours.  This seems fairly coarse in time at ~10km spatial resolution.  Are the fields interpolated in time to force NEMO and WAVEWATCH?

L266:  Suggest to clarify if Stokes drift is similarly averaged over 1m.

L295:  Specify horizontal grid, vertical grid, or both

L302:  This seems like a missed opportunity in this study, otherwise it leaves an open question if ocean circulation models need to include full wave physics or the effect can still be partially accounted for via intermediate approaches.  It seems very relevant to me to answer the second question.  Can the authors offer some comments on its potential utility in Section 5.1?

Table 2:  I find this table quite useful interpreting the definitions, it would be useful to refer to this table when defining the u_Enw and u_Ew components.

L333:  This is a practical approach, but I think this is an important point to make earlier (e.g., when discussing the decomposition).  It is important to know that it is defined as a residual and includes all the feedback that would be missed in the intermediate approach.

Figure 4:  I find the bar plot (panel a) busy and difficult to understand.  You may consider if the maps are sufficient on their own to make less effort for a reader to understand the figure.

Figure 6:  The differences between the panels are often subtle. I wonder if showing the difference from the 1st experiment instead of the value for the 2nd and 3rd experiment would make a clearer indication of the changes?

L609:  Suggest to clarify what is meant by intrinsic variability in this context.  It is an eddying model, so I did wonder if 2 years is sufficient experiment length for all the statistics?

L625:  This may be true and is a worthwhile point, but I don't know that this has really been tested in this work.  This study shows differences in the Eulerian (gridded mean) fields between coupled and non-coupled, which weren’t found in Wagner et al. (2022).  As mentioned in the general comments there are other differences between these works, particularly the horizontal grid-spacing, which could explain the different conclusions.

Review by: Brandon Reichl