GENERAL COMMENTS:

The manuscript deals with an important and debated problem, the permafrost carbon feedback. More specifically, it uses a numerical simulation approach for assessing solute transport in a permafrost catchment, in order to study the impacts of climate change on DOC export from the soil of the considered watershed. The mechanistic modelling of dissolved species in permaftrost affected areas is at the forefront of the current research in cryohydrogeology, and the used tool is a state of the art one. Thus this work is both novel and of great interest for cryospheric sciences.

Meanwhile, critical weaknesses affect its reliability. First of all, the study focuses on a specific site, while the claimed aim of the work is to draw general conclusions about lateral transfers in permafrost regions. This is contradictory, and if the authors want to produce a general study then the problems of upscaling and transposing to different biogeoclimatic permafrost contexts must be carefully dealt with. I have also a number of technical concerns about the set up of the performed numerical simulations, for instance the poorly justified pseudo-3D approach, the absence of any convergence study, and some problems in the definition of boundary conditions. Finally the discussion of the obtained simulation results is rather weak and should be strengthened. All these points are detailed in the specific comments below.

Overall I think that the manuscript cannot be published in its present form. A thorough effort is needed for making it more solid and meaningful, including possibly new computations and results depending on the answers to some of my concerns (e.g.: numerical convergence, relevance of top and bottom boundary conditions). So I recommend that a major revision of this work should be undertaken prior to reconsider whether or not it should be published in The Cryosphere.

SPECIFIC COMMENTS:

- Title and l 8-9: The title refers to ‘a permafrost catchment’, while abstract draws conclusions for ‘permafrost regions’. There is no a priori reason that conclusions made from the study of a peculiar catchment should be relevant for permafrost regions as a whole.

- l 17: Why considering mineralization of organic matter to CO₂ only, and not also to CH₄? May be due to the peculiarity of the study site. This should be explained.

- l 27-41: Other questions than PCF maybe related to lateral transfer of dissolved chemical species from the active layer under climate change (e.g.: impacts on river / ocean ecosystems).

- l 71-72: The assumption that DOC may be considered as a non reactive tracer should be discussed here, with relevant references.

- l 80-82: Most likely this behavior will be site dependent! For instance not the same for a tundra hill slope with low evapotranspiration or a boreal forest hill slope with high evapotranspiration, depending on the slope itself, etc.

- l 111: Is Qs sink/source term accounting for ionic exclusion? To say, does the used model take into account the specificity of solute transport with freeze/thaw?

- l 113: Why neglecting dispersivity? The errors associated with this simplifying assumption should be discussed.

- l 118-119, « path of highest flow accumulation »: What about lateral transfers from the slopes to the thalweg? For studying the drainage of active layer waters and the associated solute fluxes the choice of focusing on the thalweg line does not seems obvious to me. This should be discussed.

- l 121, « pseudo-3D approach »: I cannot understand the usefulness of varying y-axis width of cells along the x direction. According to Fig.2, there is only one mesh along y-axis at every x position. So this seems to me a 2D mesh (as said in l 169). The Figure 3 and the associated explanations are hard to follow and should be thoroughly improved (e.g.: a large blue area appear on Fig.2, but only green and red area are mentioned in the text.) Whatever the pseudo-3D approach really is, the statement made by the authors that it « allows us to account for thermal and hydrological balances across the entire catchment area without the need for a complex and computationally intensive full 3D mesh » should be justified. How does the y-axis width varies along the thalweg? Has any comparative study between results of a full 3D approach and results of this pseudo-3D approach been done? It should have been, prior to use the pseudo 3D approach, or at least if it is not practicable due to computation time the approximation should be discussed, as well as the associated errors.

- l 126-127, « main area of interest »: Why focusing on such a tiny area of 20 m length, while simulations are done on a 1040 m large domain? This important choice should justified and discussed.

- l 129-130, « By directing precipitation from the upper slope to the lower areas, we ensure realistic hydrological conditions with flow accumulation towards the valley bottom. »: What is meant here? The physical equations solved by the numerical model do ensure that gravity exerts a vertical descendant driving force on water, so that water flows from top to bottom when gravity dominates. This sentence seems pointless, I recommend to delete it.

- l 130-131,« This division of the mesh allows for accurate modeling of the thermal-hydrological processes in the catchment. » : Such a statement should be justified.

