In this manuscript, the authors provide an extensive description of the ice-sheet model Nix, which is an open-source, 2D thermomechanical model that solves momentum balance equations, mass conservation, and temperature evolution within ice sheets. Nix incorporates a series of approximations of the Stokes equations and includes a moving coordinate system to precisely track the grounding line at discretization points. The authors utilize the well-established MISMIP experiment to test their model, particularly to assess grounding-line migration with an overdeepened bed geometry, and extend the exercise to new cases (oceanic and atmospheric forcing) where the ice thermal model is activated. They analyze the hysteresis in each setting and find that the thermomechanics contribute to reducing the hysteresis as a result of the internal feedback between ice temperature and viscosity.

In general, I found the paper to be relatively well-structured, well-written, and easy to follow, with extensive descriptions of the physics and numerics involved. Although I did not find anything new on the methodological side, I believe the outcomes of the modified MISMIP experiments are interesting and warrant further investigation. Therefore, my main concern pertains to the type of paper and its suitability for GMD.

The fact that Nix is a 2D model significantly limits its potential as a community code, in my opinion. Although the model is well-described, which is expected for GMD, all components (physical models, numerics) employ well-known techniques that are already implemented in many other community glacier/ice sheet models in 3D (e.g., PISM, Elmer/Ice, CESM, ISSM, Bisicles, to name a few). Without offering anything distinctly new, I question the rationale behind publishing a 'model development' paper on the Nix model.

On the other hand, the MISMIP experiments incorporating active thermo-mechanics in response to oceanic and atmospheric forcings seem novel to me (has this been done before?). However, the findings are not thoroughly discussed, particularly in terms of connecting with existing literature, and the implications are barely addressed.

In conclusion, I believe the originality of the manuscript resides more in its application than in model development. Therefore, I recommend shifting the focus of the paper to highlight the results more prominently, relegating the thorough yet unoriginal model description to an extensive appendix. By providing deeper discussions and a clearer emphasis on the results, a submission to a more applied journal, such as the Journal of Glaciology, The Cryosphere, or Frontiers, would be more fitting. I hope my comments are useful.

The authors present the new ice-sheet model Nix. Nix is a 2-D ice-sheet model, which in this case means one horizontal and one vertical dimension, which are discretised on a moving grid, using scaled coordinates to ensure the ice surface and the grounding line always lie exactly at the last vertical and horizontal grid points, respectively. The authors present results of several well-established benchmark experiments as calculated by Nix, showing that it reproduces analytical solutions well, thereby verifying the model numerics. They also present some novel results where they include thermomechanical coupling in these experiments (which typically use a fixed, uniform, temperature-independent flow factor instead), showing that this has a small but significant effect on the steady-state ice geometry, and thereby on the hysteresis between geometry and thermal forcing (both atmospheric and oceanic).

In general, I think the manuscript is well-written, using clear English and well-formatted figures. I am also generally in favour of both the development of new ice-sheet models, and the publication of dedicated model description papers for these new models. However, aside from a few small comments, I have two major concerns about the applicability of the model which I believe need to be addressed before publication.

Major comments

#1: Applicability of the 2-D set-up. I am not entirely convinced of the practical use of the “flowline+vertical” set-up presented here. I understand that such a model could help provide accurate benchmark solutions for future idealised-geometry experiments, but the fact that it can never be applied to a realistic ice sheet severely limits its range of possible applications. I therefore think the manuscript would benefit from a clear statement of the intended applications of the model. I am also unsure whether the term “ice-sheet model” is fitting for a model that cannot be used to study realistic ice sheets. Perhaps something like “idealised ice-sheet model” would be more fitting.

#2: Applicability of the moving grid. As has been pointed out already in the original MISMIP publication, moving grids are very difficult to implement in two horizontal dimensions. This strongly limits the applicability of any findings produced with this model. Again, it might produce more accurate results for idealised-geometry benchmark experiments than other models do (at a similar resolution), but it cannot be used to e.g. find improved ways to represent the grounding line which could be applied in other models.

