Summary

The authors used ANSYS Fluent to simulate a Reynolds-Averaged Navier Stokes (RANS) model of seawater intrusions upstream of ice-shelf grounding lines in the presence of subglacial freshwater discharge. Running different simulations, they found that freshwater discharge velocity and conduit height are the parameters that have the strongest influence on intrusion distance. The former (resp. latter) decreases (resp. increases) intrusion distance, in qualitative agreement (but quantitative disagreement) with recent theoretical predictions. A key result of this study is that the spatial structure of the intrusion depends significantly on how the ice-water interface is modelled. A comparison of the simulation results with state-of-the-art parameterizations for melt rates is presented. The computed melt rates and parameterized predictions are quantitatively different (order(s) of magnitude different).

General comments

The authors aim to tackle a fascinating problem of practical importance. Melting near grounding lines is thought to have much more impact on glacier dynamics than melting anywhere else, yet observations are limited, and high-fidelity simulations are lacking. Thus, the work is highly novel and has significant potential for improving our understanding of ice-ocean interactions that have a high impact on climate dynamics. However, I am not sure that the ANSYS Fluent RANS solver used is appropriate for this problem. In fact, I have never seen the ANSYS Fluent RANS solver applied to environmental flows as complex as in this work. This does not mean that it cannot work, but that significant effort should be devoted to validating the code. As the governing equations solved by the codes are not clearly presented or discussed (in particular, the boundary conditions, which are yet key to evaluating the melting dynamics), I have found it difficult to assess the appropriateness of the mathematical/numerical model. I am particularly concerned with the modelling of the ice-ocean boundary: there are no salt constraints and melt is said to be activated in a way that I did not find correctly physically motivated. In the unfortunate case that the code cannot in fact solve the exact problem at hand (with temperature and salt stratification coupled through a phase change boundary) I would suggest that the authors reformulate the problem as a list of hypotheses--inspired from the full problem--and testable with simulations of a modified/simpler model (but mathematically transparant), which could be the reduced thermal driving model discussed below (point 2).

List of specific important comments

1. The governing equations that the code solves should be clearly presented, as I do not think that TC readers are familiar with RANS models (and because I have never seen RANS models used in this context). This requires presenting the Reynolds decomposition of the flow (which would help you justify the fact that 2D dynamics is expected since turbulence is not resolved) and the governing equations for the ensemble-average variables. The closure for the Reynolds stresses and turbulent fluxes should then be discussed in greater details than in the current manuscript. The kappa-eps scheme is mentioned (Eq. (2)) but important details are lacking: for instance, how are kappa and epsilon related to the resolved variables? I would like also to see more details on the so-called damping functions of the low-Re formulation that supposedly enable accurate diffusive boundary layer representation. Is it like in a wall-resolved large-eddy simulation model?
2. The decoupling of the salt dynamics from the ice-ocean boundary dynamics is not justified and a priori seems wrong. I assume that this decoupling is due to code limitations. However, if you cannot justify the decoupling, I am afraid this just means that the code is not suited for environmental flows with temperature and salt stratification and melting. In the worst case scenario you might consider reformulating the problem in terms of a single scalar variable, namely thermal driving (ref: Adrian Jenkins' papers and other people's related works). The collapse of the full dynamics onto a reduced thermal driving model is thought to be accurate in highly-turbulent environments (for which kappa-epsilon applies anyway) and small salinity variations. This latter condition is obviously problematic (which should be discussed) with regards to your problem of interest. However, I would rather see thermal driving model simulations transparantly solved by ANSYS than results from a non-transparant full temperature-salinity model.
3. The boundary conditions, especially at the ice-ocean interface, should be clearly presented and discussed. The physical motivation for the so-called melt-activated formulation is lacking. Melting produces buoyant flows even when the boundary is no slip, simply because melting acts like a sink for salinity (though this is lacking in your model). Thus, should we envision your melt-activated formulation like a velocity compensation for the lack of salt sink in the model? If so, it was not clear to me whether the velocity is prescribed vertically or horizontally, and I am not sure it is the correct way to compensate the salt sink. Estimating the velocity to enforce at the boundary to mimick the salinity sink from Eq. (4) is also not justified, i.e. why should the movement of the interface (assuming there is no immediate hydrostatic equilibrium of the ice shelf at such small scales) be directly used as a buoyancy-driven velocity input?
4. In would like to see a validation of the code, which includes the choice of turbulence closure. A simple benchmark case should be set-up, for which turbulence-resolving simulation data exist (either from direct numerical simulation (DNS) or large-eddy simulation (LES)). Several groups (with Catherine Vreugdenhil, John Taylor, Ken Zhao etc) have published such DNS and/or LES data over the past 5 years or so, making it practical. Typical configurations are channel flow configurations, which should be accessible to ANSYS. Validation could be based on quantitative comparisons of mean variable profiles (e.g. temperature, TKE) normal to the ice-ocean interface.
5. Successfull applications of ANSYS Fluent RANS solver to environmental flow configurations with temperature and salinity stratification should be cited and discussed (in particular how they validated the code).
6. The discussion of the steady state and transient runs is really confusing. You should distinguish the existence (or non-existence) of a steady state from your strategy of successive runs to achieve it. The key point that should be in the main text is that the problem has a natural steady state (for the ensemble-average variables) as there is no external variability and the ensemble-average variables do not exhibit temporal fluctuations once equilibrated. The strategy to reach it should then be discussed in an appendix.
7. I have found many typos in the appendices ("Need to list reference values, materials info, the methods and controls" found line 566), suggesting that these were not carefull reviewed by the authors. At the moment the appendices seem primarily like a draft list of comments and figures that did not make it into the main text (some in fact repeating what is already in the main text).
8. The figures are lisible but could be improved (e.g. the x-axis label of Fig. 6 is quite small).
9. Because of the many questions I had/have with respect to the simulation code, I was not able to appreciate the comparison of the simulation results with the parameterized predictions of melt rates. If the authors can validate their code, I agree that this comparison would be an important addition to the paper, but it would have to be a fair comparison. That is, if the authors end up solving a model that is distinct from the model that the parameterizations (necessarily approximate) aim to mimick (arguably the real exact model), they should discuss result differences in light of model differences.
10. Simulation snapshots and movies would really help visualize the flow.