- l 132, « Each column in this mesh area [...] varies in width in the y-direction. »: Why and how? See above the point on the pseudo-3D mesh approach.

- l 132 – 141: Here mesh cells dimensions are described, but not justified. In a modelling study relying on PDE spatio-temporal discretizations (e.g.: finite differences method, finite volume methods, finite elements methods, etc), it is mandatory to assess the truncation errors by dedicated convergence studies, designed for the simulation case under concern. The results of such a convergence study should be given here, by means of a upper bound of the truncation errors for the outputs of interest. Only in this way one can be sure that the variations discussed in the numerical results are not simply due to truncation errors. Please include here the results of such convergence study for the case under consideration.

- l 134-137: I do not understand why a buffer zone is needed. This should be explained.

- l 142-146: All of these sentences look like unjustified statements. Moreover what is exactly stated is not completely clear. For instance what means « preserving the subsurface volume representation of the catchment » and « a natural equilibrium without artificially imposed boundary conditions »? This paragraph should be either deleted or rewritten in a clearer and justified way.

- l 147-148: « The vertical sides of the model are assigned zero-flux boundaries for water » this is questionable, especially for the outlet vertical boundary. Are there any field observation that can be used to justify this choice? This should be discussed.

- l 151, « in line with borehole observation in Svalbard »: This is important, since it is likely the reason why a 40 m thickness as been chosen for the modeling domain. Please add a figure with the mentioned soil temperature profile evolution, as well as a discussion for explaining in which way these data where used for choosing not only the bottom thermal boundary conditions but also the domain thickness.

- l 159-165: I have serious concerns about the chosen methodology for building ‘present day-like’ precipitation. Why not using a day-of-the-year average like for other forcings? This is not justified. Stating that « the resulting rainfall distribution resembles the variability of natural rainfall throughout the year» is not enough for making the arbitrary, artificial precipitation forcing data set relevant. The use of arbitrary precipitation data may impair the interpretability of the simulation results, so either the ‘resemblence’ between the artificial data set and the observation should be quantitatively demonstrated, either the observation data set itself should be used as forcing data.

- l 165-166, « Soil physical properties are defined to resemble highly conductive material. »: This should be justified. Why not moderately conductive material?

- l 169-170, ‘to establish a water table at target depth »: I do not understand. What is the ‘target depth’? Are there observational evidences of a water table ‘at target depth’?

- l 173-174: Why 10 years? Please provide the criterium used for this choice.

- l 176-178: The targets of this study are simulations under climate change, while the use of a thermal boundary conditions of constant temperature equal to present day temperature at the bottom of the domain (l 150-151) is not compatible with the simulation of climate change scenarios. Or at least, it implies the assumption that the temperature at 40 m depth is not impacted by surface temperature variations at the considered time scale of 50 years. This is a major concern for the validity of the produced simulation results. In order to demonstrate that this strong approximation does not impair the discussed results, the time variations of temperature in depth, close to the bottom boundary, must be shown. If they are not negligible, then the simulations should be re-ran with an appropriate bottom boundary condition (e.g.: geothermal heat flux).

- l 184-187: Additional explanations should be added here: why this rate and total amount of injection, why this moment?

- l 204-205, «air temperature (T avg ), which is the only variable that changes over time in the respective scenarios. »: Why? Precipitation does also change along time in IPCC climate change scenario. This choice should be discussed.

- Figure 4: Why the BTCs are multimodal? For instance, 2 modes for TOL carbon in Fig4.b and 2 modes for buried carbon in Fig4.c?

- l 271-272: Is this temporal partitioning between surface and subsurface transport in agreement with field observation? At least has the ponded water before mid-June has been observed in the field? Whenever it is possible numerical results should be discussed at the light of field observations.

- l 279-281: In Figure 4.d the peak of the 13 August seems higher than the one of the 14 June, although the opposite is stated in the text. This bimodality should be better explained.

- Figure 8: Vertical peak on 30th of August, dual peaks in the 1st of September … I think these features are strange. May be due to convergence problems? The convergence study must be done for assessing it.

- l 323-324: « A substantial fraction of the initially injected tracer mass (∼ 40%) moves vertically (both upwards and downwards) within the same mesh column in which it was injected (see Fig. 9a and d and Fig. 10). » What phenomena are responsible for this vertical redistribution ? Diffusion, freeze/thaw cycles related effects? This should be explicated.