These two points together make me unsure what the added value of this model is over using other existing ice-sheet models that have a wider range of applications. E.g., if I hire a PhD candidate to do some work on the relation between thermodynamics and grounding-line migration, why would I advise them to use Nix instead of, say, CISM, PISM, Elmer, or any of the dozen other available models where such processes can be studied in both schematic and realistic cases?

Minor comments

Line 2: “Nix is a 2D thermomechanical model…”. It took me a few pages of reading to realise that you meant one horizontal plus one vertical dimension, rather than two horizontal dimensions. Please explain the “flowline+vertical” meaning of 2-D in the abstract.

Line 4: “…and shallow-ice.” None of the experiments you present use the SIA. If the aim of the model is to study grounding-line dynamics, it seems unlikely that it will ever be used at all. If so, consider removing it from the description.

Line 6: “…including those of stochastic nature.” This option is not used in any of your experiments. Either add relevant experiments or remove this statement.

Line 6: “Nix has been verified for standard test problems” Please state here already that these are the MISMIP experiments.

Line 15-16: “…Nix combines rapid computational capabilities…” You have not shown any results concerning computational capabilities. Either add relevant experiments or remove this statement.

Line 17-18: “…offering a powerful tool for advancing our insight on ice dynamics and grounding line stability” See major comments.

Line 21-23: I think the ISMIP intercomparison papers really need to be referenced here.

Line 24: “…leading a number of authors to question their stability” Which authors?

Line 42: “…concluding that moving grid models are the most reliable…” Only two of the models in the first MISMIP paper used the flux condition (“Schoofing”) approach, and if I recall, none of them used sub-grid friction scaling, both approaches that have since become commonplace in large-scale ice-sheet models (as opposed to moving grids, which I don’t think any models use). Please discuss this.

Line 44: “…a two-dimensional free-floating shelf…” Here too, this definition of “two-dimensional” deviates from what people would typically imagine when reading this term.

Line 71-72: “…similar accuracy to the Blatter-Pattyn momentum equations…” While the leading error term might be second-order with respect to the aspect ratio epsilon in both approximations, the velocity solutions produced by the DIVA are quite less accurate than those of the Blatter-Pattyn (try running the ISMIP-HOM experiments at the different length scales and you’ll see a big difference when L < 20 km). Please include this nuance.

Line 76-77: “…emerges as a clear outlier in terms of both model performance and its representation of the ice-flow physics itself.” Robinson et al. (2022) did not include a Blatter-Pattyn solver in their comparison, so this statement is slightly misleading.

Line 81-82: “…, numerical simulations of these rapidly flowing bands are a well-known difficulty.” Please elaborate on what difficulties with simulating ice streams are well known.

Line 82: “Diverse approaches are found in the literature…” What approaches in what literature?

Line 88-89: “It is a common approach to reduce the number of horizontal dimensions to the main flow direction so as to minimize computing time.” This is only true for idealised-geometry experiments, not for realistic applications.

Line 91: “Unlike previous models…” The original MISMIP paper includes both a higher-order model and a full-Stokes model. Several existing 3-D ice-sheet models can use higher-order or full-Stokes dynamics with thermomechanical coupling (e.g. Elmer/ice, PISM, ISSM, UFEMISM). Please add some nuance here.

Line 101-102: “…for efficiency and extremely fast computing…” Define “extremely fast”.

Line 103: “…NetCDF and Eigen libraries” Please include references for these libraries. Also, do any of them make use of parallel computing?

Line 113: “Our system is thought to thermodynamically evolve in time…” A strange way to phrase this.

Line 128: “…extremely high spatial resolutions…” Define “extremely high”.

Eq. 11: The basal drag coefficient beta does not appear in any other equations. Is beta an input use to calculate c_b in Eqs. 9 and 10?

Eq. 18: the squiggly rho (ratio of densities) is difficult for me to distinguish from the regular rho (density). As you simply write out rho_w/rho in e.g. Eqs. 19, 21 and 22, consider doing so here as well and removing squiggly rho altogether.

Eq. 18: what does S stand for?