Additional comments, including technical corrections: typing errors, etc.

1. line 29: could you clarify the idea of "tidally asymmetric"?
2. line 87: writing "The ice wall boundaries have a temperature boundary condition of 0◦C" really felt like a bomb! And the lack of justification or description of the full mathematical model did not help disarm it. It is really not expected that prescribing 0◦C at the ice-water interface is reasonable, especially when the salinity goes from 0 to 30 psu.
3. line 108: this equation of state may be suited for small salinity variations, but with a range of 0 to 30 psu, it may not be appropriate, unless the grounding line is beneath a lot of ice. I would recommend that you use a higher-order approximation of the true equation of state, or at least acknowledge that you know of the anomaly of the equation of state at low salinity and pressure but that you discard its effects to simplify the problem. Ideally you could discuss/speculate on how the results might change should you consider an accurate equation of state.
4. line 113: conduction, diffusion and molecular dissipation all sound like the same thing to me. Could you explain how they differ? (reading your appendices it looks like conduction might be turbulent conduction)
5. line 124: the quadratic quantity indicated is one among many, such that the sentence does not read well. Please reformulate.
6. line 126: could you recall what the Boussinesq hypothesis is?
7. Eq. (2) P2: what are the equations for kappa and epsilon? I expect that they involve the resolved ensemble-average variables.
8. Section 2.4: this section is really confusing as I already mentioned earlier, because buoyancy isn't due to some interfacial velocities but to changes in salinity and temperature at the phase-change boundary. Please reformulate.
9. line 150: I don't see the physical justification for prescribing a "horizontal ice (?!) velocity" a the horizontal ice-water interface.
10. line 179: can you provide a reference for "realistic estuarine-like mixing rates"?
11. line 193: rewrite "retrograde bed slope".
12. You refer to figure 5 before figure 4 so the two should be swapped.
13. line 253: Could we say "vertical baroclinic convective motion"?
14. Eq. (9) line 325: this equation doesn't look like a parameterization but rather like the exact expression for the melt rate as a function of the conductive heat flux at the interface.
15. Fig. 5-7: it was not clear to me whether the results you plotted were for the melt-enabled model or not.
16. Eq. (A3): should it be Re_L? Or change L into x?
17. Appendix A4: I have found many typos, you need to use capital letters for Reynolds number, kappa has become k etc Eq (A7) has signs/symbols displaced. Please discuss Eq (A7) more carefully. What is tau_eff? What is species diffusion in your case?

Review Mamer et al “High-Fidelity Modeling of Turbulent Mixing and Basal Melting in Seawater Intrusion Under Grounded Ice”

Summary:

This paper studies the estuarine-type dynamics of seawater intrusion/subglacial discharge upstream of an ice sheet grounding line in a CFD model (ANSYS fluent). The subglacial environment is two-dimensional with a large aspect ratio (5cm high and 3-5m long) and is forced by fresh subglacial discharge upstream and salty inflow from the ocean downstream. The authors investigate the structure and distance of the saltwater intrusion, and its sensitivity to various parameters (primarily freshwater discharge velocity). They find that that the height of the subglacial environment and the velocity of the freshwater discharge are strong controls on the intrusion distance, and the strength of turbulent mixing has a somewhat weaker/more ambiguous effect. The authors also investigate basal melt at the upper surface of the subglacial region and how it affects the intrusion (to do so they must change the boundary condition at the top surface away from the no-slip condition), finding that melt decreases the intrusion distance, offering a possible negative feedback.

This is an interesting study of an important process; however, I have several major concerns with the paper that may affect the conclusions. If they are addressed, I would be happy to review a revised manuscript.

Major Comments:

1. I am concerned about the turbulence modelling and the lack of inclusion of turbulence quantities in the paper. The fidelity (reported to be high!) of these simulations relies on the appropriateness of the turbulence closure, however, the reader is not given any context for the choice/s of C_mu, nor are we given any more intuitive quantities (e.g. the resulting turbulent diffusivity) to better understand the effect of varying C_mu, and how increased turbulence affects the flow. Without this information and context, I have low confidence in the modelling overall, especially since the effect of varying C_mu is non-monotonic for some cases (Fig 2B) and counter to my expectations for others (Figs 2A & C), i.e. I would have expected that increased mixing would decrease intrusion distance, whereas the authors find the opposite result. My recommendations to address this are as follows:
   • Present the turbulent diffusivity (or viscosity) alongside T, S, u. In main MS or in a supplement
   • Present the Reynolds numbers for the cases
   • Give some context for the choice of C_mu for similar flows – I presume the engineering literature can provide. Stratified plane couette flow may be a starting point in terms of a bounded, stratified flow.
2. I don’t understand why the temperature profiles have much sharper gradients than the salinity gradients (c.f. the darkest line in Fig A8A to the darkest line in Fig A8B). The authors mention the steep T gradients (not seen in S) in line 202 (Such a steep thermocline is most likely due to the temperature boundary condition we imposed on the horizontal ice boundary) however, I am concerned that the problem runs deeper than different boundary conditions. Why, for example, is temperature at 2cm depth at the GL entrance at 0.5 degrees (i.e. unmodified from ambient conditions), while salinity is <20ppt at the same location? There may be a serious issue with the mixing of these scalars which would significantly affect your results. Temperature-salinity plots may help determine if there is a problem.
3. By comparing the temperature and salinity profiles as u_f is increased (Figs 3 A,B, Figs A8 A,B,D,E), it appears to me that the greatest control on mixing is time spent in the subglacial channel. For example, at u_f=0.05 cm/s salinity is quite well-mixed over the full depth of the channel, and varies mostly with distance along-channel, implying that diffusion (rather than advection) is dominating transport. For u_f=5cm/s, the top, outflowing layer remains fresh, indicating that the transport dominated by advection. This is somewhat counter to my expectations that turbulent mixing should be stronger for the higher velocity cases. This result should be investigated and discussed. As for comment (1), plots of the turbulent diffusivities for each case may be enlightening.
4. The different BC between the melt-enabled and no-melt cases makes it impossible to attribute changes in the intrusion to the effect of melting alone. An additional case is needed with free slip and no melting to better separate these effects.
5. There are numerous different simulations included in the paper, however, few where one variable is systematically changed. This makes the paper/results hard to follow at times. A results table or bar chart showing how the intrusion distance changes across simulations would go some way to addressing this.
6. No salinity effect on melt. In the simulations, the interface is salty (Fig 3E) which will act to depress the freezing temperature and therefore the interface temperature, altering the melt rate. The authors should state why they have not included this extremely important effect (in more detail than “model limitations” line 437). In addition, some discussion of the likely effect of this simplification on the results is needed.
7. Double diffusive framework. Equation 9 is the same as Equation 4 (with marginally different thermal diffusivity) and is therefore not a parameterisation to be tested, rather a re-stating of your melt BC. I propose removing Fig 6B and the “diffusive melt” in Fig 5. Notation m_DC and m_DDC (which are interchangeably used) are inappropriate here.
8. Transfer velocities/Stanton numbers (Figure 6 )– I don’t really understand the purpose of C and D. All D shows is the weak dependence of the Stanton number on ustar at low ustar (see the denominator of (7)), and all C shows is the same but multiplied by ustar. Since the S99 parameterisation can’t accurately predict melting for your simulations, this (trivial) result becomes misleading, as readers may think that it lends support to a certain Stanton number being useful more broadly for modelling ice-ocean interactions. If you want to compare the Stanton number of your simulations to those found in the literature, you need to rearrange equation (6) less the conductive ice flux to solve for the Stanton number (and the (unitless) transfer coefficient \$Gamma_T$) using your model output melt rate, temperature etc.