- Figure 10.h: Numerical instability ? Should be corrected, or explained if it is not an artifact.

- l 346 : « This observation highlights the importance of mesh resolution in lateral transport simulations. ». I fully agree. A convergence study must be done.

- l 352, « vertical mobilization »: once again the involved mechanism must be explicated.

- l 407-411 : « We partly address this by representing a converging slope model setup, where the cell width in transverse direction varies depending on the distance in longitudinal direction. This way, the surface area of the catchment is preserved, and it is possible to accurately represent water and energy balances as well as infiltration and evaporation rates throughout the catchment. This approach has previously been applied by Gao and Coon (2022). » I do not think that using a one cell-thick discretization in the transversal direction may allow to simulate the effect of the watershed geometry, either convergent or divergent (using the terminology of Gao and Coon 2022). Nor in the present manuscript or in Gao and Coon 2022 are presented arguments for supporting the validity/usefulness of such a ‘pseudo-3D’ meshing methodology. A proper comparative study should be done for this, between results obtained with« pseudo-3D » meshes and with full 3D meshes. Of course with only one cell no lateral fluxes may be computed.The only interest I would see would be to weight the inward fluxes through the top cell face according to the area of the cell, but then why not simply apply a spatial weighting on the incoming fluxes prescribed by the boundary conditions? Including this in the meshing seems to me inappropriate and confusing. Anyway in this case the methodology used for computing the cells widths must be explicated.

- l 465, « Under the simulated environmental and soil physical conditions in this study, »: This precision hold true not only for this point, but for all the listed conclusions. The writing of the manuscript should better reflects the fact that this is a numerical study of transport in a specific site, with no possibility of automatic generalization for permafrost regions as a whole.

TECHNICAL CORRECTIONS :

- l 25: Missing ) at the end of the line.

- l 29, l30, and elsewhere: I would recommend to systematically use ‘organic carbon’ instead of just ‘carbon’ for naming the C part of the organic matter stored in permafrost.

- l 32: Missing s at ‘question’.

- l71-84: Part of this paragraph should be in the Methods section (e.g.: choice of distinguishing four carbon pools and using labelled tracer for identifying them).

- l 151, «bottom horizontal boundary »: According to Figure 3, the bottom boundary is not horizontal.

- l 251-260: These information should be included in the Methods section, along with a figure for quantitative locations of the injection points and of the measurement points within the modelling domain.

- l 258, « explicit », « continuous in space » : Odd vocabulary. A BTC represents the temporal evolution (not exactly continuous, since there is a time discretization) of concentration at a given location in space, while a plume vizualisation represents the spatial distribution of concentration at a given moment.

- Legend of table 2: The use of different concentration thresholds for TOL and buried carbon should be mentioned in the body of the text, in the Methods section.

- l 273: « when all runoff is occurring in the subsurface »: Then it is not run off, but groundwater flow.

- l 285: Do not refer specifically to a Figure in Supplementary material (Fig. S3) in the body of the text. Instead, refer to the supplementary material as a whole.

- l 286 – 287 : « However, given the specific solute transport patterns in this model, » What is meant here ? Unclear.

- l 299, « vertically transported upward »: Oddly said. Strictly speaking, the topographical effect mentioned here does not generate ascendant flow.

- l 302: Do not refer specifically to a Figure in Supplementary material (Fig. S4) in the body of the text.

- Caption of Fig. 6, « Note that the tracer mass is restricted to the uppermost subsurface cell in this snapshot and is difficult to visualize in this illustration » : True, then this figure has to be improved. May be that plotting the two variables separately would be an option? Besides, Fig. 6 is mentioned only very briefly once in the text. Either it should be deleted or more extensively commented.

- l 315, « some tracer is moved upwards »: see comment on l 299.

- l 317: 12 % does not obviously look negligible.

- Figure 7: I can’t see anything regarding the solute mass. This Figure should be thoroughly reworked so that it becomes informative. Besides, what is the « ponded depth » ?

- Figure 9: The annual cycles/peaks should be discussed. « The visual increase in mass above 100 mol in the injection columns in (a) and (d) is not a physical phenomenon but a result of aggregating and rounding errors during post-processing of the model output. » Then the post-processing should be improved.

- l 350-355: If the post-processing method does not allow to conclude, then it should be improved.

- Caption of Figure 11: « The visual increase in mass above 100 mol in the injection column in (b) is not a physical phenomenon but a result of aggregating and rounding errors during post-processing of the model output. » Once again the post-processing method should be improved.

