This paper describes the use of SPH for sea ice dynamics. In more detail, it implements VP rheology with elliptical yield curve into a SPH model. This is an interesting exercise and could lead to further work on using SPH on ice dynamics. The paper is worth publishing after some modifications. The comments by this reviewer are, mainly, related to the usefulness and applicability of the method: What is gained by using SPH when compared to FDM or DEM?

On general level:

Paper is very technical and it not easy to follow without a background in SPH. Is there a way to make it easier to read? Considering the readership of TC, effort to do this might increase the number of readers. Even if this reviewer is very familiar with numerical models, cannot go through all the equations of the paper. Authors could consider if such high level of detail needed here or could some parts rely on referencing earlier work? What is new in this description and what is from other sources?

The particle size in all simulations is of order of several kilometers. In addition, if the reviewer understands SPH correctly, all quantities in SPH become distributed over even larger area due to smoothing by kernel functions. Discussion on the following five issues in the paper is warranted:

(1)    Is your model able to describe discontinuities in the deforming ice field with higher accuracy than typical continuum models (both in the case of opening leads and formation of ridges)?

(2)    Is the resolution of your model higher than typical continuum models?

(3)    is the coarse resolution, or large particle size, due to computational burden?

(4)    Does it even make sense to decrease the particle size when VP rheology is used?

(5)    Does an individual particle in your simulation have physical meaning (do they, for example, describe ice floes – you do mention that particle collisions occur and affect your solution so the particles appear to have a physical meaning)?

Overall, do the authors consider their technique to be closer to continuum model or particle-based model?  Section 3.2 you show that you model follows a continuum solution. While this is what you appear to be aiming for, the example raises a question for the need of the approach presented. What is the advantage of using SPH in this case (or in other examples)? The authors should include a paragraph on this in the discussion; please emphasize what is gained by using SPH.

It would be beneficial for the reader if you would include information on time step lengths and simulation times into your paper so that the reader can estimate how efficient the suggested approach is.

More detailed comments:

L37-38: If ice is thought to behave like a granular material, then is there a reason to believe the emergent properties of sea ice would not depend on floe size? Please comment on this in the text.

EQ17 & 18: Maybe it is mentioned somewhere in the paper, but is it common to define ice thickness and concentration as independent parameters? Maybe this is a misunderstanding by the reviewer, but at a given point in your simulation domain these two parameters cannot be totally independent, but, for example, A=0 should imply h=0. Could the authors comment on this shortly in the paper?

L375-379: You mention particle size does not affect jamming in your simulations. This is not what one would expect for a granular media. Does this suggest that your approach is not capable to fully represent granular behavior of an ice field (if such exists)? Do the particles have a physical meaning in your simulations? Please elaborate.

L381-384 (also FIG 8): (1) Are the “tree-like” peak stress values in your simulations transient (you use word oscillating stresses) or have your reached somewhat of a steady state in your simulations? If latter, is figure 8 just showing stress waves bouncing around the ice field in your simulation domain, which does not seem physically correct? Please clarify. If this reviewer is correct, the actual arches in your simulations appear to be limited into one close by the outlet. Please comment if you would to see more arches within the deforming ice in full-scale or if you would use DEM?

L433: You again mention stress networks. Your approach adapts features from continuum models, which cannot present stress networks. Please elaborate clearly in the manuscript if you think your approach can present them reliably or not—and if yes, why does it do so if the underlying rheological model cannot present them.

Dear Oreste, Bruno, Jean-Francois and Mohammed,

First of all I want to thank you for submitting a well edited and relatively easy paper to read. I am happy to see the new methodology applied to the viscous-plastic model and believe the information presented will be of use to others contemplating using such a model for modeling situations where discontinuities exist in the ice pack.

Overall I think this is a useful paper and the new method to solve the model is well described. I particularly commend you for introducing the limitations of the method and discussing clearly where it can be used. My only concern is that you present only two idealized case studies to demonstrate the model works. The reproduction of the analytic solution for 1-D motion against a wall is a good test and a useful bench mark. Did you consider the range of test cases that you would need to do to demonstrate the model performance? An arching case is a classic example with a free ice edge that is a good test of the models ability to handle the discontinuity. Have you considered the work by Billy Ip and Hibler on the VP model representation of ice arching. The set up they use is different to yours, with a conical domain. They present the flow states involving arching with dimensionless numbers, and demonstrate the impact of yield curve shape on the flow through the channel. This might provide a framework for you to test your solution against.

References

Flato, G. M. (1993). A particle‐in‐cell sea‐ice model. Atmosphere-Ocean, 31(3), 339-358.

Hibler, W. D., Hutchings, J. K., & Ip, C. F. (2006). Sea-ice arching and multiple flow states of Arctic pack ice. Annals of Glaciology, 44, 339-344.

Ip, C. F. (1993). Numerical investigation of different rheologies on sea-ice dynamics. Dartmouth College.

Minor Suggestions

In the introduction you jump in sentence 2 stating general sea ice model architecture to the constitutive relation. For readers who are new to modeling it might help to include the information that this constitutive relation is one of the terms in the momentum balance, and how it is the continuity and momentum equations that are discretised. This is very basic, but helps guide new readers.

line 35: increase -> increased

line 45: Have you considered the work of Greg Flato under Bill Hibler? He presented a semi-lagrangian approach for solving the VP model. In his manuscript there is an example of a test case that could provide insight if used with your method. This is a free sea ice edge with a vortex forcing applied over it. In a sense you can think of this as an idealized ocean eddy at the ice edge. It was a good test case for showing how Flato's method reduced diffusion at the edge that was apparent in Hibler's solution.

line 69: throughout -> through

Please check that algorithm 1 (and the tables/figures) are all referenced in the text and in order. I note that figures 8 and 9 are referenced in the text out of order.

equations 33 and 35: The O(xx) terms are not explained in the text.

line 203: Should there be a space between Kappa and l? Or is l a subscript? Also, please check that you are not referring to two different variables with the same variable name. Kappa is used again in equation 41, with subscript n. Do you need to use the same symbol for these two different variables?

line 358: grammar is off in this sentence. I think it should be "The water drag also causes a longer  time to reach steady state, since the ice drift speed is slowed."

line 368:"than what is common" remove what.

equation A17 and A18: Why use a number 1 here in place of a variable symbol?