the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 22 Aug 2025
Effects of subgrid-scale ice topography on the ice shelf basal melting simulated in NEMO-4.2.0
Abstract. At the interface between the ocean and the ice shelf base, in the framework of the shear-controlled melt parameterisation, the ice melts due to combined actions of temperature, salinity and friction velocity. In the NEMO ocean model, the friction velocity is usually computed based on a constant drag coefficient and an ocean velocity averaged vertically within a distance from the ice, which is often referred to as the Losch layer. Instead, in this study, we use a logarithmic approach, where a constant hydrographic roughness length detetermines the drag coefficient through the law of the wall and the horizontal current speed is sampled in the 1st wet cell. The aim is to reduce the vertical resolution dependency, to homogeneise the sampling of horizontal current speed between the thermodynamic and the dynamic drag equation and to enable the use of a variable drag coefficient based on the subgrid-scale (or unresolved) ice shelf basal topography. The motivation behind a variable drag based on the topography comes from observations showing that regions with rough topographic features such as crevasses or basal melt channels experience more melts than flat ones. We compare different experiments in a configuration of Amundsen Sea at 1/12°. We find that our approach is less sensitive (6 % melt rates difference) to a coarser vertical resolution, such as the one used in global Earth System Models, than the Losch layer approach (22 % melt rates difference). We also find that it succeeds in reproducing higher melt rates in rougher regions while keeping total ice shelf melt rate within the observation range. Finally, to assess the effect of increasing ice shelf damage, we tested the sensitivity of a higher hydrographic roughness length. If the roughness of all the ice shelf grid points were to increase to the highest value currently observed, the overall ice shelf melting would increase by 16 %. This suggests the possibility of a positive feedback in which more melting leads to more ice damage and increased roughness, in turn increasing melt rates.
Competing interests: At least one of the (co-)authors is a member of the editorial board of The Cryosphere. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.-
Notice on discussion status
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
-
Preprint
(26980 KB)
-
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(26980 KB) - Metadata XML
- BibTeX
- EndNote
- Final revised paper
Journal article(s) based on this preprint
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2025-2866', Anonymous Referee #1, 25 Sep 2025
-
AC2: 'Reply on RC1', Dorothée Vallot, 15 Dec 2025
We thank the reviewer for its very good points, which help the manuscript to be clearer and provides further details and simulations to straighten the conclusions.
We noticed, however, that there has been an error in the supplementary figure which uses a logarithm of base 10 to draw the curve of coefficient of drag against the cell thickness while it should be a natural logarithm. This is an important note as part of the comments are based on this graph. In our supplementary material, we draw the same curve(s) with a natural logarithm and for different z. As expected, the result is much different.
Top-level comments
1) The main advance of the paper is in proposing the use of a formulation of the drag coefficient that is dependent upon the grid cell height (equation 3). However, this function is only very weakly dependent upon this height until the grid cell is below 2 metres or even 1 metre (see figure attached as PDF). This seems to be implying, for me, that for most purposes the assumption of a spatially constant drag coefficient should be OK. I feel that this key point should be mentioned. However, and more importantly for the paper, it leads me to wonder where the difference between the REF and LOG simulations comes from. At the moment I think that these differences are probably arising because REF uses the Losch layer of 30m to drive melting, while LOG uses the top wet cell. I propose that they are not caused by the change in the drag coefficient. Clearly, this distinction is central to the paper’s conclusions. I am reluctant to suggest any new simulations, but I think the only way to determine this is to run simulations with the melting parameterisation driven by the top wet cell, not the Losch layer, but using a constant drag coefficient (for both 75 and 121 levels). If these simulations look like REF then the paper’s conclusions stand – the change in the drag coefficient formulation is important. If these simulations look like LOG then the paper’s conclusions should be redirected towards studying the numerical stencil of the melting scheme, rather than the drag formulation. At the moment I think the latter point might be true.
This comment is partly based on the supplementary figure provided by the reviewer, which is not correct as flagged in the comment above. But it is still partly valid and we want to answer it.
This is a very good and important point that is being raised here and this is true that the manuscript lacked a deeper reasoning around it. So I thank the reviewer to give us the opportunity to give more explanation on the use of a variable roughness length.
To start with, I don’t fully agree on the first assertion that the function is only weakly dependent upon the height of the cell and a constant drag coefficient should be ok. From my experience, melt rates are very sensitive to the choice of the coefficient of drag. So using Cd = 0.0025 will provide very different results from Cd = 0.0015 in the REF experiments. By looking closer to the curve (see supplement figure 1a), for z0=3e-4 m (which is the value used in the LOG experiments), Cd = 0.0025 corresponds roughly to z = 1 m and Cd = 0.0015 to z = 15 m. One should keep in mind that Nemo uses the Arakawa C-type grid so that temperature salinity and velocity, for instance, are computed at h/2 and not at h (h being the cell thickness) so that z = h/2. In other words, z=15 m corresponds to a cell of thickness 30 m, exactly the Losch layer thickness.
This brings me to the next, and most important, point stating that the differences between both experiments come from the computing of friction velocity in the Losch layer or in the first wet cell. We actually already had run a simulation with a constant Cd = 0.0025 and the friction velocity computed in the first wet cell for 121 layers (NOLOSCH121) so we just made another simulation with 75 layers (NOLOSCH75) to have a complete picture. The results are actually different than what the reviewer was expecting. For more details, please have a look at the supplementary Table 1. With 121 layers, melt rates for LOSCH121 are in general closer to LOG121 than to REF121, particularly for Thwaites (0.6% against 55.6%) and Crosson (4.2% against 17.2%), other ice shelves having similar differences. Melt rate difference for all ice shelves is -7.7% with REF121 and -10% for LOG121 (because of Getz ice shelves). On the contrary, with 75 levels, the difference is 0.1% with REF121 and 13% for LOG121. So in summary, NOLOSCH is closer to LOG than to REF with 121 layers and much closer to REF than to LOG with 75 layers. Moreover, and most importantly, the difference between NOLOSCH75 and NOLOCH121 (32.2% for all ice shelves) is much higher than between LOG75 and LOG121 (6%) and even REF75 and REF121 (22%). This means that it is absolutely not equivalent to have a constant or variable coefficient of drag in the first wet cell.
To understand better why melt rates from NOLOSCH are closer to LOG when there is 121 levels and to REF with 75 levels, I plotted the histogram of coefficient of drag (figure 1b) and the histogram of cell half thickness (figure 1c) for all ice shelves for 121 and 75 levels. For 121 levels, cell half thickness is always smaller than 10 m and the coefficient of drag is between 0.0014 and 0.0026 while for 75 levels the cell half thickness can be up to 30 m and the coefficient of drag down to 0.001. In other words, LOG121 coefficient of drag is closer to 0.0025 than LOG75 and therefore NOLOSCH121 melt rates are closer to LOG121 although the difference is not that strong for all ice shelves. On the contrary, LOG75 coefficient of drag can be much lower than 0.0025 while the first wet cell thickness is closer to the Losch layer thickness used in REF75. For that reason, NOLOSCH75 is closer to REF75.
In a nutshell, because of the high range of cell thickness (from 1 m to 20 m for 121 levels and up to 60 m for 75 levels), having a constant of variable coefficient of drag on the first wet cell is not equivalent and the conclusion of the paper still holds.
2) There is a conceptual issue around this application of the law of the wall formulation (equation 3). This is designed for constant-stress boundary layers adjacent to an interface that is hydraulically smooth. It can cope with a variable roughness length, as used here, but only as long as that topographic roughness is not larger than the viscous boundary layer depth, such that it would start to disturb the boundary layer with lee eddies, etc. This is apparent in the paper, because the roughness lengths used are tiny - of order 0.1 mm. However, the topographic features present beneath ice shelves and being invoked here are much larger than this – e.g. order 10m in the case of terraces, and order 100m in the case of rifts, ice blocks, and basal channels. This very large topography will completely disturb the boundary layer and so should not be expected to be well represented as a sub-millimetre scale perturbation. Extreme topography on this scale would not be expected to smoothly increase the stress transfer through the boundary layer. This kind of roughness could not even be simply characterised as ‘form drag’, let alone the viscous drag assumed here. The Hughes paper cited by the authors is a nice illustration of this. This is certainly the roughness that the authors are invoking, as it is only these very large scale/large amplitude features that will be present in the surface DEMs being used here. Given this uncertainty, I am not sure what to recommend. It is possible that a constant drag coefficient is even a better choice, in the face of the extreme uncertainty around this problem. I fully support the authors’ goal of incorporating roughness into melting formulations, but the application of this law is not clear. I don’t know how the authors might like to proceed. The best approach here probably depends on whether the drag law makes a difference (point 1 above).
Again what is raised here is a very valid point, which has been a source of ongoing debates in between co-authors. We tried to show this limit in section 3.5 but this could be better flagged following the reviewer’s suggestion. One of the sentence in this section is: “The purpose of this study was mostly to use a proxy to consider the spatially variable roughness of the ice shelf base.” and a better law could be further developed and probably should.
In any case, as shown in the answer above and given the conclusion of the paper, a constant coefficient of drag is not a better choice.
3) I think the logic of the new drag formulation is this: When the top cell is thicker than the log layer it needs to be parameterised, but when the cell is thinner, it will start to resolve the log layer, and hence the parameterisation should be reduced. But, does the model resolve the log layer when the cells become thin? I am not sure we know this. I think it depends upon the vertical viscosity parameterisation in the model, and the implementation of the momentum boundary conditions. In real fluids the log layer profile arises because the eddy viscosity increases linearly with distance from the interface. I don’t believe this eddy viscosity parameterisation is usually used in ocean models. There is also an important corollary here: when the vertical resolution is varied, is the vertical viscosity varied as well? With a higher resolution it would be typical to reduce the viscosity. That would lead to a smaller boundary layer, which is harder for the model to resolve. So it is not clear that a finer grid resolution will better resolve the log layer, as assumed in this new parameterisation.
What you're mentioning is probably true. NEMO, and our configuration doesn't resolve the log-layer (cell thickness is usually more than 1m). It never becomes fine enough to fully resolve the log layer. It could be simulated by varying the eddy viscosity, as you mention, it could be improved by playing with VVL, or it could be resolved totally differently. This studies focuses on one parametrization of the log-layer. We think that the improvement this parametrization is still valuable, and while it could be improved, this is outside the scope of this paper.
4) The reduced sensitivity to vertical grid resolution changes is a key feature of the new approach. However, I wonder if the higher sensitivity of the old approach is due to the way in which the ‘Losch’ layer is implemented in NEMO. In the Losch (2008) paper, this approach was only used for temperature and salinity, not velocities. More importantly, the Losch layer thickness was equal to the top grid cell thickness, not the fixed 30 m that is used here. I wonder if the 22% sensitivity of the old approach here is caused because NEMO uses a fixed 30m Losch layer thickness? Clearly this is going to mean that as the cell thickness increases, fewer cells will be in the Losch layer. But with the original Losch approach, which I think is used in other models, there will always be 2 cells in the Losch layer. Given that the vertical viscosity may change as the grid resolution is changed (?), those 2 cells may contain a fixed fraction of the changing boundary layer thickness. Thus it could be speculated that the original Losch approach would be less grid dependent than the approach apparently used in NEMO. This implies that the new parameterisation may not be as beneficial in models that follow the original Losch approach. Gwyther et al (2020) document substantial variety in the way the models approach this. I think the only way to address this would be with additional simulations following the standard Losch approach, though again I am reluctant to suggest further simulations.
This may be true and would need to be further studied with yet another parameterisation of the Losch layer. We believe this should be a study on its own and may consider it in another paper. Nonetheless, this is important to mention it in the current article: the comparison is only valid in the current version of the Losch layer approximation in NEMO and would need further work to establish the vertical sensitivity of the original Losch layer approach.
5) A key conclusion is that increasing the roughness increases the melting, which could mean there is a feedback between increased ice shelf damage and increased ocean melting. I like this result, but isn’t the top-level result already obvious? If we increase the drag coefficient, melting will increase, as shown by previous studies, and so this conclusion should be stated in relation to this existing knowledge. Where this study could offer a new finding is in the precise sensitivity of melting to changes in roughness – i.e. the curve of melting vs z0. The authors could consider plotting a curve of drag coefficient versus z0 to show how melting might change as the roughness changes? I assume it is highly nonlinear in some way, which would be an interesting finding. However, there is of course the high uncertainty in whether this is the right drag law to use (see above), so perhaps the authors may feel this is not useful. There are also lots of other considerations in terms of how melting might change, such as how the velocity changes with roughness.
Yes, higher roughness leading to higher melt is a tautologism given the law of the wall and the 3-equation melt parameterisation. But as mentioned, having a proxy for the roughness offers a better understanding of its sensitivity, of course with the reserve that it is the right law to use but this can be flagged. In supplementary figure 2, the coefficient of drag is plotted against ln(z0) for different cell half thickness z. This is only giving a partial picture of what the melt rates would look like because of its dependency on velocity, temperature, salinity. Nonetheless, one may see it as a coarse way to speculate on the roughness snesitivity.
Smaller points
General: With partial cells, do the authors apply a minimum cell thickness? I believe that is standard, to avoid numerical instability as the cells get very thin. If, for example, a cell cannot go below 1m thickness, this will reduce the influence of the new drag law even further (see above)
In NEMO, the partial step thickness has a minimum. In our configuration, it is set to be the minimum between 25 m and 20% of the level thickness. For both 121 and 75 layers, the minimum cell half thickness is around 0.5 m so a cell thickness of 1 m. As discussed above, the influence of the law is still important.
Line 95: ‘upper and lower’ is confusing
for u, we can use the east and west side of the cell and for v, we can use the north and south side of the cell.
Line 100: why the subscript 3 ?
This is the way it is called in NEMO. The subscript 3 refers to the vertical component while 1 and 2 refer to the horizontal component. It seems to be worth mentioning in the text.
Figure 2a: I’m afraid I found this really confusing because the cells are all in different places. I think this is because the z* coordinate means that all of the different levels are at different heights and that is randomly applied here, though I am not sure. Could it possibly be redrawn for z coordinates, since that expresses the key difference here (Losch vs 1st wet cell), and then just note that actually z* is used? Also, the cells next to the ice only have 1 open horizontal face (right), with their left boundary being ice. I think it is far more common that they would have two open horizontal faces, with ice only on the top, so maybe better to draw that? This emphasises the need to average over the 4 u and v values on the cell edges.
The purpose of the figure was double: understand what a partial cell of the z*-coordinate system was and show the difference between the Losch layer and the first wet cell. It might be too much info in one picture so we could have instead two schemes, each showing one of the concepts.
Concerning the ice boundary, the purpose was to show the need to average different levels of u and v but we could add a case with flat ice boundary to show the possibility.
Figure 3a: Similarly, I found this really hard to interpret. Maybe select the depth range of ice shelves and plot cells on a linear vertical axis? Maybe that is not possible. Give a table with cell heights at chosen depths (10m, 100m, 500m, 10000m)? Also it looks like the 121 levels were deliberately chosen to represent ice shelf cavities better – is that true? i.e. they are not a uniform increase in resolution, relative to 75 cells.
Yes it is true, the 121 vertical resolution has a fairly constant nominal vertical resolution between 100 and 1000m depth (ie level thickness between 20 and 30m). This should be mentioned.
The depth range of ice shelves is in plain color while other depths are more transparent. We can add a table as suggested so it is clearer.
Line 127: I’m afraid I didn’t understand how this roughness measure was calculated, and this is important. I assume that in the NEMO cell there are many DEM cells. How are all of these cells compared to arrive at R?
Line 127 states:
“To calculate the “topographic roughness“, R, at the base of the ice shelf, we use a simple definition: the largest inter-cell difference between a central pixel and its surrounding cell (Margaret F. J. Wilson Brian O’Connell and Grehan, 2007).”
It means that for each resolution and for each grid cell we compare the difference between the depth of this cell and the ones of the surrounding cells and take the largest difference. Afterwards, to get a value in the Nemo grid, we use a conservative regridding.
Line 134: The mention of 2 m resolution here is confusing.
This is the highest resolution but could be confusing with respect to the following part of the paragraph. So it can be removed as it is not necessary.
Line 170/Table 1: I was confused about this roughness length, because it gives a drag coefficient much higher than 0.0025 for all grid cells (see figure above) and yet it gives the same melting overall. I think this confirms that it is not the changed drag coefficient formulation that is important, it is the switch to using the top cell velocity.
This comment is based on the supplementary figure provided by the reviewer, which we flagged to be incorrect.
Line 221: this explanation does not seem to work for Thwaites?
Yes this is true in the Eastern sector of Thwaites ice shelf close to the grounding line. This should be mentioned.
Paragraph 1 on page 14: I’m sorry, I was not able to follow the 3 points of explanation offered in this paragraph.
Could be rephrased like this:
“This is due to the fact that, in the current version of NEMO, using the Losch layer approach, the mean horizontal current speed, used in the calculation of the friction velocity, is the magnitude of the components averages and not the average of the magnitudes”
Figures 9 and 10 and sections 3.1 and 3.2: I could not find a physical explanation of why the melt rates were different in these cases. There are interesting patterns of difference, but these were not explained in terms of the underlying velocities in the different regions, and how well they were represented by the different drag laws, and the different vertical resolutions.
Melt rates depend of velocity, salinity and temperature following non-linear equations and the differences of melt rates between LOG and REF (figure 9) is also dependent on the coefficient of drag and the resolution in the Losch layer and in the first wet cell. So differences can come from different factors at the same time. What we can see is that there is the difference pattern of melt rates is similar to the one of friction velocity, being the main contributor. By looking at the treatment of the mean horizontal current speed, UM, used in the calculation of the friction velocity, we see that the difference pattern is similar in general (not everywhere) so we attribute most of the difference to the tretment of UM.
The difference between vertical resolution should now be more complete with the NOLOSCH experiments mentioned in the first comment.
Equation (B4): bottom equation has an extra term C, which I don’t follow. Is this a constant background velocity?
It was a mistake, this will be removed.
-
AC2: 'Reply on RC1', Dorothée Vallot, 15 Dec 2025
-
RC2: 'Comment on egusphere-2025-2866', Robert Larter, 03 Oct 2025
I consider this to be an interesting study addressing a fundamental problem in cryospheric science, the parameterisation of ice shelf basal melting. Having said that, I am not a numerical modeller or mathematician and think it is essential that the manuscript is also reviewed by someone with expertise in these fields (I have not yet looked at the other review I was notified about because I want to provide completely independent comments).
Recent studies have revealed that ice shelf basal melting is more spatially and temporally variable than previously appreciated, so new ways of modelling the melting, particularly for areas where ice shelf basal topography is complex, need to be explored. Therefore, I am pleased to see this account of such a study.
In general, the manuscript appears to me to clearly articulate results of a thorough exploration of different approaches to parameterising basal melt. I have no major comments, but I have annotated a copy of the manuscript with many minor ones. I have also included a few of the more important of these comments below.
A general comment is that annotations on several of the figures are in a font that is too small if they are going to be published at the same size as in the review copy. Colour scale bars in several figures are annotated with values expressed with exponents, whereas I think converting these to more simple numbers would make them easier to understand for a greater number of readers.
Labels need to be added to identify the different parts of figures 2, 3 and 4 that are referred to in the captions.
After studying Fig. 2a and the related text it is still not completely clear to me where on the boundaries of cells various parameters are being computed. I infer that the blue box on the left of the figure is depicting a plan view of a cell since it is labelled with x and y axes, but does the position of the T-point that is depicted mean that is on the upper or lower face of the cell, or both, or in the middle of the cell? The V-point is described in the text as being “at the center of the upper and lower grid face”, but if the blue box labelled with x and y axes is indeed a plan view, then the V-point appears to be on one of the lateral faces. I think it would be helpful if the authors could try to further clarify their explanation of where these different points lie by modifying the text or the figure, or both.
Fig. 2b is the first place in the paper that maps of the different ice shelves appear. It would be useful to include a regional map in the figure, or in an earlier figure, showing where the areas represented by these local maps lie. A similar map to the one in Figure 4b could be used for this purpose. It would also be helpful to add distance scale bars to each of the local maps to make clear that they are at different scales.
The detail I am most puzzled about is the apparent near random distribution of first wet cell heights shown in Fig. 2b. Fig. 2a and Fig. 3a indicate that first wet cell height increases with depth of the ice shelf base. Therefore, I would expect to see a trend towards more purple and black cells, indicating greater cell height, nearer the grounding lines. Unless an error has been made in the production of the figure, I think a further explanation of the apparent near random distribution of first wet cell heights is needed.
I am also puzzled about the difference between the values of the constants given in the “roughness length” column in Table 1 and the values of the constant α given in the caption to Figure 5. Perhaps I have misunderstood something, but it seems to me that these constant values for a particular experiment should be the same in both places.
There are a couple of places where variable terminology is used to describe the same thing, e.g. “forcing” and “driving”, and “hydrographic roughness” and “hydrographic roughness length”. More consistent terminology would make the paper more accessible to non-specialists.
R D Larter
3rd October 2025
-
AC1: 'Reply on RC2', Dorothée Vallot, 15 Dec 2025
We would like to thank the reviewer for its pertinent comments and we reply to them in the text (review comments in bold and reply in plain).
I consider this to be an interesting study addressing a fundamental problem in cryospheric science, the parameterisation of ice shelf basal melting. Having said that, I am not a numerical modeller or mathematician and think it is essential that the manuscript is also reviewed by someone with expertise in these fields (I have not yet looked at the other review I was notified about because I want to provide completely independent comments).
Recent studies have revealed that ice shelf basal melting is more spatially and temporally variable than previously appreciated, so new ways of modelling the melting, particularly for areas where ice shelf basal topography is complex, need to be explored. Therefore, I am pleased to see this account of such a study.
In general, the manuscript appears to me to clearly articulate results of a thorough exploration of different approaches to parameterising basal melt. I have no major comments, but I have annotated a copy of the manuscript with many minor ones. I have also included a few of the more important of these comments below.
We agree with all comments in the annotated file and are correcting the manuscript accordingly.
A general comment is that annotations on several of the figures are in a font that is too small if they are going to be published at the same size as in the review copy. Colour scale bars in several figures are annotated with values expressed with exponents, whereas I think converting these to more simple numbers would make them easier to understand for a greater number of readers.
Good remark that is taken into account.
Labels need to be added to identify the different parts of figures 2, 3 and 4 that are referred to in the captions.
True. These will be added in the next version.
After studying Fig. 2a and the related text it is still not completely clear to me where on the boundaries of cells various parameters are being computed. I infer that the blue box on the left of the figure is depicting a plan view of a cell since it is labelled with x and y axes, but does the position of the T-point that is depicted mean that is on the upper or lower face of the cell, or both, or in the middle of the cell? The V-point is described in the text as being “at the center of the upper and lower grid face”, but if the blue box labelled with x and y axes is indeed a plan view, then the V-point appears to be on one of the lateral faces. I think it would be helpful if the authors could try to further clarify their explanation of where these different points lie by modifying the text or the figure, or both.
The other reviewer also had a comment on this figure and a new version we hope more understandable is in the supplementary.
We use now a 3D version of the cell. The T point is actually in the middle of the cell. For u, we can use the east and west side of the cell and for v, we can use the north and south side of the cell to be clearer.
Fig. 2b is the first place in the paper that maps of the different ice shelves appear. It would be useful to include a regional map in the figure, or in an earlier figure, showing where the areas represented by these local maps lie. A similar map to the one in Figure 4b could be used for this purpose. It would also be helpful to add distance scale bars to each of the local maps to make clear that they are at different scales.
The detail I am most puzzled about is the apparent near random distribution of first wet cell heights shown in Fig. 2b. Fig. 2a and Fig. 3a indicate that first wet cell height increases with depth of the ice shelf base. Therefore, I would expect to see a trend towards more purple and black cells, indicating greater cell height, nearer the grounding lines. Unless an error has been made in the production of the figure, I think a further explanation of the apparent near random distribution of first wet cell heights is needed.
This is because we use the z*-coordinate system with partial steps. In the paragraph 2.1.1, it is explained in the first sentence: “We use a curvilinear z*-coordinate system with partial steps to adjust the thickness of grid cells adjacent to the sea floor (index kb ) or ice shelf draft (index kt ).” The thickness of the cell depends on the level position but also where the boundary is. This is what we try to show in figure 2a.
We can add an explanation in the figure caption: “First wet cell half thickness below ice shelves in the ASE configuration with 121 layers: Pine Island, Thwaites, Dotson-Crosson and Getz from left to right. The cell thickness is adjusted to the ice shelf draft through partial step.
I am also puzzled about the difference between the values of the constants given in the “roughness length” column in Table 1 and the values of the constant α given in the caption to Figure 5. Perhaps I have misunderstood something, but it seems to me that these constant values for a particular experiment should be the same in both places.
This is a mistake from our side. It is from an older version where the alpha value was wrong.
There are a couple of places where variable terminology is used to describe the same thing, e.g. “forcing” and “driving”, and “hydrographic roughness” and “hydrographic roughness length”. More consistent terminology would make the paper more accessible to non-specialists.
Good point. We use hydrographic roughness length everywhere.
R D Larter
3rd October 2025
-
AC1: 'Reply on RC2', Dorothée Vallot, 15 Dec 2025
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2025-2866', Anonymous Referee #1, 25 Sep 2025
This paper presents a modification to the way that melting of glacial ice shelves is implemented in the NEMO ocean model, with the intention of capturing the effects of variable ice base roughness on melting. This is certainly an important aim, since current melting parameterisations assume that the ice base is uniformly rough, which is known to be inaccurate. The current models also assume that the boundary layer is evenly resolved by the model, which is also inaccurate. The recommended approach aims to address these points by implementing a logarithmic boundary layer speed profile into the model.
The paper concludes that i) the new implementation is less sensitive to changes in vertical grid resolution in the model, and ii) variable spatial roughness is an important control over ice shelf melting patterns, and iii) any temporal change in roughness will cause a related change in melting.
My personal opinion is that the new melting implementation may be better than the current NEMO approach, which is based on an ad-hoc boundary layer thickness. Also, it is clearly desirable that the stress formulations used in the melting parameterisation and momentum solver are consistent with each other. However, this is quite subjective and arguable in several areas, and I do have some substantial points to make about the current manuscript, requiring major revisions. Ideally, these would require new simulations to address, though I am reluctant to suggest this. I think this work deserves to be published, ultimately. I think this new approach raises some discussion points that will be very valuable to the community. I would like the authors to consider the points below and I am happy for them to disagree with my view points, as long as the discussion is reflected in the paper.
Top-level comments
1) The main advance of the paper is in proposing the use of a formulation of the drag coefficient that is dependent upon the grid cell height (equation 3). However, this function is only very weakly dependent upon this height until the grid cell is below 2 metres or even 1 metre (see figure attached as PDF). This seems to be implying, for me, that for most purposes the assumption of a spatially constant drag coefficient should be OK. I feel that this key point should be mentioned. However, and more importantly for the paper, it leads me to wonder where the difference between the REF and LOG simulations comes from. At the moment I think that these differences are probably arising because REF uses the Losch layer of 30m to drive melting, while LOG uses the top wet cell. I propose that they are not caused by the change in the drag coefficient. Clearly, this distinction is central to the paper’s conclusions. I am reluctant to suggest any new simulations, but I think the only way to determine this is to run simulations with the melting parameterisation driven by the top wet cell, not the Losch layer, but using a constant drag coefficient (for both 75 and 121 levels). If these simulations look like REF then the paper’s conclusions stand – the change in the drag coefficient formulation is important. If these simulations look like LOG then the paper’s conclusions should be redirected towards studying the numerical stencil of the melting scheme, rather than the drag formulation. At the moment I think the latter point might be true.
2) There is a conceptual issue around this application of the law of the wall formulation (equation 3). This is designed for constant-stress boundary layers adjacent to an interface that is hydraulically smooth. It can cope with a variable roughness length, as used here, but only as long as that topographic roughness is not larger than the viscous boundary layer depth, such that it would start to disturb the boundary layer with lee eddies, etc. This is apparent in the paper, because the roughness lengths used are tiny - of order 0.1 mm. However, the topographic features present beneath ice shelves and being invoked here are much larger than this – e.g. order 10m in the case of terraces, and order 100m in the case of rifts, ice blocks, and basal channels. This very large topography will completely disturb the boundary layer and so should not be expected to be well represented as a sub-millimetre scale perturbation. Extreme topography on this scale would not be expected to smoothly increase the stress transfer through the boundary layer. This kind of roughness could not even be simply characterised as ‘form drag’, let alone the viscous drag assumed here. The Hughes paper cited by the authors is a nice illustration of this. This is certainly the roughness that the authors are invoking, as it is only these very large scale/large amplitude features that will be present in the surface DEMs being used here. Given this uncertainty, I am not sure what to recommend. It is possible that a constant drag coefficient is even a better choice, in the face of the extreme uncertainty around this problem. I fully support the authors’ goal of incorporating roughness into melting formulations, but the application of this law is not clear. I don’t know how the authors might like to proceed. The best approach here probably depends on whether the drag law makes a difference (point 1 above).
3) I think the logic of the new drag formulation is this: When the top cell is thicker than the log layer it needs to be parameterised, but when the cell is thinner, it will start to resolve the log layer, and hence the parameterisation should be reduced. But, does the model resolve the log layer when the cells become thin? I am not sure we know this. I think it depends upon the vertical viscosity parameterisation in the model, and the implementation of the momentum boundary conditions. In real fluids the log layer profile arises because the eddy viscosity increases linearly with distance from the interface. I don’t believe this eddy viscosity parameterisation is usually used in ocean models. There is also an important corollary here: when the vertical resolution is varied, is the vertical viscosity varied as well? With a higher resolution it would be typical to reduce the viscosity. That would lead to a smaller boundary layer, which is harder for the model to resolve. So it is not clear that a finer grid resolution will better resolve the log layer, as assumed in this new parameterisation.
4) The reduced sensitivity to vertical grid resolution changes is a key feature of the new approach. However, I wonder if the higher sensitivity of the old approach is due to the way in which the ‘Losch’ layer is implemented in NEMO. In the Losch (2008) paper, this approach was only used for temperature and salinity, not velocities. More importantly, the Losch layer thickness was equal to the top grid cell thickness, not the fixed 30 m that is used here. I wonder if the 22% sensitivity of the old approach here is caused because NEMO uses a fixed 30m Losch layer thickness? Clearly this is going to mean that as the cell thickness increases, fewer cells will be in the Losch layer. But with the original Losch approach, which I think is used in other models, there will always be 2 cells in the Losch layer. Given that the vertical viscosity may change as the grid resolution is changed (?), those 2 cells may contain a fixed fraction of the changing boundary layer thickness. Thus it could be speculated that the original Losch approach would be less grid dependent than the approach apparently used in NEMO. This implies that the new parameterisation may not be as beneficial in models that follow the original Losch approach. Gwyther et al (2020) document substantial variety in the way the models approach this. I think the only way to address this would be with additional simulations following the standard Losch approach, though again I am reluctant to suggest further simulations.
5) A key conclusion is that increasing the roughness increases the melting, which could mean there is a feedback between increased ice shelf damage and increased ocean melting. I like this result, but isn’t the top-level result already obvious? If we increase the drag coefficient, melting will increase, as shown by previous studies, and so this conclusion should be stated in relation to this existing knowledge. Where this study could offer a new finding is in the precise sensitivity of melting to changes in roughness – i.e. the curve of melting vs z0. The authors could consider plotting a curve of drag coefficient versus z0 to show how melting might change as the roughness changes? I assume it is highly nonlinear in some way, which would be an interesting finding. However, there is of course the high uncertainty in whether this is the right drag law to use (see above), so perhaps the authors may feel this is not useful. There are also lots of other considerations in terms of how melting might change, such as how the velocity changes with roughness.
Smaller points
General: With partial cells, do the authors apply a minimum cell thickness? I believe that is standard, to avoid numerical instability as the cells get very thin. If, for example, a cell cannot go below 1m thickness, this will reduce the influence of the new drag law even further (see above)
Line 95: ‘upper and lower’ is confusing
Line 100: why the subscript 3 ?
Figure 2a: I’m afraid I found this really confusing because the cells are all in different places. I think this is because the z* coordinate means that all of the different levels are at different heights and that is randomly applied here, though I am not sure. Could it possibly be redrawn for z coordinates, since that expresses the key difference here (Losch vs 1st wet cell), and then just note that actually z* is used? Also, the cells next to the ice only have 1 open horizontal face (right), with their left boundary being ice. I think it is far more common that they would have two open horizontal faces, with ice only on the top, so maybe better to draw that? This emphasises the need to average over the 4 u and v values on the cell edges.
Figure 3a: Similarly, I found this really hard to interpret. Maybe select the depth range of ice shelves and plot cells on a linear vertical axis? Maybe that is not possible. Give a table with cell heights at chosen depths (10m, 100m, 500m, 10000m)? Also it looks like the 121 levels were deliberately chosen to represent ice shelf cavities better – is that true? i.e. they are not a uniform increase in resolution, relative to 75 cells.
Line 127: I’m afraid I didn’t understand how this roughness measure was calculated, and this is important. I assume that in the NEMO cell there are many DEM cells. How are all of these cells compared to arrive at R?
Line 134: The mention of 2 m resolution here is confusing.
Line 170/Table 1: I was confused about this roughness length, because it gives a drag coefficient much higher than 0.0025 for all grid cells (see figure above) and yet it gives the same melting overall. I think this confirms that it is not the changed drag coefficient formulation that is important, it is the switch to using the top cell velocity.
Line 221: this explanation does not seem to work for Thwaites?
Paragraph 1 on page 14: I’m sorry, I was not able to follow the 3 points of explanation offered in this paragraph.
Figures 9 and 10 and sections 3.1 and 3.2: I could not find a physical explanation of why the melt rates were different in these cases. There are interesting patterns of difference, but these were not explained in terms of the underlying velocities in the different regions, and how well they were represented by the different drag laws, and the different vertical resolutions.
Equation (B4): bottom equation has an extra term C, which I don’t follow. Is this a constant background velocity?
-
AC2: 'Reply on RC1', Dorothée Vallot, 15 Dec 2025
We thank the reviewer for its very good points, which help the manuscript to be clearer and provides further details and simulations to straighten the conclusions.
We noticed, however, that there has been an error in the supplementary figure which uses a logarithm of base 10 to draw the curve of coefficient of drag against the cell thickness while it should be a natural logarithm. This is an important note as part of the comments are based on this graph. In our supplementary material, we draw the same curve(s) with a natural logarithm and for different z. As expected, the result is much different.
Top-level comments
1) The main advance of the paper is in proposing the use of a formulation of the drag coefficient that is dependent upon the grid cell height (equation 3). However, this function is only very weakly dependent upon this height until the grid cell is below 2 metres or even 1 metre (see figure attached as PDF). This seems to be implying, for me, that for most purposes the assumption of a spatially constant drag coefficient should be OK. I feel that this key point should be mentioned. However, and more importantly for the paper, it leads me to wonder where the difference between the REF and LOG simulations comes from. At the moment I think that these differences are probably arising because REF uses the Losch layer of 30m to drive melting, while LOG uses the top wet cell. I propose that they are not caused by the change in the drag coefficient. Clearly, this distinction is central to the paper’s conclusions. I am reluctant to suggest any new simulations, but I think the only way to determine this is to run simulations with the melting parameterisation driven by the top wet cell, not the Losch layer, but using a constant drag coefficient (for both 75 and 121 levels). If these simulations look like REF then the paper’s conclusions stand – the change in the drag coefficient formulation is important. If these simulations look like LOG then the paper’s conclusions should be redirected towards studying the numerical stencil of the melting scheme, rather than the drag formulation. At the moment I think the latter point might be true.
This comment is partly based on the supplementary figure provided by the reviewer, which is not correct as flagged in the comment above. But it is still partly valid and we want to answer it.
This is a very good and important point that is being raised here and this is true that the manuscript lacked a deeper reasoning around it. So I thank the reviewer to give us the opportunity to give more explanation on the use of a variable roughness length.
To start with, I don’t fully agree on the first assertion that the function is only weakly dependent upon the height of the cell and a constant drag coefficient should be ok. From my experience, melt rates are very sensitive to the choice of the coefficient of drag. So using Cd = 0.0025 will provide very different results from Cd = 0.0015 in the REF experiments. By looking closer to the curve (see supplement figure 1a), for z0=3e-4 m (which is the value used in the LOG experiments), Cd = 0.0025 corresponds roughly to z = 1 m and Cd = 0.0015 to z = 15 m. One should keep in mind that Nemo uses the Arakawa C-type grid so that temperature salinity and velocity, for instance, are computed at h/2 and not at h (h being the cell thickness) so that z = h/2. In other words, z=15 m corresponds to a cell of thickness 30 m, exactly the Losch layer thickness.
This brings me to the next, and most important, point stating that the differences between both experiments come from the computing of friction velocity in the Losch layer or in the first wet cell. We actually already had run a simulation with a constant Cd = 0.0025 and the friction velocity computed in the first wet cell for 121 layers (NOLOSCH121) so we just made another simulation with 75 layers (NOLOSCH75) to have a complete picture. The results are actually different than what the reviewer was expecting. For more details, please have a look at the supplementary Table 1. With 121 layers, melt rates for LOSCH121 are in general closer to LOG121 than to REF121, particularly for Thwaites (0.6% against 55.6%) and Crosson (4.2% against 17.2%), other ice shelves having similar differences. Melt rate difference for all ice shelves is -7.7% with REF121 and -10% for LOG121 (because of Getz ice shelves). On the contrary, with 75 levels, the difference is 0.1% with REF121 and 13% for LOG121. So in summary, NOLOSCH is closer to LOG than to REF with 121 layers and much closer to REF than to LOG with 75 layers. Moreover, and most importantly, the difference between NOLOSCH75 and NOLOCH121 (32.2% for all ice shelves) is much higher than between LOG75 and LOG121 (6%) and even REF75 and REF121 (22%). This means that it is absolutely not equivalent to have a constant or variable coefficient of drag in the first wet cell.
To understand better why melt rates from NOLOSCH are closer to LOG when there is 121 levels and to REF with 75 levels, I plotted the histogram of coefficient of drag (figure 1b) and the histogram of cell half thickness (figure 1c) for all ice shelves for 121 and 75 levels. For 121 levels, cell half thickness is always smaller than 10 m and the coefficient of drag is between 0.0014 and 0.0026 while for 75 levels the cell half thickness can be up to 30 m and the coefficient of drag down to 0.001. In other words, LOG121 coefficient of drag is closer to 0.0025 than LOG75 and therefore NOLOSCH121 melt rates are closer to LOG121 although the difference is not that strong for all ice shelves. On the contrary, LOG75 coefficient of drag can be much lower than 0.0025 while the first wet cell thickness is closer to the Losch layer thickness used in REF75. For that reason, NOLOSCH75 is closer to REF75.
In a nutshell, because of the high range of cell thickness (from 1 m to 20 m for 121 levels and up to 60 m for 75 levels), having a constant of variable coefficient of drag on the first wet cell is not equivalent and the conclusion of the paper still holds.
2) There is a conceptual issue around this application of the law of the wall formulation (equation 3). This is designed for constant-stress boundary layers adjacent to an interface that is hydraulically smooth. It can cope with a variable roughness length, as used here, but only as long as that topographic roughness is not larger than the viscous boundary layer depth, such that it would start to disturb the boundary layer with lee eddies, etc. This is apparent in the paper, because the roughness lengths used are tiny - of order 0.1 mm. However, the topographic features present beneath ice shelves and being invoked here are much larger than this – e.g. order 10m in the case of terraces, and order 100m in the case of rifts, ice blocks, and basal channels. This very large topography will completely disturb the boundary layer and so should not be expected to be well represented as a sub-millimetre scale perturbation. Extreme topography on this scale would not be expected to smoothly increase the stress transfer through the boundary layer. This kind of roughness could not even be simply characterised as ‘form drag’, let alone the viscous drag assumed here. The Hughes paper cited by the authors is a nice illustration of this. This is certainly the roughness that the authors are invoking, as it is only these very large scale/large amplitude features that will be present in the surface DEMs being used here. Given this uncertainty, I am not sure what to recommend. It is possible that a constant drag coefficient is even a better choice, in the face of the extreme uncertainty around this problem. I fully support the authors’ goal of incorporating roughness into melting formulations, but the application of this law is not clear. I don’t know how the authors might like to proceed. The best approach here probably depends on whether the drag law makes a difference (point 1 above).
Again what is raised here is a very valid point, which has been a source of ongoing debates in between co-authors. We tried to show this limit in section 3.5 but this could be better flagged following the reviewer’s suggestion. One of the sentence in this section is: “The purpose of this study was mostly to use a proxy to consider the spatially variable roughness of the ice shelf base.” and a better law could be further developed and probably should.
In any case, as shown in the answer above and given the conclusion of the paper, a constant coefficient of drag is not a better choice.
3) I think the logic of the new drag formulation is this: When the top cell is thicker than the log layer it needs to be parameterised, but when the cell is thinner, it will start to resolve the log layer, and hence the parameterisation should be reduced. But, does the model resolve the log layer when the cells become thin? I am not sure we know this. I think it depends upon the vertical viscosity parameterisation in the model, and the implementation of the momentum boundary conditions. In real fluids the log layer profile arises because the eddy viscosity increases linearly with distance from the interface. I don’t believe this eddy viscosity parameterisation is usually used in ocean models. There is also an important corollary here: when the vertical resolution is varied, is the vertical viscosity varied as well? With a higher resolution it would be typical to reduce the viscosity. That would lead to a smaller boundary layer, which is harder for the model to resolve. So it is not clear that a finer grid resolution will better resolve the log layer, as assumed in this new parameterisation.
What you're mentioning is probably true. NEMO, and our configuration doesn't resolve the log-layer (cell thickness is usually more than 1m). It never becomes fine enough to fully resolve the log layer. It could be simulated by varying the eddy viscosity, as you mention, it could be improved by playing with VVL, or it could be resolved totally differently. This studies focuses on one parametrization of the log-layer. We think that the improvement this parametrization is still valuable, and while it could be improved, this is outside the scope of this paper.
4) The reduced sensitivity to vertical grid resolution changes is a key feature of the new approach. However, I wonder if the higher sensitivity of the old approach is due to the way in which the ‘Losch’ layer is implemented in NEMO. In the Losch (2008) paper, this approach was only used for temperature and salinity, not velocities. More importantly, the Losch layer thickness was equal to the top grid cell thickness, not the fixed 30 m that is used here. I wonder if the 22% sensitivity of the old approach here is caused because NEMO uses a fixed 30m Losch layer thickness? Clearly this is going to mean that as the cell thickness increases, fewer cells will be in the Losch layer. But with the original Losch approach, which I think is used in other models, there will always be 2 cells in the Losch layer. Given that the vertical viscosity may change as the grid resolution is changed (?), those 2 cells may contain a fixed fraction of the changing boundary layer thickness. Thus it could be speculated that the original Losch approach would be less grid dependent than the approach apparently used in NEMO. This implies that the new parameterisation may not be as beneficial in models that follow the original Losch approach. Gwyther et al (2020) document substantial variety in the way the models approach this. I think the only way to address this would be with additional simulations following the standard Losch approach, though again I am reluctant to suggest further simulations.
This may be true and would need to be further studied with yet another parameterisation of the Losch layer. We believe this should be a study on its own and may consider it in another paper. Nonetheless, this is important to mention it in the current article: the comparison is only valid in the current version of the Losch layer approximation in NEMO and would need further work to establish the vertical sensitivity of the original Losch layer approach.
5) A key conclusion is that increasing the roughness increases the melting, which could mean there is a feedback between increased ice shelf damage and increased ocean melting. I like this result, but isn’t the top-level result already obvious? If we increase the drag coefficient, melting will increase, as shown by previous studies, and so this conclusion should be stated in relation to this existing knowledge. Where this study could offer a new finding is in the precise sensitivity of melting to changes in roughness – i.e. the curve of melting vs z0. The authors could consider plotting a curve of drag coefficient versus z0 to show how melting might change as the roughness changes? I assume it is highly nonlinear in some way, which would be an interesting finding. However, there is of course the high uncertainty in whether this is the right drag law to use (see above), so perhaps the authors may feel this is not useful. There are also lots of other considerations in terms of how melting might change, such as how the velocity changes with roughness.
Yes, higher roughness leading to higher melt is a tautologism given the law of the wall and the 3-equation melt parameterisation. But as mentioned, having a proxy for the roughness offers a better understanding of its sensitivity, of course with the reserve that it is the right law to use but this can be flagged. In supplementary figure 2, the coefficient of drag is plotted against ln(z0) for different cell half thickness z. This is only giving a partial picture of what the melt rates would look like because of its dependency on velocity, temperature, salinity. Nonetheless, one may see it as a coarse way to speculate on the roughness snesitivity.
Smaller points
General: With partial cells, do the authors apply a minimum cell thickness? I believe that is standard, to avoid numerical instability as the cells get very thin. If, for example, a cell cannot go below 1m thickness, this will reduce the influence of the new drag law even further (see above)
In NEMO, the partial step thickness has a minimum. In our configuration, it is set to be the minimum between 25 m and 20% of the level thickness. For both 121 and 75 layers, the minimum cell half thickness is around 0.5 m so a cell thickness of 1 m. As discussed above, the influence of the law is still important.
Line 95: ‘upper and lower’ is confusing
for u, we can use the east and west side of the cell and for v, we can use the north and south side of the cell.
Line 100: why the subscript 3 ?
This is the way it is called in NEMO. The subscript 3 refers to the vertical component while 1 and 2 refer to the horizontal component. It seems to be worth mentioning in the text.
Figure 2a: I’m afraid I found this really confusing because the cells are all in different places. I think this is because the z* coordinate means that all of the different levels are at different heights and that is randomly applied here, though I am not sure. Could it possibly be redrawn for z coordinates, since that expresses the key difference here (Losch vs 1st wet cell), and then just note that actually z* is used? Also, the cells next to the ice only have 1 open horizontal face (right), with their left boundary being ice. I think it is far more common that they would have two open horizontal faces, with ice only on the top, so maybe better to draw that? This emphasises the need to average over the 4 u and v values on the cell edges.
The purpose of the figure was double: understand what a partial cell of the z*-coordinate system was and show the difference between the Losch layer and the first wet cell. It might be too much info in one picture so we could have instead two schemes, each showing one of the concepts.
Concerning the ice boundary, the purpose was to show the need to average different levels of u and v but we could add a case with flat ice boundary to show the possibility.
Figure 3a: Similarly, I found this really hard to interpret. Maybe select the depth range of ice shelves and plot cells on a linear vertical axis? Maybe that is not possible. Give a table with cell heights at chosen depths (10m, 100m, 500m, 10000m)? Also it looks like the 121 levels were deliberately chosen to represent ice shelf cavities better – is that true? i.e. they are not a uniform increase in resolution, relative to 75 cells.
Yes it is true, the 121 vertical resolution has a fairly constant nominal vertical resolution between 100 and 1000m depth (ie level thickness between 20 and 30m). This should be mentioned.
The depth range of ice shelves is in plain color while other depths are more transparent. We can add a table as suggested so it is clearer.
Line 127: I’m afraid I didn’t understand how this roughness measure was calculated, and this is important. I assume that in the NEMO cell there are many DEM cells. How are all of these cells compared to arrive at R?
Line 127 states:
“To calculate the “topographic roughness“, R, at the base of the ice shelf, we use a simple definition: the largest inter-cell difference between a central pixel and its surrounding cell (Margaret F. J. Wilson Brian O’Connell and Grehan, 2007).”
It means that for each resolution and for each grid cell we compare the difference between the depth of this cell and the ones of the surrounding cells and take the largest difference. Afterwards, to get a value in the Nemo grid, we use a conservative regridding.
Line 134: The mention of 2 m resolution here is confusing.
This is the highest resolution but could be confusing with respect to the following part of the paragraph. So it can be removed as it is not necessary.
Line 170/Table 1: I was confused about this roughness length, because it gives a drag coefficient much higher than 0.0025 for all grid cells (see figure above) and yet it gives the same melting overall. I think this confirms that it is not the changed drag coefficient formulation that is important, it is the switch to using the top cell velocity.
This comment is based on the supplementary figure provided by the reviewer, which we flagged to be incorrect.
Line 221: this explanation does not seem to work for Thwaites?
Yes this is true in the Eastern sector of Thwaites ice shelf close to the grounding line. This should be mentioned.
Paragraph 1 on page 14: I’m sorry, I was not able to follow the 3 points of explanation offered in this paragraph.
Could be rephrased like this:
“This is due to the fact that, in the current version of NEMO, using the Losch layer approach, the mean horizontal current speed, used in the calculation of the friction velocity, is the magnitude of the components averages and not the average of the magnitudes”
Figures 9 and 10 and sections 3.1 and 3.2: I could not find a physical explanation of why the melt rates were different in these cases. There are interesting patterns of difference, but these were not explained in terms of the underlying velocities in the different regions, and how well they were represented by the different drag laws, and the different vertical resolutions.
Melt rates depend of velocity, salinity and temperature following non-linear equations and the differences of melt rates between LOG and REF (figure 9) is also dependent on the coefficient of drag and the resolution in the Losch layer and in the first wet cell. So differences can come from different factors at the same time. What we can see is that there is the difference pattern of melt rates is similar to the one of friction velocity, being the main contributor. By looking at the treatment of the mean horizontal current speed, UM, used in the calculation of the friction velocity, we see that the difference pattern is similar in general (not everywhere) so we attribute most of the difference to the tretment of UM.
The difference between vertical resolution should now be more complete with the NOLOSCH experiments mentioned in the first comment.
Equation (B4): bottom equation has an extra term C, which I don’t follow. Is this a constant background velocity?
It was a mistake, this will be removed.
-
AC2: 'Reply on RC1', Dorothée Vallot, 15 Dec 2025
-
RC2: 'Comment on egusphere-2025-2866', Robert Larter, 03 Oct 2025
I consider this to be an interesting study addressing a fundamental problem in cryospheric science, the parameterisation of ice shelf basal melting. Having said that, I am not a numerical modeller or mathematician and think it is essential that the manuscript is also reviewed by someone with expertise in these fields (I have not yet looked at the other review I was notified about because I want to provide completely independent comments).
Recent studies have revealed that ice shelf basal melting is more spatially and temporally variable than previously appreciated, so new ways of modelling the melting, particularly for areas where ice shelf basal topography is complex, need to be explored. Therefore, I am pleased to see this account of such a study.
In general, the manuscript appears to me to clearly articulate results of a thorough exploration of different approaches to parameterising basal melt. I have no major comments, but I have annotated a copy of the manuscript with many minor ones. I have also included a few of the more important of these comments below.
A general comment is that annotations on several of the figures are in a font that is too small if they are going to be published at the same size as in the review copy. Colour scale bars in several figures are annotated with values expressed with exponents, whereas I think converting these to more simple numbers would make them easier to understand for a greater number of readers.
Labels need to be added to identify the different parts of figures 2, 3 and 4 that are referred to in the captions.
After studying Fig. 2a and the related text it is still not completely clear to me where on the boundaries of cells various parameters are being computed. I infer that the blue box on the left of the figure is depicting a plan view of a cell since it is labelled with x and y axes, but does the position of the T-point that is depicted mean that is on the upper or lower face of the cell, or both, or in the middle of the cell? The V-point is described in the text as being “at the center of the upper and lower grid face”, but if the blue box labelled with x and y axes is indeed a plan view, then the V-point appears to be on one of the lateral faces. I think it would be helpful if the authors could try to further clarify their explanation of where these different points lie by modifying the text or the figure, or both.
Fig. 2b is the first place in the paper that maps of the different ice shelves appear. It would be useful to include a regional map in the figure, or in an earlier figure, showing where the areas represented by these local maps lie. A similar map to the one in Figure 4b could be used for this purpose. It would also be helpful to add distance scale bars to each of the local maps to make clear that they are at different scales.
The detail I am most puzzled about is the apparent near random distribution of first wet cell heights shown in Fig. 2b. Fig. 2a and Fig. 3a indicate that first wet cell height increases with depth of the ice shelf base. Therefore, I would expect to see a trend towards more purple and black cells, indicating greater cell height, nearer the grounding lines. Unless an error has been made in the production of the figure, I think a further explanation of the apparent near random distribution of first wet cell heights is needed.
I am also puzzled about the difference between the values of the constants given in the “roughness length” column in Table 1 and the values of the constant α given in the caption to Figure 5. Perhaps I have misunderstood something, but it seems to me that these constant values for a particular experiment should be the same in both places.
There are a couple of places where variable terminology is used to describe the same thing, e.g. “forcing” and “driving”, and “hydrographic roughness” and “hydrographic roughness length”. More consistent terminology would make the paper more accessible to non-specialists.
R D Larter
3rd October 2025
-
AC1: 'Reply on RC2', Dorothée Vallot, 15 Dec 2025
We would like to thank the reviewer for its pertinent comments and we reply to them in the text (review comments in bold and reply in plain).
I consider this to be an interesting study addressing a fundamental problem in cryospheric science, the parameterisation of ice shelf basal melting. Having said that, I am not a numerical modeller or mathematician and think it is essential that the manuscript is also reviewed by someone with expertise in these fields (I have not yet looked at the other review I was notified about because I want to provide completely independent comments).
Recent studies have revealed that ice shelf basal melting is more spatially and temporally variable than previously appreciated, so new ways of modelling the melting, particularly for areas where ice shelf basal topography is complex, need to be explored. Therefore, I am pleased to see this account of such a study.
In general, the manuscript appears to me to clearly articulate results of a thorough exploration of different approaches to parameterising basal melt. I have no major comments, but I have annotated a copy of the manuscript with many minor ones. I have also included a few of the more important of these comments below.
We agree with all comments in the annotated file and are correcting the manuscript accordingly.
A general comment is that annotations on several of the figures are in a font that is too small if they are going to be published at the same size as in the review copy. Colour scale bars in several figures are annotated with values expressed with exponents, whereas I think converting these to more simple numbers would make them easier to understand for a greater number of readers.
Good remark that is taken into account.
Labels need to be added to identify the different parts of figures 2, 3 and 4 that are referred to in the captions.
True. These will be added in the next version.
After studying Fig. 2a and the related text it is still not completely clear to me where on the boundaries of cells various parameters are being computed. I infer that the blue box on the left of the figure is depicting a plan view of a cell since it is labelled with x and y axes, but does the position of the T-point that is depicted mean that is on the upper or lower face of the cell, or both, or in the middle of the cell? The V-point is described in the text as being “at the center of the upper and lower grid face”, but if the blue box labelled with x and y axes is indeed a plan view, then the V-point appears to be on one of the lateral faces. I think it would be helpful if the authors could try to further clarify their explanation of where these different points lie by modifying the text or the figure, or both.
The other reviewer also had a comment on this figure and a new version we hope more understandable is in the supplementary.
We use now a 3D version of the cell. The T point is actually in the middle of the cell. For u, we can use the east and west side of the cell and for v, we can use the north and south side of the cell to be clearer.
Fig. 2b is the first place in the paper that maps of the different ice shelves appear. It would be useful to include a regional map in the figure, or in an earlier figure, showing where the areas represented by these local maps lie. A similar map to the one in Figure 4b could be used for this purpose. It would also be helpful to add distance scale bars to each of the local maps to make clear that they are at different scales.
The detail I am most puzzled about is the apparent near random distribution of first wet cell heights shown in Fig. 2b. Fig. 2a and Fig. 3a indicate that first wet cell height increases with depth of the ice shelf base. Therefore, I would expect to see a trend towards more purple and black cells, indicating greater cell height, nearer the grounding lines. Unless an error has been made in the production of the figure, I think a further explanation of the apparent near random distribution of first wet cell heights is needed.
This is because we use the z*-coordinate system with partial steps. In the paragraph 2.1.1, it is explained in the first sentence: “We use a curvilinear z*-coordinate system with partial steps to adjust the thickness of grid cells adjacent to the sea floor (index kb ) or ice shelf draft (index kt ).” The thickness of the cell depends on the level position but also where the boundary is. This is what we try to show in figure 2a.
We can add an explanation in the figure caption: “First wet cell half thickness below ice shelves in the ASE configuration with 121 layers: Pine Island, Thwaites, Dotson-Crosson and Getz from left to right. The cell thickness is adjusted to the ice shelf draft through partial step.
I am also puzzled about the difference between the values of the constants given in the “roughness length” column in Table 1 and the values of the constant α given in the caption to Figure 5. Perhaps I have misunderstood something, but it seems to me that these constant values for a particular experiment should be the same in both places.
This is a mistake from our side. It is from an older version where the alpha value was wrong.
There are a couple of places where variable terminology is used to describe the same thing, e.g. “forcing” and “driving”, and “hydrographic roughness” and “hydrographic roughness length”. More consistent terminology would make the paper more accessible to non-specialists.
Good point. We use hydrographic roughness length everywhere.
R D Larter
3rd October 2025
-
AC1: 'Reply on RC2', Dorothée Vallot, 15 Dec 2025
Peer review completion
Journal article(s) based on this preprint
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 2,329 | 370 | 37 | 2,736 | 72 | 85 |
- HTML: 2,329
- PDF: 370
- XML: 37
- Total: 2,736
- BibTeX: 72
- EndNote: 85
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Dorothée Vallot
Nicolas C. Jourdain
Pierre Mathiot
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(26980 KB) - Metadata XML
This paper presents a modification to the way that melting of glacial ice shelves is implemented in the NEMO ocean model, with the intention of capturing the effects of variable ice base roughness on melting. This is certainly an important aim, since current melting parameterisations assume that the ice base is uniformly rough, which is known to be inaccurate. The current models also assume that the boundary layer is evenly resolved by the model, which is also inaccurate. The recommended approach aims to address these points by implementing a logarithmic boundary layer speed profile into the model.
The paper concludes that i) the new implementation is less sensitive to changes in vertical grid resolution in the model, and ii) variable spatial roughness is an important control over ice shelf melting patterns, and iii) any temporal change in roughness will cause a related change in melting.
My personal opinion is that the new melting implementation may be better than the current NEMO approach, which is based on an ad-hoc boundary layer thickness. Also, it is clearly desirable that the stress formulations used in the melting parameterisation and momentum solver are consistent with each other. However, this is quite subjective and arguable in several areas, and I do have some substantial points to make about the current manuscript, requiring major revisions. Ideally, these would require new simulations to address, though I am reluctant to suggest this. I think this work deserves to be published, ultimately. I think this new approach raises some discussion points that will be very valuable to the community. I would like the authors to consider the points below and I am happy for them to disagree with my view points, as long as the discussion is reflected in the paper.
Top-level comments
1) The main advance of the paper is in proposing the use of a formulation of the drag coefficient that is dependent upon the grid cell height (equation 3). However, this function is only very weakly dependent upon this height until the grid cell is below 2 metres or even 1 metre (see figure attached as PDF). This seems to be implying, for me, that for most purposes the assumption of a spatially constant drag coefficient should be OK. I feel that this key point should be mentioned. However, and more importantly for the paper, it leads me to wonder where the difference between the REF and LOG simulations comes from. At the moment I think that these differences are probably arising because REF uses the Losch layer of 30m to drive melting, while LOG uses the top wet cell. I propose that they are not caused by the change in the drag coefficient. Clearly, this distinction is central to the paper’s conclusions. I am reluctant to suggest any new simulations, but I think the only way to determine this is to run simulations with the melting parameterisation driven by the top wet cell, not the Losch layer, but using a constant drag coefficient (for both 75 and 121 levels). If these simulations look like REF then the paper’s conclusions stand – the change in the drag coefficient formulation is important. If these simulations look like LOG then the paper’s conclusions should be redirected towards studying the numerical stencil of the melting scheme, rather than the drag formulation. At the moment I think the latter point might be true.
2) There is a conceptual issue around this application of the law of the wall formulation (equation 3). This is designed for constant-stress boundary layers adjacent to an interface that is hydraulically smooth. It can cope with a variable roughness length, as used here, but only as long as that topographic roughness is not larger than the viscous boundary layer depth, such that it would start to disturb the boundary layer with lee eddies, etc. This is apparent in the paper, because the roughness lengths used are tiny - of order 0.1 mm. However, the topographic features present beneath ice shelves and being invoked here are much larger than this – e.g. order 10m in the case of terraces, and order 100m in the case of rifts, ice blocks, and basal channels. This very large topography will completely disturb the boundary layer and so should not be expected to be well represented as a sub-millimetre scale perturbation. Extreme topography on this scale would not be expected to smoothly increase the stress transfer through the boundary layer. This kind of roughness could not even be simply characterised as ‘form drag’, let alone the viscous drag assumed here. The Hughes paper cited by the authors is a nice illustration of this. This is certainly the roughness that the authors are invoking, as it is only these very large scale/large amplitude features that will be present in the surface DEMs being used here. Given this uncertainty, I am not sure what to recommend. It is possible that a constant drag coefficient is even a better choice, in the face of the extreme uncertainty around this problem. I fully support the authors’ goal of incorporating roughness into melting formulations, but the application of this law is not clear. I don’t know how the authors might like to proceed. The best approach here probably depends on whether the drag law makes a difference (point 1 above).
3) I think the logic of the new drag formulation is this: When the top cell is thicker than the log layer it needs to be parameterised, but when the cell is thinner, it will start to resolve the log layer, and hence the parameterisation should be reduced. But, does the model resolve the log layer when the cells become thin? I am not sure we know this. I think it depends upon the vertical viscosity parameterisation in the model, and the implementation of the momentum boundary conditions. In real fluids the log layer profile arises because the eddy viscosity increases linearly with distance from the interface. I don’t believe this eddy viscosity parameterisation is usually used in ocean models. There is also an important corollary here: when the vertical resolution is varied, is the vertical viscosity varied as well? With a higher resolution it would be typical to reduce the viscosity. That would lead to a smaller boundary layer, which is harder for the model to resolve. So it is not clear that a finer grid resolution will better resolve the log layer, as assumed in this new parameterisation.
4) The reduced sensitivity to vertical grid resolution changes is a key feature of the new approach. However, I wonder if the higher sensitivity of the old approach is due to the way in which the ‘Losch’ layer is implemented in NEMO. In the Losch (2008) paper, this approach was only used for temperature and salinity, not velocities. More importantly, the Losch layer thickness was equal to the top grid cell thickness, not the fixed 30 m that is used here. I wonder if the 22% sensitivity of the old approach here is caused because NEMO uses a fixed 30m Losch layer thickness? Clearly this is going to mean that as the cell thickness increases, fewer cells will be in the Losch layer. But with the original Losch approach, which I think is used in other models, there will always be 2 cells in the Losch layer. Given that the vertical viscosity may change as the grid resolution is changed (?), those 2 cells may contain a fixed fraction of the changing boundary layer thickness. Thus it could be speculated that the original Losch approach would be less grid dependent than the approach apparently used in NEMO. This implies that the new parameterisation may not be as beneficial in models that follow the original Losch approach. Gwyther et al (2020) document substantial variety in the way the models approach this. I think the only way to address this would be with additional simulations following the standard Losch approach, though again I am reluctant to suggest further simulations.
5) A key conclusion is that increasing the roughness increases the melting, which could mean there is a feedback between increased ice shelf damage and increased ocean melting. I like this result, but isn’t the top-level result already obvious? If we increase the drag coefficient, melting will increase, as shown by previous studies, and so this conclusion should be stated in relation to this existing knowledge. Where this study could offer a new finding is in the precise sensitivity of melting to changes in roughness – i.e. the curve of melting vs z0. The authors could consider plotting a curve of drag coefficient versus z0 to show how melting might change as the roughness changes? I assume it is highly nonlinear in some way, which would be an interesting finding. However, there is of course the high uncertainty in whether this is the right drag law to use (see above), so perhaps the authors may feel this is not useful. There are also lots of other considerations in terms of how melting might change, such as how the velocity changes with roughness.
Smaller points
General: With partial cells, do the authors apply a minimum cell thickness? I believe that is standard, to avoid numerical instability as the cells get very thin. If, for example, a cell cannot go below 1m thickness, this will reduce the influence of the new drag law even further (see above)
Line 95: ‘upper and lower’ is confusing
Line 100: why the subscript 3 ?
Figure 2a: I’m afraid I found this really confusing because the cells are all in different places. I think this is because the z* coordinate means that all of the different levels are at different heights and that is randomly applied here, though I am not sure. Could it possibly be redrawn for z coordinates, since that expresses the key difference here (Losch vs 1st wet cell), and then just note that actually z* is used? Also, the cells next to the ice only have 1 open horizontal face (right), with their left boundary being ice. I think it is far more common that they would have two open horizontal faces, with ice only on the top, so maybe better to draw that? This emphasises the need to average over the 4 u and v values on the cell edges.
Figure 3a: Similarly, I found this really hard to interpret. Maybe select the depth range of ice shelves and plot cells on a linear vertical axis? Maybe that is not possible. Give a table with cell heights at chosen depths (10m, 100m, 500m, 10000m)? Also it looks like the 121 levels were deliberately chosen to represent ice shelf cavities better – is that true? i.e. they are not a uniform increase in resolution, relative to 75 cells.
Line 127: I’m afraid I didn’t understand how this roughness measure was calculated, and this is important. I assume that in the NEMO cell there are many DEM cells. How are all of these cells compared to arrive at R?
Line 134: The mention of 2 m resolution here is confusing.
Line 170/Table 1: I was confused about this roughness length, because it gives a drag coefficient much higher than 0.0025 for all grid cells (see figure above) and yet it gives the same melting overall. I think this confirms that it is not the changed drag coefficient formulation that is important, it is the switch to using the top cell velocity.
Line 221: this explanation does not seem to work for Thwaites?
Paragraph 1 on page 14: I’m sorry, I was not able to follow the 3 points of explanation offered in this paragraph.
Figures 9 and 10 and sections 3.1 and 3.2: I could not find a physical explanation of why the melt rates were different in these cases. There are interesting patterns of difference, but these were not explained in terms of the underlying velocities in the different regions, and how well they were represented by the different drag laws, and the different vertical resolutions.
Equation (B4): bottom equation has an extra term C, which I don’t follow. Is this a constant background velocity?