This paper describes a model of ice melange. It employs the Hibler 1979 sea-ice rheology and modifies it by introducing tensile strength and by selectively enhancing strength for grid cells labeled as icebergs. The authors then compare convergence of Picard, Newton and modified Newton solvers, concluding that modified Newton solver has the best performance. 

Major comments summary:

As far as I can tell the only new thing in this paper is the implementation of spatially variable strength within the Hibler viscous-plastic rheology.

The representation of icebergs as thick and compact pieces of sea ice and the concept of spatial variable strength depending on sea ice or iceberg label is identical to that of Vankova and Holland. The only differences are the advective scheme and the shape of the yield curve in the implementation. 

The yield curve the authors use is identical to Konig and Holland.

The authors emphasize throughout the paper the implementation of the strength modifications within viscous-plastic rheology, but it is important to realize that the Flato model, where the same was implemented is also viscous-plastic in its 2-D form.

In terms of the advection scheme and labeling icebergs, the authors do not show whether their scheme is able to handle separation of two icebergs that had been joined by convergent winds. I think this is a crucial feature to show. As of now there have been no iceberg-iceberg interactions tested and its unclear how the scheme handles the presence of multiple iceberg portions within a single grid cell. The authors also do not show (although it is claimed they do) that the scheme handles shearing - in test case 3 (which claims to test that) wind shear is applied away from the icebergs, the wind over icebergs is uniform and so the ability to rotate instead of deforming has not been shown.

Finally, I think the second selling point, after Hibler rheology, is the improved computational efficiency. The conclusion of the authors is that modified Newton leads to 40% improvement over Picard. However, there have been much more efficient solvers developed for the same system of equations. See for example Lemieux et al 2012, they show that Jacobian-Free Newton–Krylov (JFNK) solver is 3–7 times faster than the Picard solver for the same system of equations. The thing that makes the equations less likely to converge to a solution is whenever there is large stress buildup, which is what happens over the thick and compact bits of melange as material strength scales with thickness per the parameterization (that is typically an issue for thick and compact sea ice cover in continuum sea-ice only models as well, but the presence of thicker icebergs this is made worse). The test case on which the authors calculate the convergence properties is not one that experiences large stresses - it is icebergs moving through very thin and not very compact sea ice, so they are essentially freely drifting. In summary, I don't think there has been made any improvement towards computationally efficient tractability of melange as a continuum.

Other comments are below but there are many more in the attached pdf:

Good chunk of the abstract sounds like copy pasted abstract from Vankova and Holland 2017. It would be good to separate out clearly what it is that has been done before and what is the new development and lay that out unambiguously. Arguably the new thing is that the existing sea-ice/iceberg prototype model (continuum model with spatially variable rheology and icebergs represented as thick and compact pieces of sea ice held together by large shear and tensile strength) was implemented within the much more frequently used Hibler rheology instead of the original Flato cavitating fluid one.

The references are bit out of date - there have been some more new observational and modeling papers out that should probably be mentioned in this paper (I made a few notes in the pdf, but there is more out there). Some of the provided references are inaccurate (e.g. observational reference for a modeling claim). I didn't get to check all, but a few random ones were wrong so please check and fix all references.

There is an imbalance between the description of the equations which is quite comprehensive although already present in the literature, and the numerics section which is really thin and barely touched on. We don't even learn what grid the authors use.

figs 2-4: arrows are not clear - it was those figures that made me thing the grid might be unstructured, but then I realized it is just arrows overlapping which makes it difficult to compare their relative size.

Fig6: why is 48th time step interesting? Show all time steps or some statistics from them so that we can see the spread. Here Picard actually does better than Newton in this case, it seems, which is not typical.

Conservation of shape rectangle vs cylinder - this is a big misconception. The preservation of the shape of a moving object has to do with the advection scheme properties, not with the shape of choice.

The manuscript "Modelling ice mélange based on the viscous-plastic sea-ice rheology" describes a new method for modeling icebergs within sea ice (collectively called "ice melange") in existing viscous-plastic sea ice formulations within climate models. The study shows that this method shows improved accuracy and performance over a prior effort to incorporate melange into climate models. This study should provide a path forward for improving iceberg representation within ice sheet models.

The study is technically sound, but below I outline some ways it could be improved in order to be of interest to the audience of The Cryosphere and be in more direct dialogue with the iceberg melange modeling literature.

Major concerns:

1. In its current form, the manuscript is a straightforward model description paper, which looks at how to get from existing sea ice models to this new model, and certain improvements in accuracy and computational perfomance by adding tensile stress and using a modified Newton Solver. At that level, this paper is more appropriate for a journal like Geoscientific Model Development than The Cryosphere. To be of interest to the broader readership of The Cryosphere, the authors could take this wonderful new tool and use it to answer certain oustanding questions about ice melange that current tools cannot answer. One is how ice melange responds differently to forcing than just sea ice. Why does it matter that you have bits of thick ice which can resist tensile stress embedded within sea ice? The answer would likely be from a ocean/atmosphere perspective, as the manuscript does not currently argue that this tool would necessarily be used in conjunction with ice sheet models.

