| 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.
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.
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?