Dear referee, 

We thank you for your valuable comments about our manuscript. The modifications for taking into account these comments will be included in the updated manuscript to be submitted for the next steps of the review process. Besides, please find below first answers to the points you raised.

Best regards,

The authors

General comments

    1. Methods 

The methods section is fairly extensive and contains a lot of information about the site as well as the model. However, from subsection 2.3 onward, it also contains results that are better placed in the result section. It would be good if the authors could clearly separate methods and results to make it easier for the reader to follow.

The proposed re-organization will be applied.

Specific comments

    1. L53-54: It would be helpful to add a number and reference here as to how much of the permafrost area is covered by boreal forest. 

About 55 % of the permafrost affected area is covered by boreal forest (Stuenzi et al., 2021). This information will be added to the updtaed manuscript.

    2. L75-78: Here the objectives of the paper can be highlighted more clearly. Naming the actual quantities that are being investigated (soil temperature, moisture, active layer thickness, southward shift) can help forming specific objectives that then are answered in a very concise bullet point list in the conclusion section of the current paper. Linking these two will make it easier for the reader to understand the full picture.

The proposed improvement will be applied.

    3. L101: Also provide an average active layer thickness either for each, NAS and SAS, or overall. This can also be done after the NAS and SAS had been introduced in one of the following sentences. 

The proposed improvement will be applied.

    4. L103-105: The description of the catchment is a bit abstract. A figure with a map of the catchment(s) would be extremely helpful here. It could be part of Figure 1. 

The proposed improvement will be applied.

    5. L168-171: This sentence is very long and a bit confusing, please clarify.  

The sentence will be splitted and rewritten for the sake of clarity.

    6. L178-188: A plot with the forcing data described in this sentence would be a valuable addition to the methods. It can be part of the supplement.

This information will be added in the supplementary material.

    7. L200 and throughout the entire manuscript: “coldest” and “hottest” scenarios do not adequately describe the conditions in the SSP scenarios. Terms like “high-end forcing pathway” (for SSP5-8.5) and “low-forcing sustainable pathway” (for SSP1-2.6) are more adequate. 

The proposed terminology will be used in the updated manuscript.

    8. L207-209: I do not fully understand this sentence. What is meant by “the means haven been summed”? 

What we call « current climatic conditions » consist in a synthetical daily climatic forcing dataset for a virtual representative year, constructed from observations available between 1999 and 2014. A multi-annual average is made day by day for each meteorological variable (air temperature, precipitation) for all the years of the considered observation period to build this synthetic dataset. For instance, the precipitation of the first January of the synthetic year is the average of the precipitations of the first Januaries from 01/01/1999 to 01/01/2014. Then this current climatic conditions synthetic year is used for performing the initial spin-up of the permafrost simulations. To build the climate scenarios, we add to this synthetic year the trends of yearly averaged air temperature and precipitation changes from 2015 to 2100, trends obtained from CMIP6 scenarios available for East Siberian region (IPCC). In this way we derived the used climatic projections for air temperature and precipitation in Kulingdakan (Fig. 2). For instance, the air temperature signal for year 2050 for a given climatic scenario is the sum of the air temperature of the current climatic conditions synthetic year and of the offset computed from the considered climate scenario. Explanations will be added in supplementary material for clarification.  

    9. L215: Is the rate +1.9°C/100yr interpolated to a full 100 year period? From my understanding you are only looking at the time frame 2014-2100, which is not exactly 100 years. Please clarify. 

These rates are calculated for the 2014-2100 timeframe. For instance, the +1.9°C/100 yr corresponds to a  +1.65 °C/87 yr increase for the considered simulation period. We chose to rescale the rates to century time scale for the sake of readability. This will be explicited in the updated manuscript.

    10. L216: What is the global temperature increase rate? Please provide and cite this information.  

Global air temperature increases over the 21st century from Fan et al.,2020 are given below and will be added to the manuscripted. 

SSP1-2.6 : +1.18 °C/100 yr 

SSP2-4.5 : +3.22 °C/100 yr 

SSP3-7.0 : +5.50 °C/100 yr 

SSP5-8.5 : +7.20 °C/100 yr 