2. The motivation for this work, as described in this manuscript, it to be able to efficiently include ice melange into climate models. However, the test cases here exclusively look at cases where the resolution is much, much finer (0.1 km) than it is in any global climate model (10-100 km). What do the author's anticipate being the specific use case for such a formulation? Should it be used in global simulations? Should it only be used for regional simulations at very high resolution? Is the indicator function Phi (Equation 20) strictly necessary to have a discrete division between grid boxes with sea ice and icebergs, or could one have a continuous function (phi) between 0 and 1 that can capture the scenario where you have a large climate model grid box with some icebergs in it (but not completely filled with them)?

3. It remains unclear how the results shown here with the new method specifically improve upon the method of Vankova. The comparisons made in section 3 are between a standard VP rheology versus a VP rheology with tensile strength, whereas Vankova uses an entirely different rheology with tensile strength. Section 3.1 makes the most direct comparison by looking at a similar setup, but there needs to be a specific argument made about why the inaccuracies in Vankova's method are problems (i.e. related to error in representing aspects of the physical climate system that we care about). In general, this comparison with Vankova is critical to arguing why this new method is necessary, and it should be explored in more detail throughout the results and discusison section.

Minor concerns:

Line 1: influence on sea-ice-ocean interactions. Ice melange currently is not represented in climate model as numerically efficienty realizations do not exist. This motivates our development

Line 11: Newton's method

Line 12: melange into climate models

Line 14: Fjord with marine-terminating glaciers can be filled with sea ice into...

Line 17: the ice melange at some glaciers disintegrates [melange is year-round at some glaicers]

Line 17: ...the seasonal cycle of glacier flow is in phase with the presence of ice melange, leading to the...

Line 22: First sentence of this paragraph can be deleted as its redundant

Line 24: through releasing their freshwater storage from icebergs.

Line 25: delete "also"

Line 28: Delete "This is the more the case as"

Line 29: and the challenges of taking in-situ measurements within densely packed melange or near claving fronts

Line 39: several studies have implemented a particle-based approach for representing icebergs into large-scale ocean models, including those by Alon Stern and Anders Damsgaard.

Line 41: to commonly used continuum sea-ice formulations.

Line 43: model using a continuum approach.

Line 45: Using such an approach

Line 46: cavitating

Line 50-53: Amundson and Burton (2018) don't do a direct comparison between observations and a VP rheology. They do conclude that continuum granular rheologies based on non-local or friction dependent on inertial number work well in describing granular flow through channels, but they don't directly compare to observations. My suggestion is to rewrite this sentence to make it a more accurate reflection of that paper. Even if there isn't observational support for using a VP rheology for ice melange, there are good modeling reasons, and for now, that is sufficient.

Equation (in general): typically div is written as "del dot" ($\nabla \cdot$ in LaTeX) since the term "divergence" is not the same in all languages, whereas "del dot" is more "universal" notation 

Line 75: describes the divergence of the two-dimensional...

Equation 2: I'm a bit confused by how forces and stresses are being mixed here since they do not have the same units. Either some notation is mixed up, or some constants are missing.

Equations 3-4: why does the melange velocity enter into equation 4, but not 3? (In the most general sense)

Line 88-89: this sentence is a bit odd. Ice melange is melted from below by the ocean, this is an important sink term. I think the point is that you are just considering idealized cases where melting is ignored. If so, then just state that.

Equation (6): It would help the exposition to explain what e and P are right after this equation

Line 92: expressed in terms of the principal components...

Line 105: It would help if equation 14 came before equations 12 and 13 as a way to explain that there are two regimes with a smooth transition between them before giving equations for the two regimes.

Line 107: what is meant by the "viscous closure"?

Line 118: strength lead to a...

Line 126: what is meant by "normal flow rule"?

Equation 19: Why is the indicator function Phi needed if you already have phi? Also, what is the meaning of phi? It seems like a concentration of icebergs. The reason for this question relates to major point #2 above. If you are allowed to have a continuous concentration of icebergs between o and 1, then the model could make sense at grid resolution compared to global climate model, where a grid box might have some icebergs in it (while not being fully filled with icebergs).

Equation 19: What is c and how is it set?

Line 141: Is this the same thing as operator splitting? If so, that seems like the more common terminology to me.

Line 151: prove not proof

Figure 2-4: it is somewhat confusing to me why P-T is the correct contour to plot here (as opposed to concentration), and also that there are maybe three different quantities plotted? (P-T, arrows and a black contour for the iceberg, but not clear what quantity that represents). There should probably be a legend for these plots and a better explanation of what is plotted in the caption. It is also unclear why the quiver arrows indicating melange velocity have such variable velocity.

Line 200: Vankova and Holland (2017) stated the need for efficient

Line 204: what about the third setup makes it "difficult"

Line 213: residual threshold for convergence

Line 218: It would be useful to have some kind of explanation as to why the modified Newton solver appears to work so much better than the others.

Line 222: icebergs

Line 231: Most climate models represents the

Line 232: Delete "An"

Line 233: we recommend implementation