Other Comments

• Line 83 – move references (carter…) after “non-summer months”
• Line 97 - “Our vertical domain size is at the upper bound of the viscous sublayer length scale that could exist between a well-mixed boundary layer and the ice”– I don’t know what this means, please clarify.
• Line 154 – should be dT/dy in your coordinates
• Line 177 – “vary by a factor of 1000 in response to the range of input velocities tested here”.
• Line 179 – insert comma before “increased”
• Line 180 – “For the middle freshwater velocity (Figure 2 C),” - should this be 2A?
• Line 182 – “To contrast the effects of turbulent mixing, we tested a laminar flow case with no turbulent mixing (green line Figure 2 A) and saw no meaningful difference in intrusion distance.” It would be good to compare and contrast this case more, i.e. how different is the T/S structure? If it’s not different, then presumably turbulent mixing is not occurring in the channel.
• Line 197/8 – refs should not be in parentheses.
• Line 212 – How is Cd calculated? i.e. at what value of y are ustar and u evaluated? Also, as mentioned earlier I think Cd is a main result and this figure should come to the main text.
• Line 241 – “Reduction in velocity gradients arises from an increase in stratification, suppressing turbulence, and the kinematic boundary condition being a velocity inlet and not a no-slip wall.” – as per major comment 4, actually you can’t isolate these effects currently and a no-slip no-melt case is needed.
• Line 256 –Predominantly horizontal motion?
• Line 261 – again, you can’t attribute this to increased stratification (which MAY decrease shear/drag) because you haven’t isolated the effect of the no-slip BC (which WILL decrease shear/drag)
• Line 286 – Diffusive convective melting is also involves convection driven by cooling and is different from diffusive melting. See Martin & Kauffman (1977) section 3 for diffusive melting and Martin & Kauffman (1977) section 4 for diffusive-convective melting.
• Fig 3 – Is any averaging (time or space) done to obtain these profiles?
• Fig 3 - It would be great to see the systematic change in the intrusion as u_f is increased with a series of side-by-side plots, rather than having to move back and forward from figs 3 to A8.
• Line 319 – I would not say that the flow is weak at 5,10 cm/s. However, the height of the channel is very small, so the Reynolds number will be small. Again, Re or other turbulence metrics are needed.
• Line 397 – how is the reduced gravity calculated here? Based on the density difference between the freshwater and saltwater?
• Line 427 – “If we anticipate viscous effects to dominate in seawater intrusions under grounded ice, then using the thermal and haline molecular diffusivities as so-called “transport velocities” would be appropriate, similar to the diffusive-convective framework presented above” – The units (m^2/s vs m/s) are not consistent. To turn the molecular diffusivity into a transport velocity, you need a lengthscale, i.e. the width of the diffusive sublayer. That’s the hard part which is not addressed here!
• Figure 6 caption – The Stanton number should be γT /U, not γT /Cd
• Figure 6 caption last line – are both the dashed and two solid lines from Washam et al 2023?
• Line 490 – “bolstering the idea that grounding zones are subglacial estuaries (Horgan et al., 2013).
• Line 566 – seems like a note-to-self.
• Line 594 – either “given by (3)” or “given by equation 3”
• The appendix is quite bloated. I think some of the figures could be relegated to SI and some should go to the main text. For example, I think Fig A5 is a key result and should go in the main text. Fig A1 could go in SI.
• Table A1 – inconsistent unit formatting (italics/roman). Look to TC style guide, or use roman which is typical. Units for theta should just be degrees.
• Table A2 – again, unit formatting. H [cm]
• Figure A5 – the colours are too hard to tell apart.
• Figures A9 & A10. The figures are labelled “law of the wall” but no interpretation is offered. What is the black vertical line and what does it mean? What would the profiles be expected to look like if a log layer was present? What portion of the flow is being shown? Is y+=0 at the top or bottom of the domain? Do we see a viscous boundary layer, i.e. u+=y+? in addition, it would be much more helpful if y+ was on the y axis, since that’s how the model is set up and how all your other profile plots are oriented.
• Reference 1 (Adusumilli) seems incomplete