Fan, X., Duan, Q., Shen, C., Wu, Y., and Xing, C., Global surface air temperatures in CMIP6: historical performance and future changes, Environmental Research Letters, vol. 15, no. 10, IOP, 2020. doi:10.1088/1748-9326/abb051.

    11. L238 and throughout the methods section: It is not entirely clear to me with which temporal resolution you are working. The snow cover plots seem to be interpolated over the year using monthly values (Fig 3.) but for soil surface temperature you are providing annual means. Please clarify this in the method section.

For snow cover estimation and soil surface temperature, the solver works at daily timestep. The smoothness in the SWE curves presented in figure 3b is due to the smooth shape of atmospheric condition signals used as input for the model (the synthetic year of current climatic conditions, see answer to point 8).  However, for vizualisation purpose when addressing the century timescale, we use a rolling annual average of the displayed quantity. This provides the climatic trend obtained, while avoiding the display of intra-annual variability, which would make the figure unreadable on a century-scale plot. Beside, it should be mentioned that soil surface temperature and water fluxes are provided to the permaFoam solver with a daily frequency. However, the permaFoam solver makes use of adaptative timestepping, with here minimal timesteps of 1 second. OpenFOAM proposes natively the linear interpolation of the daily boundary conditions to the time step sequence required by the permafrost simulations.

    12. L286-287: This is still part of the methods section but you are introducing an important discussion point, which would be better placed in the discussion section

The proposed re-organization will be applied.

    13. Section 2.4: I am unsure about how useful this information is for the main text of the paper. It is certainly interesting but might be better placed in the supplementary information as it is distracting from the main text. Also, I am curious about the amount of energy that was required to perform this modeling using 1.8 million CPU hours. Is it possible to convert this computational effort into an energy equivalent, offering a rough estimate of the energy required for permafrost modeling and comparing it to a more relatable context? For instance: The energy consumed by our numerical modeling process is equivalent to the energy required to power a typical 4-person household for [X] days/months/years, assuming an average household energy consumption rate of [Y] kWh per day/month/year." Information like this will be increasingly important with computer codes increasing in complexity and ever higher performing supercomputers.  

The section regarding High Performance Computing requirements and methodologies will be moved to the Supplementary Material. Regarding the energy consumption associated with the performed simulations, we propose the following estimate. On IRENE-ROME supercomputer, the power consumption is estimated to 5,02734 W/core (personnal communication from the operating team). Thus the energy consumption of our 1.8 millions of CPU hours simulation campain could be roughly estimated to 9MWh. For comparing this consumption with the one of a typical 4-person household in the European Union, we propose to consider the final energy consumption in households (all end uses, including water heating, space heating and cooling, cooking and electrical appliance) available in the Eurostat database for the year 2021 (1584677 terajoules*, equal to 440188055 MWh). Then, this total consumption may be divided by the census population in 2021 (445649015 inhabitants**) to obtain an estimate of the average energy consumption per person and per year in the European Union (0,987 MWh/person/year). According to this estimate, the energy consumed by our numerical modeling process is equivalent to the energy required to power a typical 4-person household for about 27 months. This comparison will be added to the supplementary material.

*  https://doi.org/10.2908/NRG_D_HHQ

**  https://doi.org/10.2908/CENS_21AG

    14. Figure 7: Is it possible to add the current soil temperature profile to these plots? 

The current soil temperature profil will be added.

    15. L383-387: This feels like an afterthought. It would be helpful if the authors would introduce the idea of simulating 30 extra years to test for thermal equilibrium section to explain what it is used for. Then this part could also be shortened and becomes more clear. 

The proposed improvement will be applied.

    16. Section 3.3: The first part of section 3.3 and Figure 10 are hard to understand for me. The quantity of “total water content averaged over the first 2m” is very abstract and does not provide much information. Could the authors please either highlight which processes are depending on changes in the total water content (vegetation is mentioned here, but for this aspect Figure 11 would be enough) or potentially remove part of this section and Figure 10? Figure 11 provides more information at a glance and is more easily understandable. 

Figure 10 presents temporal evolution until 2100 of the mean annual total water content averaged on the first two meters of the soil – this thickness of the considered surface soil layer is chosen so that it encompasses all the area with water content evolution under the considered climate change scenarios (see vertical profiles in Figure 11). The two main points of this figure are the following :

- illustrating the evolutions of moisture content of surface soils in the slopes of Kulingdakan watershed, stable in NAS and slightly decreasing in SAS ;

- illustrating the evolutions in the partition of total soil water between liquid water and ice.

