This is a thorough investigation of numerics for simulation of coupled processes in snow. As such, it is not an easy read and is very long. I wonder if it might be shortened beneficially. For example, the description of FEM will be unnecessary for those already familiar with it while not providing everything needed for understanding by those who are not. Otherwise, I have just reproduced editorial corrections that I noted in the manuscript; these will be of no interest for readers of the discussion, but should help the authors to improve the paper.

12

delete “to get”

33

“echoes on” is a curious choice of words

51

“too large relative to”

72

None of these references are to firn intercomparison projects. The statement is, I think, true for snow intercomparison projects. Mass and energy conservation was checked in earlier land surface model intercomparison projects (e.g. https://www.sciencedirect.com/science/article/pii/S0921818198000447); it was a surprise at the time that models did not all balance. Non-conservation and consequent drift is of more concern and more apparent in coupled climate models.

90

“is a an opening” – delete “a”

108

“the macroscopic system of equations is”

123

“dependency of the effective parameters on T”

129

https://doi.org/10.1029/2010EO450004

139

delete “which then writes”

141

“consists of a single prognostic equation for T”

158

“lies in the fact”

207

“we code every step of the method internally”

227

“that can be cast”

237

“the matrix form of the resulting algebraic system is”

250

“while every case … is first-order in time”

253

“combined with its stability”

264

“it is possible to diagnose the fluxes”

270

“One of the main goals”

271

“Details of the various approaches … to reach this conclusion”

285

delete “what regards”

334

“we refer to this as the proper form of the mass matrix”

348

“this method enables removing spurious oscillations”

383

delete “imperatively”

398

“obtained at the previous time step, is equivalent to performing”

399

“thus does not necessarily”

404

“without any reference to how these fluxes”

425

“this allows eliminating … from the continuity equation, which is then”

447

“the mass contained within each layer remains the same”

452

“would enable switching”

457

“consists of two numerical operations”

471

“without any explicit reference to how these fields should vary between nodes”

484

“an inconsistency that hampers mass conservation”

487

delete “which writes”

490

delete “comprised”

515

no need to abbreviate “respectively”

Figure 3

If the dashed blue “No Lump.; Prop.” line were plotted on top of the solid orange “Lump.; Prop.”, then both could be made visible (like “FEniCS” and “No Lump.; Improp.”).

543

“treated in detail”

556

“deviation in c”

565

delete “of”

613

delete “what regards”

625

“associated with stronger deposition rates”

629

“implies solving”

674

“equivalent to imposing”

685

“This oscillatory pattern propagates inwards over a number of nodes that grows in time.”

686

“but does not trigger oscillations”

708

“into an ice”

733

“associated with continuous sublimation”

737

“the dependence of k_eff and D_eff on ice volume fraction”

744

delete “of”

753

delete “of”, twice

719

“This is a generic result”

772

“might require considerably decreasing the time step”

778

“coming up with”

780

“This required knowing”

783

“associated with”

824

“the reader’s attention to”

825

“Even when one claims to be adopting”

826

“always requires discretization of the domain”

828

“must necessarily be made”

829

“This is equivalent to saying that there must”

833

“with a robust numerical treatment”

868

“enables closing”

971

A link to https://zenodo.org/record/5588308 would be useful here in addition to reference to Simson et al. (2021).

979

“consists of”

Minor comments/typos:

1. #90 – a an -> an
2. #251, 253, 256 - Cranck -> Crank
3. #526 - Full-implicit -> fully implicit
4. #261 – l.h.s -> l.h.s.
5. #263 – r.h.s -> r.h.s.