Response to Reviewer #1 comments.

First of all we would like to thank the Reviewers for taking the time to read our article and to prepare a thorough review. The corrections and improvements which we are going to introduce to the revised version of the manuscript will certainly increase the clarity and quality of the paper. Please note, that, at this stage according to ESD rules of submission, we do not provide the revised version of the manuscript and the intended changes are appropriately discussed.

Regarding the wave generation:

It is not clear why the authors are using wavemaker theory here. In a numerical model, more elegant methods are available to intruduce waves. While "monochromatic wavemaker motion" gives a clear number of cases to be examined, a more realistic scenario would be to consider a wave spectrum

The issue raised in this comment corresponds to the second comment of Reviewer #2. This is why, we decided to provide the answer, which is repeated in the reply to Reviewer #2.

The paper presents the methods of analysis of the problem of modelling of wave-induced vertical mixing from the wavemaker perspective to be verified in the wave flume. We focus on the side effects of the waves generated in a hydraulic laboratory using the presented mathematical model to increase awareness of the influence of these effects on the experimental outcome and correct interpretation of the results. This problem was already mentioned in our previous studies on wave-induced mixing (please see, Sulisz and Paprota, 2015, Math. Probl. Eng.). We have discussed this problem in closer detail in section 3.4. We will also add this explanation in the last paragraph of introduction for more clarity in the revised version of the manuscript.