These information are important for heat and water transfers in the soil, due to the couplings and non-linearities between these transfers. For instance, decreasing total water content induce decreasing soil thermal inertia, while decreasing share of ice vs liquid water induce a decrease in apparent thermal conductivity. These information are also important for vegetation dynamics, since vegetation needs soil water uptakes for transpiration, and may only uptake liquid water. These elements will be added to section 3.3.

    17. L428-431: What is meant by “upward water movements” here? Figure 11 is very illustrative and important, but I cannot follow its description very well. Partly because both SAS and NAS are being mentioned in alternating sequence and partly because the processes leading towards the shift in soil moisture are not well described. Hence, the second part of section 3.3 can be extended and a bit more information can be provided on the processes causing the soil moisture gradients to change.

Figure 11 presents the comparison of vertical profiles of the annual mean of total water content under current climatic conditions and in 2100 for the four considered climate scenarios. The processes driving these evolutions of vertical moisture profiles are complex, that intertwin coupled and non-linear heat and water transfers as well as changing evapotranspiration fluxes. Then the comment of Figure 11 is largely descriptive: according to the simulations, we observe changes in the moisture gradients within the soil. Then these moisture gradients can be interpreted in terms of water movements according to the generalized Darcy’s law. The second part of section 3.3 will be rewritten so that the limit of the proposed interpretations appear in a clearer way.

    18. L467 and onwards: While evaporation is discussed briefly here, I think it deserves more attention as it is a crucial variable of future climate change in the Arctic (see e.g., Clark et al. 2023). Its potential effect on soil temperatures and runoff would be interesting to discuss in a boreal forest context.

For now in permaFoam, evapotranspiration is assumed to be solely constituted by transpiration, while the evaporation within the soil is neglected (Orgogozo et al., 2019). This assumption is made in the context of the study of boreal forest areas, in which transpiration largely dominates over evaporation in the hydrological budget (e.g., Park et al., 2021). Meanwhile, evaporation may be dominant in tundra environments (Clark et al., 2023), and may increase in the future in forested environement under climate change. Since considering soil evaporation would add another coupling between heat and water transfers through exchanges of latent heat, it could directly affect soil temperature evolution. These points will be the scope fo future modelling works, and a related preliminary discussion will be added to this manuscript.  

Clark, J.A., Tape, K.D. and Young‐Robertson, J.M., 2023. Quantifying evapotranspiration from dominant Arctic vegetation types using lysimeters. Ecohydrology, 16(1), p.e2484.

Park, H., Tanoue, M., Sugimoto, A., Ichiyanagi, K., Iwahana, G., & Hiyama, T. (2021). Quantitative separation of precipitation and permafrost waters used for evapotranspiration in a boreal forest: A numerical study using tracer model. Journal of Geophysical Research: Biogeosciences, 126, e2021JG006645. https://doi.org/10.1029/2021JG006645 

    19. Figure 13: Panel a is is very dense in information and an attempt is being made at describing what the gray shaded area and the black line are supposed to represent, but the methodology behind this approach is described very briefly and raises questions such as 1) is the 1°x1° polygon comparable in terms of vegetation type to the catchment in question? 2) Does it make sense to look at the entire time frame from 1997-2019 when looking at ALT? Active layer deepening over the last years is likely to skew the representation of the current state of the Active Layer. An average over the e.g., last 5 years might be more representative of the current state of the active layer. Since the dataset is only available until 2019, the years 2015-2019 might be suitable in this case.

First, it should be emphasized that in boreal forest environments active layer thickness can significantly vary even at scales way smaller than 1°: exposition, vegetation cover, soil profile, topography can affect the active layer. In the Kulingdakan catchment for instance, strong differences exist between the two slopes of the watershed, at the km scale. These spatial heterogeneities can be seen at the permafrost_cci product resolution (1 pixel = 926.63m).