Eqs. 21-22: the melt rate M does not appear anywhere in the continuity equation. Also, T0 in these equations in my understanding indicates the pressure+salinity-corrected (and therefore depth-dependent) ocean freezing temperature, is this also depth-dependent in your model?

Line 263: “Oceanic melting beneath ice shelves…” Earlier, you explained that in the flowline case (which assumes no buttressing), the geometry of the shelf does not affect the flow of the grounded ice, and that therefore you do not need to model a shelf. How then can sub-shelf melt affect your model? Is this included as a negative mass balance term at the last grid point of your grounded domain, as a sort of frontal melt rate? If so, how is this derived from the sub-shelf melt rate?

Line 280-281: “…the grid points distribution yields higher resolution near the grounding line following a polynomial or an exponential law.” Please provide this law.

Line 283: “…moving grid models are presumably the best choice from a numerical perspective…” Please clarify that this does not apply to models with two horizontal dimensions.

Figure 2: It is unclear on which grid points you define velocities, and on which you define ice thicknesses/temperatures. Also, the caption states that “…, the position of the last horizontal point (r−1/2) explicitly tracks the grounding line L(t)”. Does that mean that grid point r is floating?

Line 321-322: “…and the surface mass balance S(x, t) (Eq. 7).” Ah, S is surface mass balance.

Section 6: at what horizontal and vertical resolutions do you run your simulations? What kind of error would you expect based on the convergence tests you mention (but do not show) later on?

Figure 3b: does this include both the advancing and retreating phases of the experiment?

Line 405: “…near the melting point (Fig. 4d).” This should be Fig 5d.

Line 408: “…minimum temperature of −23ºC (Fig. 4c).” This should be Fig. 5c.

Figure 5: panels A and B are not referenced in the text.

Line 424: “Figure 7 illustrates the high sensitivity to that stems from the heat exchange velocity parameter gamma.” I find it difficult to understand the goal of this experiment. If I understand correctly that you use the melt rate M as a sort of frontal melt/calving rate, then M is a scalar number which scales linearly with gamma. So if you were to put M on the horizontal axis, then the curves of the four experiments should overlap, correct?

Figure 6b: how much time is there between the steps in ocean temperature? It looks like about 30,000 years, is that enough for the temperature to reach a steady state? In my experience the ice geometry itself equilibrates quite a lot faster. What would you expect to see if you reduce the time between the steps, so you deliberately prevent the temperature from reaching equilibrium? This is where the added value of thermomechanical coupling would really appear.

Line 439: “…a sensitivity test to spatial resolution (not shown)…” Why do you not show this? I’m actually quite interested. The stress-free boundary condition to the Blatter-Pattyn approximation at the ice surface is very tricky to implement, and I’m curious to see what order of convergence you get.

Line 459: “Setting the vertical advection at the surface equal to the accumulation (0.3 m/yr) is a standard choice…” What do you mean by this? It is my understanding that most large-scale ice-sheet models these days explicitly solve the heat equation in three dimensions, with the vertical velocity that appears in the advection term resulting from vertically integrating the 3-D horizontal ice velocity field (i.e. conservation of mass for incompressible ice). This is not so much a choice as it is simply a direct consequence of basic physics.

Line 473: “…stochastic boundary conditions capability…” You have not shown anything relating to his.

Line 491-493: “More generally … and viscosity” What does this imply for models of realistic ice sheets? It is still not uncommon for models in the ISMIP6 ensemble to neglect evolving ice temperatures in future projections of the Antarctic ice sheet. Do your findings imply that this introduces a significant error/bias in the results of these models?

Appendix A: there are two Appendix A2’s, please fix this.

Line 510: “The position in the spatial coordinates is then given by…” These expressions are not correct when the grid is irregular (and delta_Sigma_i is indeed a function of i). Also, time is not a spatial coordinate.

Line 522: “We thus have a linear system of 6 × r × p unknowns…” Don’t you mean r*p unknowns, interrelated by a matrix with 6*r*p non-zero coefficients?

Appendix A1: Please also provide the discretisation scheme you used for the boundary conditions.

Appendix A5: adaptative = adaptive

Appendix B: Either include some experiments with stochastic forcing, or remove this text.