It is a well written paper documenting a finite-element based approach to modeling sea ice dynamics in the framework of viscous plastic rheology. The approach is based on quadrilateral elements and explores both bilinear and bi-quadratic representation for sea ice velocity. The tracers are represented using discontinuous elements, but the authors find that the sensitivity to the accuracy of tracer advection is not very high. I recommend publishing this manuscript after minor revision.

I would recommend to briefly discuss the case of nonconforming (CR) element. According to earlier studies, it provides a very high resolution, and will similarly perform very close to the bi-quadratic case. The computational load is in this case seemingly lower. Also, I would suggest to add numbers allowing one to judge about the computational time required by bilinear and bi-quadratic cases. The first one largely corresponds to the approach taken by CICE, so the numbers will be helpful to judge on how the common type discretization is related to the  higher-order one. 

Minor points:     

line 14  'It is, for example, an important part ....' -- It contributes importantly to the global ...

15 circulation

62 'to limit diffusion'? Incremental remapping helps to ensure positivity. It relies on limiting, so it is not clear to what an extent it is limiting diffusion. 

83 It is the mean height (volume per unit area)

117 remove b

120 It is worthwhile to mention that this is never achieved in practice and mEVP commonly deals with non-converged solutions.

Section 3.1. The discussion of limitations due to advection is not really relevant. First, \Delta t is generally governed by the ocean model, where velocities are larger, so sea ice CFL will not be an issue. Second, there are internal stresses that depend on the mean thickness and concentration, and this may lead to a limitation

(for coupled equations) that is more demanding.

130 Sea ice was sometimes run with larger time step than ocean, but I am not aware about the situation described here.

133 Why this range? It will obviously depend on mesh resolution.

135 This is the limitations due to advection, it may matter for high-order in principle, but as I wrote, it is commonly ocean that matters most, and there is a wave-type limitation due to plastic response.

140 To be consistent ... -- MPAS, ICON and FESOM communities occupy a substantial part of climate modeling, and they rely on different meshes

141 unstructured --> distorted

144 for a depiction --> for an illustration         

Expressions (12) and (13): the lower indices start from 0 in (13), and 1 in (12)

186 Which order of RK is used? 

212 The Babushka-Brezzi condition is a subject of many publications and is well know, so at least include (see, e.g. Ern ...)  

Remark 1 -- Why 'Optimality...'?  Say before that the discretization on quads involves a lot of DoFs and computations, and a question arises on its optimality with respect to triangular meshes where functional spaces seem smaller.

256 by Gaussian

325 a triangle -- an element

Table 1  Ice concentration? 

450 is strongly affected ? I do not think it is the case. As soon as  stability is ensured, sensitivity is generally very moderate.

466 THis

5.2.3. I do not think this should go in the main text, move it to an appendix

Figure 14, caption  '10 momentum steps' --- 10 or 100?

5.3.1.   The efficiency on shared memory level is perhaps not surprising given the measures described. The real challenge, however, is the MPI parallelism, because of the large number of substeps in mEVP.  

Sergey Danilov

The authors present a finite-element model of sea-ice under elastic-visco-plastic approximation. Specifically, they chose a dG advection scheme for sea-ice concentration, and a mixed approach for the momentum equation, with continuous Galerkin method for all terms except the stress tensor, which was discretized with dG method. In the numerical tests, the authors focus on how different order cG and dG discretizations affect the sea-ice shear deformation, and study computational performance of their code parallelized using shared memory paradigm. Despite the increased computational cost, higher order discretization seem to provide superior results and much smaller diffusion. I am glad that authors consider moving to even higher basis polynomial orders, as the increased computational per-element intensity can lead to better parallel scalability results, offsetting the increased cost.

The paper is well written and contains sufficient math background, numerical tests and physics validation to merit publication. I do not have recommendation for revisions.