Moreover, interannual variability also results in strong changes in ALT. In the current state of Fig 13(a), the black line describes the multi-annual (1997-2019) temporal average of the spatial average of the active layer thickness  over a 1°-1° polygon centered on a moving latitude ; the gray shaded area represents the min/max obtained for this spatial average during the considered period.  We can observe that the interannual variabilty can reach 40% of the average value. Therefore, in a given 1°-1° polygon at a given latitude, active layer thickness is strongly variable both spatially and temporally.  However, in order to propose a southward shift equivalent estimation based on a latitudinal variability, the 1°-1° polygon was considered big enough to smooth the small-scale inhomogeneities (~km) and small enough to capture the latidunal effect, including biome transitions (~hundreds of km) (see Anisimov et al., 2015 for vegetation zone). Regarding the timeframe considered (1997-2019), the choice of the maximum data available was made in order to reduce the interannual variability effect. Meanwhile, we agree with the referee that this timeframe does not represent the most accurately the current state of the active layer, in a context of a rapid climate change. An updated figure considering only the 2015-2019 timeframe will be proposed with the updated version of our manuscript, along with a more detailled description of the methodology used for producing it.

Anisimov, O. A., Zhiltcova, Y. L., and Razzhivin, V. Y., Predictive modeling of plant productivity in the Russian Arctic using satellite data, Izvestiya Atmospheric and Oceanic Physics, vol. 51, no. 9, Springer, pp. 1051–1059, 2015. doi:10.1134/S0001433815090042.

This contribution to the open discussion of our TCD manuscript ‘Future permafrost degradation under climate change in a headwater catchment of Central Siberia: quantitative assessment with a mechanistic modelling approach’ is constituted of two parts. The first part contains the answers to the referee 2. The second part presents a self-motivated correction of our manuscript.

FIRST PART: ANSWERS TO REFEREE 2

Dear referee,

Thank you for your comments on our manuscript. The changes to address these comments will be included in the updated manuscript to be submitted for the next steps of the review process. In addition, please find below initial responses to the points you have raised.

Best regards,

The authors

Major comments

1 - The study lacks a comparison of important thermos-hydrologic components against observations such as active layer thickness, evapotranspiration, and watershed runoff. The authors compared SWE and surface temperature against observations, but that is not sufficient to understand if the model can accurately simulate permafrost conditions, and importantly ALT.

Numerical modelling of permafrost conditions under current climate conditions in the study site has been presented in previous papers (Orgogozo et. al, 2019, Orgogozo et al., 2023). These simulation results has been obtained with permaFoam, the same cryohydrogeological solver used for building the centennial projections presented in the submitted manuscript. These results were in good agreement with the available observations of active layer thickness and soil temperature profiles. This point will be better highlighted in the updated manuscript.

2 - The authors focus on boreal forests, but vegetation dynamics are not included. This looks to me incomplete and inconsistent. Vegetations play a critical role in regulating ALT and ET. Is this due to model limitations?

The role of vegetation in regulating ALT and ET is taken into account in permaFoam through a mechanistic approach, by computing a sink term in the root layer accounting for tree water uptake on the basis of the potential evapotranspiration and of the soil water content (Orgogozo, 2015, Orgogozo et al., 2019). The main conclusion of the previously published numerical study of permafrost conditions under current climate conditions is that the transpiration by vegetation is a key parameter for active layer dynamics in the study site (Orgogozo et al., 2019). Then vegetation is taken into account in the used modelling approach. What is not taken into account is the evolution of vegetation cover in response to climate change (e.g.: change in root layer thickness). Coupling a vegetation dynamics model with the cryohydrogeological model used here would allow to study the impact of the climate warming-induced changes of the vegetation cover on permafrost conditions. This is beyond the scope of the present study and will be the focus of future works. This point will be better highlighted in the updated manuscript.

3 - The abstract needs to be revised. In lines 24-25, the authors mention quantifying soil temperature, soil moisture, active layer thickness, and water fluxes, however, the abstract reports an analysis of active layer and evapotranspiration only. In a modeling paper, more explicit quantification should be provided about important thermo-hydrologic components.

The active layer thickness is a variable that strongly integrates heat and water transfers in permafrost affected soils. This is the reason why we put it forward in the report of results in the abstract. For keeping the abstract concise enough we will not include analysis for all the variable considered in the study, but we will rephrase it for highlighting the integrative nature of ALT.

4 - There are many grammatical mistakes in the manuscripts that require careful consideration.

A proofreading service will be requested to review the updated version of the manuscript.

5 - The domain depth is 10 m, in such shallow domains, the bottom boundary conditions could impact the results as the surface thermal signal could penetrate deeper than 10 m, especially during projection. How sensitive are these results to the bottom boundary conditions?

In this study we do not impose a fixed temperature at the bottom of the domain, but a fixed gradient equal to the geothermal flux. Thus this bottom boundary condition should not be surface conditions / climate change dependent. In order to illustrate this, the figure below shows the computed vertical profile of mean annual vertical temperature gradient in the middle of both slopes, for current climatic conditions and in 2100 for the four considered scenarios of climate change.