- l 414-420: this should be part of the introduction, not the discussion.

- l 420-425: this should be part of the Methods section, not of the Discussion.

- l 432-447: This should be part of a perspective section, not of the Discussion.

- Section 4 Discussion: Given the parts that should not be included in the discussion (see the three comments above), the Discussion section is rather short (a bit more than one page), and lack of in depth analyses of the produced results. I recommend to strengthen this section, including for instance elements for linking the simulated behavior and the site specificity.

- l 452, « mechanical transport »: Sounds weird. Passive transport would be more relevant.

Major comments

The manuscript describe a study using a numerical model to investigate the factors that influence mobilization and transport of carbon in soils in permafrost regions. Warming and permafrost thaw are likely to cause the transport of carbon from soils through aquatic systems, so a better understanding of physical mechanisms is crucial to the development of more accurate models and climate simulations.

The study is designed as a highly idealized case for a small area in Svalbard. The grid setup is not standard. The time periods noted in the text are variable and confusing. The short time period of the abrupt that complicates interpretation. This application of a single model in this fashion makes extension of the study findings to broad regions of the Arctic questionable. The results would be more meaningful by applying the model for at least one or two other configurations based on observations of soil texture and organic carbon content, weather data, elevation gradients, which are available at select sites in Europe, Siberia, Canada, and northern Alaska. Doing so would allow the authors to have more confidence that what they are seeing represents a robust response and not an artifact of a unique model setup. The authors should also state the choice of application to the Endalen valley on Svalbard where observations are sparse. Is there a rich history of field study there?

The 0.25 km2 catchment is relatively small. Why not use a traditional uniform horizontal grid and with vertical soil layers? How many grid cells would be needed at resolution of 1 m? Would the computational expense be prohibitive? How does the setup result in “a natural equilibrium without artificially imposed boundary conditions”, and what does that phrase mean? Boundary conditions and zero-flux boundaries in a typical 3-D model setup could be the same as those in this study. In other words, I don’t see the advantages of this current grid mesh, nor the implications for the interpretation of results. While the authors refer to computational expense, more justification should be presented.

The wide variation in time steps mentioned is awkward and makes interpretation of the influence on results difficult. There is a rate for tracer injection of per second and a time step of monthly for the warming experiments. If it were implicitly daily, could mineralization in the future simulations like the present day simulation? The time step of per second is confusing in light of a monthly step in the future simulations.

The term seasonal variability in the title is odd and unnecessary, as seasonal variability in a study like this is essential. Having the term ‘sensitivity study’ in the title would be more meaningful.

Minor comments

Line 6: grammar, “Here, we analyze of solute transport…”

Line 43: Unclear. What is “they” in the statement “With permafrost acting as a largely impermeable layer between the two, they are…”

Figure 3 caption: It appears from the graphic that width (y-direction) decreases with distance in the x-direction. As x increases from 0 to 1000, the width gets smaller. Please clarify.

Line 159 states that the model is run for an average yearly cycle. Then at line 172 there is reference to year-to-year differences. How can year-to-year differences in spinup be determined? Further clarify the spinup time period and protocol.

Line 171: typo “a the”

Line 173: What ten years? Is weather data not for an average yearly cycle (climatology)? Be specific about the ten year period.

Line 186: Clarify how numerically does a tracer become injected at a rate with a time step of a second. What is the model time step? If it is daily, how can something occur every second? Does the mass enter as a total each day? Intrinsic time step of the model is not clear.

Line 235: typo, for for

Figure S3: component mass overlay is difficult to see in the graphic. Also add label for elevation for the y axis to make clearer what is being shown.

Line 18 in supplement: units of %C ¿ 10 are unclear.

Line 337: typo, “have lead to”

Line 357: Do the magnitudes (+25 and +60 cm) of the abrupt active layer deepening have any grounding in reality? That is, do any studies report abrupt thaw that is some meaningful fraction of the amount modeled here? Some perspective on the magnitudes would aid in interpretation of the results.

Line 368: “...is released from both the injection column…” grammar and unclear meaning.

Line 391: At the mention of ancient carbon, it would be helpful to remind the reader that this is frozen carbon. For readers that skim the abstract and discussion, the terms ancient and buried will have unclear meanings. Making clear that ancient means ancient currently frozen carbon would help.