We start our study from simple cases to be able to investigate the influence of laboratory side effects for regular waves. We then are able to separate the Stokes drift from Lagrangian and Eulerian mean velocities using the method presented in our previous works (Paprota et al. 2016, J. Hydraul. Res., 54, 423–434). Since Stokes drift is derived based on weakly-nonlinear assumptions, our higher-order approximation gives improved estimate of the drift for steeper waves (relevant to more severe conditions). Simple models basing on Stokes drift may be also applied to random waves characterised by the whole spectrum as presented by e. g. Myrhaug et al. (Myrhaug et al. 2018, Oceanologia, 60, 305-311, https://doi.org/10.1016/j.oceano.2017.12.004), which is more relevant to ocean waves. We will add this explanation to the revised manuscript with the above-mentioned reference in section 3.4.

About the particle tracking:

The authors state that  “the improvements to the method of evaluation of mass transport velocity based on the Lagrangian particle tracking (Paprota and Sulisz, 2018) are introduced” (line 168-169). It is not quite clear what they mean by this. Please be more specific about the improvements. 

In the last paragraph of section 2.1 we state that:

“In the present study, the improvements to the method of evaluation of mass transport velocity based on the Lagrangian particle tracking (Paprota and Sulisz, 2018) are introduced. In order to get better estimation of the time-independent velocity field, the two hydrodynamic states corresponding to both zero up- and down-crossings of the regular wave are used to start-up the tracking procedure - contrary to the previous method based only on the zero down-crossing initial position (e.g. Paprota et al., 2016; Paprota and Sulisz, 2018).”

In order to make this paragraph more clear we are going to rephrase it. First, we will provide more details on the procedure from our previous works (Paprota et al. 2016, J. Hydraul. Res., 54, 423–434. Paprota and Sulisz, 2018, Phys. Fluids, 30, 102 101). It was based on tracking particles initially evenly distributed along the depth for zero-up crossing phase of regular waves. Now, we use both down and up crossing phases, which gives more accurate estimation on resultant mean velocity field. We will add this description to the revised manuscript.

It is also not clear why we need Lagrangian particle tracking siden for example eq. (29) only uses the Eulerian velocities.

We need the Lagrangian particle tracking to obtain mean transport velocity field, which in our case is a sum of Stokes drift and return current. The inclusion of Stokes drift (without return current) to advection equation was already discussed explicitly for the case of ocean waves by (McWilliams and Sullivan, 2000, Spill Sci. Techn. Bull., 6(3/4), 225-237, https://doi.org/10.1016/S1353-2561(01)00041-X) and is also mentioned in (van den Bremer and Breivik, 2018, Proc. R. Soc. A, 376, 2111, https://doi.org/10.1098/rsta.2017.0104). We will address this issue with relevant references at the beginning of section 2.2 in the revised manuscript.

About the mixing: What is the relative size of the mixing efficiencies kappa(m) and kappa(v) ? The authors state the the dimensionaless parameter alpha has been measured ..... Can you be more specific? 

More information on the values of \kappa_v can be found in the previous work (Sulisz and Paprota, 2015, Math. Probl. Eng.), where additional analysis is presented. In order to be more informative, we will add one more column to Table 1 in the revised version, which will present the maximum of \kappa_v in relation to \kappa_m for each of wave cases. For example for the first case (kh = 0.5, Ak = 0.0125) max(\kappa_v) is approx. 16 times lower than \kappa_m, while for the last case (kh = 2.0, Ak = 0.2) max(\kappa_v) is approx. 2000 times larger.

Also the procedure of estimating (not exactly measuring) of parameter \alpha is presented in the work (Sulisz and Paprota, 2015, Math. Probl. Eng.). Its value based on measurements from the work (Sulisz and Paprota, 2015, Math. Probl. Eng.) is equal to 0.002. We will present this explanation together with relevant referencing and \alpha value in the revised manuscript.

In the abstract, the authors state that this work may lead "to improved estimates of subsurface mixing intensity and ocean surface temperature." Do the authors mean in the nearshor ocean. Note: in the caption of Table 1, the authors state that they "wave-induced vertical mixing processes in offshore conditions." However, the parameters given in this table appear to be mostly relevant for surfzone Dynamics. Two questions: 1. Does this study apply to nearshore or offshore or both? Please address this issue. 2. There could be other effects of equal or greater importance on mixing (either shallow or deep water). Please address this issue.

The model may be applied to deep-water as well as shallow water conditions, which was already presented in (Sulisz and Paprota, 2015, Math. Probl. Eng.). Here, we present the cases corresponding to kh = 0.5 … 2.0 for two reasons. First of all we want to directly compare the cases presented in the previous work (Sulisz and Paprota, 2019, Ocean Engineering, 194, 106622), in which only weakly-nonlinear model was presented to showcase differences and importance of including higher-order terms (strong nonlinearity). The second issue is that in an analysed region, for deep-water cases corresponding to kh = \pi, the mixing would completely swept away the hotter water from the domain. We will address this issue as presented here in the revised manuscript.

The arrows in Figure 2, 3 and 4 are not extremely informative. Please add more explanations of what can be seen and learned from these figures.

The arrows represent the vectors of a phase-averaged velocity corresponding to mass-transport induced by waves. Black arrows correspond to weakly-nonlinear solution, while green and blue arrows correspond to higher-order solution and two methods of averaging – EMTV and LMTV, respectively. This explanation will be added to the second paragraph of section 3.2 for clarity in the revised manuscript.



In section 3.4, the authors mention "laboratory experiments" and "wave flume" many times, giving the impression this study was conducted in a laboratory. Together with the introduction (for example third paragraph) this leads to confusing the reader, and deflects the focus from the numerical work. Please be very clear what this article covers, and what it does not cover, both in introduction and in discusssion. For example, line 46: "First, the problem of the generation of waves in a laboratory flume is formulated and solved" This is misleading.

We will modify the text with this respect. Namely, we will change laboratory flume to numerical one, wherever it is possible (as for example in line 46).

Reply to Reviewer #2 comments

First of all we would like to thank the Reviewers for taking the time to read our article and to prepare a thorough review. The corrections and improvements which we are going to introduce to the revised version of the manuscript will certainly increase the clarity and quality of the paper. Please note, that according to ESD rules of submission, at this stage we do not provide the revised version of the manuscript and the intended changes are appropriately discussed.

1. Nonlinear vs weakly nonlinear: Is this problem truly nonlinear? It seems to me that a major feature of nonlinear waves, namely wave breaking, could not be captured in this model. I suspect this case would correspond to where the expansions for phi ceased to converge. While wave-breaking is not relevant here, I feel that the requirement of convergence means the solution given is just a higher order weakly nonlinear model (as stated) and as such, the use of the term 'strong nonlinearity' in the title is misleading. Can you discuss this point and convince the reader one way or another whether this methods fully captures nonlinearity.

The problem is truly nonlinear, while we solve Laplace equation with nonlinear boundary conditions more accurately and beyond the applicability of weakly-nonlinear approaches. Our higher-order method allows modelling of waves with strong nonlinearity indicated by a high value of Ursell parameter (see Table 1) and admits amplitude dispersion, nonlinear wave-wave interactions in deep and intermediate waters as well as solitary waves propagation. Corresponding methods basing on pseudo-spectral approach also refer to nonlinear or fully nonlinear waves although they consider only non-breaking waves (see e.g. Paprota and Sulisz , 2019, J. Hydro-Eniron. Res., 22, 38-49 for review). We will add this reference with short explanation in the last paragraph of introduction.

2. Problem setup and oceanographic relevance: The problem is pitched as being of relevance to oceanographers though the setup modelled is a labatory one. While there are results which may be relevant to the ocean, these are mixed with results which are not. For example, the flow around the wavemaker paddle is unlikely to be of interest to anyone who isn't an experimentalist. It feels as if the authors plan to compare with experiments at a later date. I would suggest reframing the paper as an engineering problem with oceanographic relevance rather than the other way around. Results could be discussed in terms of the problem studied and a new section could be included which transfers the relvant results to an oceanographic context.

The issue raised in this comment corresponds to the first comment of Reviewer #1. This is why, we decided to provide the answer, which is repeated in the reply to Reviewer #1.

The paper presents the methods of analysis of the problem of modelling of wave-induced vertical mixing from the wavemaker perspective to be verified in the wave flume. We focus on the side effects of the waves generated in a hydraulic laboratory using the presented mathematical model to increase awareness of the influence of these effects on the experimental outcome and correct interpretation of the results. This problem was already mentioned in our previous studies on wave-induced mixing (please see, Sulisz and Paprota, 2015, Math. Probl. Eng.). We have discussed this problem in closer detail in section 3.4. We will also add this explanation in the last paragraph of introduction for more clarity in the revised version of the manuscript.

We start our study from simple cases to be able to investigate the influence of laboratory side effects for regular waves. We then are able to separate the Stokes drift from Lagrangian and Eulerian mean velocities using the method presented in our previous works (Paprota et al. 2016, J. Hydraul. Res., 54, 423–434). Since Stokes drift is derived based on weakly-nonlinear assumptions, our higher-order approximation gives improved estimate of the drift for steeper waves (relevant to more severe conditions). Simple models basing on Stokes drift may be also applied to random waves characterised by the whole spectrum as presented by e. g. Myrhaug et al. (Myrhaug et al. 2018, Oceanologia, 60, 305-311, https://doi.org/10.1016/j.oceano.2017.12.004), which is more relevant to ocean waves. We will add this explanation to the revised manuscript with the above-mentioned reference in section 3.4.

3. Discussion of new material: this work is building on previous studies by the same authors and it is not entirely clear what is new. In some places old results are repeated and in others technical details are skipped and it is not immediately clear if they're covered elsewhere. I think the authors should clarify their new contributions and give a clearer exposition of their previous work, either leaning entirely on another reference or repeating enough (clearly labelled) content that the paper stand alone.

In the present paper, the introduction of higher nonlinearities allowed for describing the enhanced sub-surface streaming when compared to previously published studies relying on weakly-nonlinear theory. The second important contribution is presentation of two methods of averaging of the wave velocity field either basing on Eulerian averaging or Lagrangian particle tracking, to demonstrate better accuracy of the latter especially in immediate proximity of the wavemaker. The third novelty is the improving of the Lagrangian method of averaging which will be more thoroughly discussed in the revised manuscript. We are going to add the relevant (above mentioned) information on differences to previous studies especially (Sulisz and Paprota, 2019, Ocean Eng.) and new contribution to introduction section last paragraph and conclusions to the revised version of the manuscript.

I have made various comments on the attached PDF document.

Thank you for your additional comments improving the quality of the presentation.

• Revise grammar

The first sentence of the abstract will be revised.

• Why? There's various setups you could choose which are more relevant to the ocean. Are you planning to directly compare results with experiments?

This is already clarified in query no 2. We will change the word “laboratory” to “numerical”

• Clearer to write \phi_{n+1} ?

Thank you for pointing this out. We will work on the precise formatting in the revised version.

• Is this approach fully non-linear? or just a higher order weakly non-linear theory?

As discussed previously, we use another approach which relies on a spectral method and is different to weakly-nonlinear solution basing on perturbation expansions with respect to small steepness parameter discussed in section 2.1.1. We will address this issue as previously stated to the revised manuscript.

• What is this procedure? Are there expressions for A_i? If this solution has been given in previous work, what is new?

The wavemaker solution was already reported in (Paprota and Sulisz, 2018, Phys. Fluids). We will refer to that in the revised text. We will also provide explicit formulas for A_i and a_i, which are calculated as coefficients of the Fourier expansion in the revised manuscript.

• What theory is this based on? Weakly nonlinear, fully nonlinear etc? When is it valid and do you expect it to align with your higher order theory presented above?

The derivation of \kappa_v is based on the weakly nonlinear theory. More general equation is provided in (Sulisz and Paprota, 2015, Math. Probl. Eng.), where there is a comparison between weakly-nonlinear and more general form. The parameter \alpha was estimated based on the experiments only for the presented form of \kappa_v and this is the form used in the present study. In order to use more general formula, we need to determine \alpha using experimental data from a greater range of wave conditions. We will add this explanation to the revised version of the manuscript just after the \kappa_v formula.

• Which measurements? What is a typical value?

Again, the \alpha is determined based on experimental data. This will be clarified in the revised manuscript in this particular line/paragraph.

• I would suggest splitting this section into a discussion of the behaviour near the wavemaker and far away. The far field behaviour is presumably more relevant to the ocean?

We have addressed this issue in query no 2. If you do not mind, we would like to keep the structure, but instead add a relevant information on near-field and far-field kinematics at the beginning of section 3.2.

• These differences are not immediately clear from the arrows on the figures given. The three cases (WNL vs EMTV vs LMTV) look qualitatively very similar. Can you plot a more quantitative diagnostic that makes these differences more apparent?

 Indeed, this point needs additional clarification. We will precisely refer to particular location of the velocity field in the text to clarify this. The differences are highest near the wavemaker.

• So what is new?

The differences will be indicated in this part as well, thank you for pointing this out.

• Originating

This will be corrected according to your suggestion. Thanks.

• Is this return current likely to be relevant for the ocean?

Since the Stokes drift moves water in the direction of wave propagation it must be somehow compensated by the opposing flow. In the flume it is well recognized phenomenon called return current, but in open ocean conditions this process is not fully understood and is an open question due to complexity of the problem of ocean dynamics.