See Figure R2.1 in the pdf file associated with this reply.

One can see that as expected thermal gradients are much more important close to the surface than at 10 m depth. Besides, the thermal gradients vary only slightly below 5 m depth. Then assuming a fixed thermal gradient equal to the geothermal flux at the bottom of the modelling domain at 10 m depth is a reasonable approximation. It does not mean that temperature does not vary at 10 m depth. In order to illustrate the changes of temperature at the bottom of the domain we added plots for 10 m depth in Figure 6.

Minor comments

L18: thermos-hydric? This is not a common terminology in the Arctic hydrology literature.

The ‘thermo-hydric’ adjective has been changed into ‘thermal and hydrological’.

L27: I assume the plus sign in +46% and +29% indicates increases compared to the current climate. Please clarify it.

The interpretation of the sentence is correct. The sentence has been clarified.

L28: this line is confusing without proper context.

The idea here is to propose a ‘space for time’ illustrative approach: to which southward spatial shift in current climatic conditions would correspond the future increase of active layer thickness simulated under climate change projections? The sentence has been reformulated for making this idea clearer.

L43: please provide the area percent covered by boreal.

This information has been added.

L120: Fig1. Please add a DEM image of the watershed, if available.

We added to Fig1. a vizualisation of the DEM of the Kulingdakan watershed.

L133: please clarify what “freeze-thaw of pore water” means here. Does it not freeze/thaw dry soil?Around 0°C phase changes occur only for water, not for the solid part of the soil.

The permaFoam solver does take this into account in its description of variably saturated and variably frozen porous media as four phase (soil grains, ice, liquid water, air) porous media. That is the reason why we specified « freeze-thaw of pore water ».

L164-167: Since the authors are discussing the speed and efficiency of the OpenFOAM, it would make sense to include further details for the sake of completion, such as the physical dimensions of the domain that was discretized into 1 billion mesh cells and how is this related to the current study?

The sentence pointed out by the referee does not pretend to discuss the capabilities of OpenFOAM, but only of permaFoam, the permafrost dynamics solver developed within the framework of OpenFOAM. The assessment of the numerical capabilities of permaFoam is the main subject of a previously published paper, ‘Permafrost modelling with OpenFOAM®: New advancements of the permaFoam solver‘ (Orgogozo et al., 2023). The test cases used for assessing the numerical performance of permaFoam in this previous numerical study are not directly related to the study site of the present manuscript, and this has been specified in the updated version. - Where else permaFORM is used in the Arctic, more references The works using permaFoam are all listed in the bibliography: Orgogozo et al., 2019, Orgogozo et al., 2023. The study site of the present work is the main Arctic site for which permaFoam simulations have been published (Orgogozo et al., 2019, Orgogozo et al., 2023). Simulations with permaFoam have also been performed for the Syrdakh watershed in Eastern Siberia (Orgogozo et al., 2023).

L128: I would suggest splitting this subsection into two sections, adding, for instance, Model domain. Also, it would help to add a schematic of the model domain with boundary conditions, etc. For the 2D domain, adding a transect can help better understand what exactly the authors are simulating.

The proposed improvement will be applied. A detailed description of the simulation domains is available in Orgogozo et al. (2019).

L212: Fig 2. “Annual air temperature” should be “Mean annual air temperature” and “Precipitation” should be “Mean annual precipitation” for clarity.

The proposed improvement will be applied.

L260: provide some metrics on the plot or in the text as well. Also, please explain why the model failed to capture the dynamics/fluctuation in the observed data.

Metrics will be added to the text of the updated version. The related discussion in Supplementary Material A will be more explicitly referred to.

L275-276: please explain what does “observed increases in air temperature (per 100)” mean here. If this refers to CMIP6 data, that is not observed data.

In the sentence (L274-276) «These rates of increase […] are lower than the observed increases in air temperature (+1.9°C/100yr for SSP1-2.6 and +7.8°C/100 yr for SSP5-8.5)», the second part refers to the mean annual air temperature increase available in CMIP6 data. Then the word ‘observed’ will be removed and replaced by ‘projected’.

L290: Without more details about the topography and model domain, this section does not help in understanding the complexity. The authors used 525K cells, what is the resolution? And why was such a high resolution needed, could the results change with 525K/2 cells – twice coarser resolution?

The resolution is of 1.2 m laterally, and between 0.25 cm (top) and 16.5 cm (bottom) vertically, since we use a vertically graded mesh in order to save computation time. A convergence study has been made to estimate the resolutions of the spatial discretization required for these calculation. This study is mentioned L300 and is detailed in Supplementary Material B. 525K/4, 525K and 525Kx4 cells meshes have been evaluated. The conclusion of this convergence study is that the 525K cells mesh is needed to obtained results accurate enough for discussing climate change induced variations of Active Layer Thickness. This convergence study will be introduced earlier in the updated version, as part as section 2.2.

L390: How long did one scenario take – real clock time?

One scenario took between 20 days and 30 days to be completed.

L315: 2.5 km wide and 10 m thick… are these measurements of the 2D model domain that was divided into 512K cells? This is not that big domain, if I am understanding correctly, this is probably over-discretized. How was this discretization chosen?

These are the dimensions of the 2D numerical domain, representing one slope (either NAS or SAS) of the Kulingdakan watershed. The discretization is chosen according to the convergence study discussed above, with a convergence criterium based on active layer thickness. Three meshes were used for this convergence study : 1024x128, 2048x256 and 4096x512, labelled as "coarse", "medium", "large" respectively. The differences in computed active layer thickness between medium and large mesh cases were small, with maximum differences of 2.2% for NAS and 1.3% for SAS. These criterium information will be added to the updated supplementary material.

L370: Please provide relative change/increase for both slopes, which would be more meaningful when comparing the North slope (cold) vs the South slope (warm) permafrost.

Plots with the relative changes compared to the current value of ALT will be added for both slopes.

L390: (left figure) In general, the change is less than 2% in all scenarios, however, I would expect the red curve on top of the green curve, as green is a much colder scenario and should be closer to thermal equilibrium as is the case for the SAS (right figure). Please explain.

This complex behaviour is due to the fact that in its current state the study of thermal equilibrium in 2100 is based on further simulations based on the repetition of the 2100 year, thus encompassing inter-annual variability (e.g.: year 2100 significantly hotter and dryer than previous ones in scenario SSP1-2.6). We will update the thermal equilibrium study by considering 5 years average of the period 2096-2100 as final climatic conditions rather than the year 2100 solely.

L419: Figure 10. The water would only exist in the active layer which is shallower than the 2 m. It would be more meaningful to plot a time series of the water content in the active layer only. This water content would also depend on the evapotranspiration, and Figure 12 shows ET does not vary significantly in both cases.

Quantifying the liquid water available on the active layer thickness could lead to interpretation difficulties, since the active layer thickness changes over time. This motivated the choice to averaged liquid water and ice volumetric content over a constant depth superior to the maximum active layer thickness obtained in all the scenarios/slopes. In this way the temporal evolutions of averaged water/ice contents are defined alike and are thus comparable between the different climate scenarios. Figure 12 shows that actual evapotranspiration experiences an increase between current conditions and the 2100 state, quantified between +5 % and +29 % depending on the scenario and slope (Table S3-S4, Supplementary Material C).

SECOND PART: SELF-MOTIVATED CORRECTION

An error of parameterization has been identified for the simulations presented in our manuscript submitted to TC. A clipping operation that limited the amplitude of the simulated soil temperature variations between the minimum and maximum monthly values of top soil temperature under current climate was erroneously implemented. So we reran the whole set of simulations with a corrected parameterization, by removing the wrong clipping operation. The comparison between the ill projections (with amplitude limitation) and the correct ones (without amplitude limitation) are presented below.

Using a correct parameterization, i.e. without imposed limitation of the amplitude of simulated soil temperature variations, we obtained detectable but limited changes in the simulation results compared to those obtained with the ill-parameterized simulations included in the currently submitted manuscript. For instance, in the case of the scenario SSP5-8.5, the simulations with correct parameterization leads to an increase of 9.9 % of active layer thickness in NAS and +14.1 % of active layer thickness in SAS in 2100, compared to the simulations with ill parameterization. Meanwhile, no changes in the trends we discussed in our manuscript occur. We will replace all the wrong results and figures with the corrected ones in our updated manuscript. Once again, these modifications do not change the conclusions of our study.

Comparison between ill-results and correct results are available in the pdf file associated with this reply.