Answer to anonymous referee # 1:

We appreciate the reviewer’s thoughtful comments and the attention given to our manuscript. While several important points were raised, we respectfully disagree with some of the interpretations presented. Specifically, we do not recognize the characterization in the summary that our paper “concluded that incorporating the fate of harvested wood products would shift the results from sink to source for forest peatland carbon accounting.” In our abstract, we clearly state that the net carbon balance (NCB) is negative (i.e., more exports than uptake, indicating a carbon source) for all system boundaries considered, except at the end of the simulation for the system boundary defined as Ecosystem + harvested wood products (HWP). Thus, our findings are contrary to the statement provided in the reviewer’s summary.

Additionally, we highlight that even when the NCB is positive for one system boundary, the integrated carbon sink (ICS) remains negative, underscoring the benefits of using ICS over NCB for evaluating these systems. Previous studies (e.g., He et al., 2016, BG; Kasimir et al., 2018, GCB) have used typically the NCB.

Here we will address the different comments from the anonymous referee #1:

FIRST MAIN COMMENT:

“First, a full forest rotation approach has been undertaken earlier in He et al. 2016 Biogeosicences study for the same site, although use another model but the same conclusion was already made there. Moreover, Kasimir et al 2018 GCB further presented full rotational GHG balance (not only CO2 but also CH4, N2O) of their studied site for several land use scenarios, including final harvest of the forest. Thus, I would argue their findings on the system boundaries are already known. The advance of ForSAFE has a few interesting developments in it (e.g. simple dynamic volume for peat, although they did not discuss much about it), but put it in the community of peatland models, the novelty is rather minor.”

ANSWER:

We appreciate the reviewer’s comments regarding the novelty of our work and its relation to previous studies. While we acknowledge the significant contributions of He et al. (2016) and Kasimir et al. (2018), we believe our paper adds valuable insights to the ongoing discussion around northern peatland management, particularly in terms of system boundaries and model dynamics.

First, our study differs in the modeling approach and process representation. He et al. (2016) and Kasimir et al. (2018) noted that heterotrophic respiration decreased in the final 30 years of their simulated forest rotation, attributing this to a shallower water table and reduced peat “concentration” due to a static soil volume coupled with high decomposition flux. We address this issue by implementing a dynamic volume approach, which we believe better represents the coupling between soil organic matter decomposition and water table dynamics.

Second, the calibration strategies differ significantly. Both He et al. (2016) and Kasimir et al. (2018) utilized extensive site-specific calibration with approximately 30 parameters, which can achieve good fits to observational data but may do so for reasons that are not universally applicable. In contrast, our study emphasizes in mechanistic default parameterization, minimizing site-specific calibration to ensure broader applicability. For instance, we adopted the peat decomposition rate constant from long-standing models like ORCHIDEE (Qiu et al., 2019) and LPJ-GUESS (Chaudhary et al., 2017), as well as the spruce carbon assimilation from PnET model for (Aber et al., 1996). This approach allows us to reflect on assumptions used in widely accepted models and provides a foundation for testing the model under diverse conditions.

These methodological differences yield distinct results. For example, He et al. (2016) reported an average NEP of 217 gC/m² for 1980–2011, closely aligning with observed NEE during the dry year of 2008 (204 gC/m²). However, by using a mechanistic model and avoiding the need to calibrate our simulations, we observe that 2008 was an outlier in terms of carbon uptake. This elasticity in response to environmental drivers is what allows us to simulate future NNE, but come with the tradeoff of a weaker fit to observations on high temporal resolution. Interestingly, in He et al. (2016), spruce NPP peaked in 1972 despite rising atmospheric CO₂, temperature while in our model this better environmental conditions manage to modulate the effect of stand age. Further comparison is highly interesting but model outputs of He et al. (2016) are not publicly available. We plan to make our model outputs available for future comparisons

Third, our definition and analysis of system boundaries differ markedly from those of Kasimir et al. (2018). While they present final values for various boundaries, the analysis in Kasimir et al. (2018) lacks a temporal perspective on how these metrics evolve. Moreover, the grouping of fluxes and therefore their definitions of system boundaries differ from ours in several key aspects:

• Their system boundary defined as soil includes N₂O and do not include carbon loses via water exports, in our work we focus on carbon and include carbon loses via water.
• We do not consider any system boundary that is the defined exclusively by fluxes included in the NEE.
• They do not incorporate the slow decay of a big proportion of harvested wood products (HWP), which we explicitly include in our “Ecosystem + HWP” boundary. This last point is particularly relevant for understanding post-clearcut carbon dynamics and has not been incorporated in the previous analysis mentioned by the reviewer.

Additionally, our study simulates a second full forest rotation under climate change conditions—an aspect not explored in previous work. We agree that the historical management history of the site was not entirely conventional, yet we chose to reproduce for a second generation in order to isolate the effect of a changing climate. This addition provides a more realistic “business-as-usual” scenario, knowing that the forest management plans are not always follow to the letter, and enabling us to test the model’s performance under changing environmental conditions and assess the impacts on carbon dynamics. It also lays the groundwork for evaluating alternative land-use scenarios, such as peatland restoration through rewetting.

While we acknowledge that our model developments are incremental rather than revolutionary, the incorporation of a dynamic soil volume marks a significant improvement. In He et al. (2016), peat decomposition decreased over time despite rising temperatures, due to a static volume assumption. By allowing soil volume to change dynamically, we enhance the coupling between SOM decomposition and the water table, providing a more realistic representation of peatland processes.

We would like to take the opportunity to address the reviewer’s concern about novelty to highlight that the dynamic volume approach is not included in previous studies, we would like to include the simulated change in physical properties (height of the soil profile and bulk density) in an appendix and discuss them as part of the representativeness of the model comparing with published data.

SECOND MAIN COMMENT:

“Second, I have many concerns on their model evaluation, specifically: do three years of NEE data (one year under mature forest, two years after clear cutting) enough for constraining (or support) their full carbon budget analysis for the two-rotation periods? One-year of 2008 (note even a dry year thus not a normal climate year) NEE data to evaluate, or precisely speaking benchmark, 80 year (first rotation) of modelled C results, with a very detailed processed-based approach (multiple factors control the C flows) like ForSAFE.  Not saying many of its simulated C flows did not benchmark with any field data at all. Their model data comparison (Fig. 5) shows large derivates between the measured and model daily data, suggesting these controls of the CO2 exchange are poorly captured by current model structure. The uncertainties in these (at least partly) also reflected by their simulated striking increase (double) growth over the second rotation period compared to the first rotation period (Fig. 6c). So, what cause these mismatches? if ForSAFE can not adequately simulate the underlaying process control and how these controls respond to clear cutting, as current results suggest. How can we trust the long-term predictions made by the model. This is a crucial question for the authors to think it over.”

ANSWER:

We acknowledge the reviewer’s concerns regarding the evaluation of our model and the limitations of using short-term NEE data. Indeed, long-term carbon flux data from forested drained peatlands are scarce, with only one known site providing more than five years of eddy covariance (EC)-based NEE data. This scarcity of long-term measurements underscores the necessity of employing process-based models, as noted in our introduction, to explore carbon dynamics over extended periods.

We agree that the available NEE data is insufficient to constrain the full carbon budget across the two-rotation period. As explained in the previous point, our modelling approach does not require constraining the model through calibration, but rather to test the underlying mechanistic assumptions. We did not perform extensive calibration based on the NEE limited dataset. Instead, we dedicated a section of the discussion to comparing our model’s outputs with published fluxes from relevant literature, demonstrating that the model still produces reasonable results. We argue that the strength of our paper lies in how we analyze and group these outputs to compare the relative importance of carbon gain and loss over time, providing new insights into carbon dynamics across different system boundaries.

Although the data is limited, it covers contrasting periods, including post-clearcut conditions, which offer valuable insights into soil responses with minimal vegetation influence. While the model does not perfectly capture daily fluxes, it is crucial to highlight that we deliberately avoided extensive calibration. The model involves numerous parameters to describe complex interactions within the system, and excessive calibration would risk overfitting, potentially obscuring structural deficiencies.

For the purpose of the paper, achieving a perfect fit to daily flux data was not a prioritized objective. Fine-tuning parameters such as decomposition rates, carbon use efficiency, or initial bulk density might improve fit but could mask the need for structural improvements, such as incorporating explicit microbial pools. Despite the mismatches at the daily scale, the model captures seasonal and yearly fluxes reasonably well. It also reproduces the magnitude difference between periods with and without tree cover, as shown in the results.

Given the long-term focus of our analysis, our priority was to produce reasonable yearly carbon flux estimates rather than optimize short-term flux predictions through extensive parameter adjustments. It is also important to note that uncertainties in the climatic time series, derived from reconstructions based on nearby weather stations, may contribute to discrepancies in daily flux values while still capturing overall trends.

Finally, the observed increase in growth during the second rotation is primarily driven by changing climate conditions. Rainfall increased by 23%, temperature by 40%, and CO₂ concentration by 57% during the second rotation. The enhanced growth reflects the influence of elevated CO₂ and temperature on photosynthetic rates, particularly under conditions with sufficient water and highlighting the fact that nutrient availability will remain high given the nutrient richness of the site. This effect is modeled using the carbon assimilation response function from the PnET model (Ollinger et al., 2002) which compiled empirical data of photosynthesis response to elevated CO₂. We discuss these dynamics in lines 330–335 and 503–506. This scenario likely represents an upper limit for carbon fixation, providing valuable insights into system behavior under high carbon uptake conditions.

In order to address the reviewer comment, we will more explicitly clarify the rationale behind the lack of calibration and its implication for the performance of the model compared to observations in the methods section. We will clarify the we are modelling a high tree growth scenario in the method section describing more thoroughly how the climate variables affect photosynthesis in the model linking to the specific formulations describe in the supplementary material, we will explain that we aim for a high growth scenario and we will expand over the effect of the CO₂ modifier function effect in the discussion.

THIRD MAIN COMMENT:

“Third, the authors commented the first rotation period was non-conventional management, however, still design their second rotation periods with the same management. Why not have a simulation design with conventional management that additionally evaluate the sensitivity of the results to varying and more realistically forest management, because the chance of having a late 72% thinning and minor storm harvest two rotation in a row, is just way too low if not impossible. A model simulation with more realistic management would make more sense for forest management implications.”

ANSWER:

We acknowledge the reviewer’s concern regarding the use of non-conventional management practices for both rotations in our simulations. Our rationale for maintaining the same management approach was to facilitate a direct comparison between the two rotations. This approach allows us to better evaluate the model’s mechanistic representation of carbon and nutrient dynamics under similar conditions except for climate. Introducing a different management regime for the second rotation would have reduced the comparability between the two periods, making it more challenging to isolate and assess the effects of climate change and nutrient cycling dynamics.

That said, we recognize the importance of simulating more realistic management scenarios to enhance the practical applicability of our findings for forest management. Future studies could explore the sensitivity of the model to varying management practices, providing insights into how these systems might respond under more conventional regimes.

OTHER COMMENTS:

“Please also note there are important issues concerning your data file shared on Zenodo: 1) the measured GWL data only contains the first two years; 2) the measured soil temperature data at 0.05 and 0.15m differ with what reported in your paper and clearly are wrong.”

R: Zenodo file has been fixed and is attached to this response.

“Line 16 negative meaning uptake?”

R: Metric fluxes and directions are explained thoroughly in the method section and figure 2.0, but further clarification will be provided in the abstract.

“Line 37 gC”

R: We are using another convention to be consistent with the way variables are described in the supplementary where the way of adding more information to units is with subscript; for example, cubic meter of water is m³_w not m³W.

“Line 120 The site is agricultural used before tree planting how is the initial conditions for the soil consider those land use, does it make any difference or not. He et al 2016 clearly show how the initial conditions influence the results of the simulated GHG balance.”

R: The initial conditions of the soil reflect its historical use as agricultural land, showing a relatively high bulk density compared to pristine peat and a very low CN ratio. These characteristics are indicative of the organic matter mineralization that occurred during the site’s agricultural phase. The initial conditions used in the model were derived from on-site measurements of bulk density, ensuring that the land-use history of the site is appropriately accounted for in our simulations.

To address the reviewer’s comment, we will include an appendix containing a sensitivity analysis of carbon metrics, exploring the influence of variations in the initial CN ratio and bulk density on the results. This addition will provide a deeper understanding of how initial conditions impact the simulated greenhouse gas balance and enhance the robustness of our findings

Line 140 Site description for the clear-cut site will also help the reader to understand the relevance of the understory vegetation in the C exchange for the forest few years of clear cutting

R: We appreciate the suggestion and will include the following edited paragraph in the methods section to provide a detailed site description and highlight the relevance of understory vegetation in carbon exchange during the years following clear cutting.

Edited Paragraph:

We simulated two forest rotations from 1951 to 2088 at a drained afforested peatland located at Skogaryd Research Station (https://meta.fieldsites.se/resources/stations/Skogaryd) in southwest Sweden (58°23’N, 12°09’E). This site experiences a hemiboreal climate, with high nutrient content organic soil, substantial peat depth, and effective drainage, managed under conventional forestry practices. Originally a fen valley, the site was drained in the late 19th century for agriculture before being converted to forestry in 1951. Until clear cutting in 2019, the site was dominated by Norway spruce. The clear-cut area covered approximately 0.16 km², with logging debris left on most of the site (4/5). Norway spruce was replanted on 2/3 of the site following the clear cut. In 2023, a barrier was constructed in the main ditch to raise the water level in the northern third of the site. Visual inspections revealed that vegetation cover increased in the years following the clear cut. By 2022, significant portions of the area near the eddy covariance towers remained covered by logging residues, while grasses and sedges, particularly in areas without logging debris, reached heights of 90 cm during midsummer.

Line 115 section 2.2 site description. I generally lack the clear cutting description, is there any rewetting effects for the site, like blocking drainage ditches, do they replanting spruce trees or they left small trees to regrow? How is the hydrological state of the clear cutting site, also I am lacking the descriptions of the eddy tower data, since the measured EC data is not published before so a brief description of the set up and measure and process processing is of need.

R: We will add a paragraph to explain the EC setup:

Added Paragraph:

“The 2020-2021 high-frequency data needed for flux calculations was acquired with an ultrasonic anemometer (USA-1, METEK GmbH, Germany) and a LI-7200RS gas analyzer (LI-COR Biosciences, NE, USA) mounted at 2.15 m height above low vegetation. Data acquisition frequency was 10 Hz and half-hourly average fluxes of CO₂ was calculated with the EddyPro software, version 7.0.7 (LI-COR Biosciences, NE, USA) following ICOS methodology (Sabbatini et al., 2018) when applicable. In short, turbulent fluctuations from the mean were derived using block averaging, de-spiking and statistical tests were according to Vickers and Marth (1997). Spectral corrections according to Moncrieff et al. (2004) and Horst (1997) were applied. Final fluxes were flagged for quality according to Foken et al. (2004) and fluxes flaged “2” were removed from further analysis. A threshold was estimated according to Papale et al. (2006) and data below this threshold was removed. Gaps in the dataset was subsequently filled using the REddyProCWeb online tool (Wutzler et al., 2018) that also partitioned CO2 fluxes into GPP and Reco following the nighttime method described by Reichstein et al. (2005).”

Line 173 why only these two parameters calibrated, even though we know the limKsat quite well. Plus, how these parameter uncertainties contribute to the uncertainties in the short-term model-data mismatch and long-term upscaling results?

R: We calibrated only two parameters to avoid the clear risk of overfitting associated with adjusting a large number of model parameters. Our approach prioritizes ensuring that the main controls on carbon dynamics are not overly specific to the unique characteristics of this site, thereby supporting broader applicability of the model.

Regarding the comment on how well limKsat is known, we are unclear about its intended context. If the reviewer could provide further clarification, we would be happy to address this point in more detail.

Figure 2, why not show the simulated fluxes with the width of the arrows

R: The figure was designed to illustrate the conceptual differences in system boundaries, as it appears in the methods section before any results are introduced. Using arrow widths to represent simulated fluxes would shift the focus toward presenting results, which are discussed later in the manuscript. Keeping the figure schematic ensures clarity in explaining the methodology without prematurely referencing outcomes.

Figure 3, there are clearly more than measured WTD time series, need to add that in the figure captions

R: If the reviewer is referring to the inclusion of several water table depth (WTD) measurement locations in the figure, we will clarify this in the caption. The revised caption will specify that the figure includes data from multiple WTD locations to ensure accurate interpretation.

I do have many concerns on how ForSAFE-peat handles hydrology which is arguably the most important variable for peatlands development and regulates C dynamics over this transitional clear-cutting period. The hydrological changes (including ET, Runoff etc) over the clear-cutting period and compare to before would be interesting to test the model in detail.

R: We agree with the reviewer on the critical role of hydrology in peatland carbon dynamics, particularly during transitional periods like clear-cutting. Our modeling objective is to simulate soil saturation and temperature conditions that are reasonable enough to test the coupled carbon and nitrogen dynamics in the system. Among hydrological variables, the position of the water table is arguably the most critical for peatland processes (Evans et al., 2021, Nature). Our results indicate that the model successfully simulates observed water table dynamics at the site, including the rise in water table levels following clear-cutting in the absence of ditch maintenance.

We acknowledge that the exclusion of ground vegetation in the model may influence evapotranspiration values. However, the reasonable agreement in water table dynamics suggests that the model captures similar water balances, provided that precipitation data are accurate and lateral inflow is negligible. While some processes, such as swelling and shrinking, anisotropic conductivity, and pronounced hysteresis in water conductivity, are not included in ForSAFE-Peat or other models previously applied to this site, we believe that achieving reasonable water table simulations provides a valuable context for testing carbon dynamics.

To address the reviewer’s comment, we will include an appendix comparing the water balance during the five years before clear-cutting with the five years following clear-cutting. This analysis will demonstrate the model’s response to hydrological changes associated with clear-cutting disturbances.

Figure 5 what is cold and warm period, define it and also why not show the daily measured vs daily modelled numbers

R: We chose to focus on seasonal periods (warm and cold) rather than daily values to highlight that, despite daily discrepancies, the model successfully captures seasonality. This focus aligns with the time dimension and scope of the study, where capturing broader seasonal trends is more relevant to long-term carbon dynamics.

The definitions of warm and cold periods are provided in the text (lines 247–248) and have now been added to the figure caption for clarity.

Answer to anonymous referee # 2:

We appreciate the reviewer’s comments on the paper, we found them very constructive and useful to improve the manuscript. Herewith we address the comments and propose ways to adapt the manuscript where possible.

FIRST MAIN COMMENT:

“Throughout the manuscript, the authors write about "afforested drained peatlands". Indeed, their site is afforested (with agricultural history) and their results primarily tell about an afforested peatland. But almost all the more general parts of the paper (and majority of the references also) actually address more generally forested drained peatlands. Some of those peatlands have been afforested, but many have been forested even in their natural state before drainage. I suggest to generally use "forested drained peatland", and only use "afforested drained peatland" only when truly addressing afforested peatland. "Afforestation" is well defined in forest-related fields of science (and also in common English), and it means turning a previously permanently open site into a forest”

ANSWER:

We appreciate the reviewer's comment regarding the terminology used in the manuscript. We agree that distinguishing between "afforested" and "forested" drained peatlands will enhance clarity. To address this, we will make the following changes:

• Throughout the manuscript, we will use the term "forested drained peatland" when discussing peatlands more generally, and reserve the term "afforested drained peatland" specifically for instances directly related to our study site.
• In the introduction, we will clarify that tree cover in northern peatlands increases with continentality. This will include noting that while peatlands with significant tree cover are rare in the UK, spruce and pine mires are more common in Sweden and even more prevalent in Finland.
• We will also update the title of the manuscript from "Evaluation of long-term carbon dynamics in afforested drained peatlands: Insights from using the ForSAFE-Peat Model" to "Evaluation of long-term carbon dynamics in forested drained peatlands: Insights from using the ForSAFE-Peat Model."
• Additionally, we will add a comment in the discussion sub-chapter titled "On the representativeness of simulated carbon dynamics and the abiotic context" to clarify the differences between our afforested peatland study site and other drained forested peatlands, particularly those found in Finland that are referenced in the paper.

These revisions aim to ensure that our terminology accurately reflects the context of our study and that our readers can clearly differentiate between afforested and naturally forested drained peatlands.

SECOND MAIN COMMENT:

“The model quite nicely follows the measured C balance variables. This is quite remarkable! This is very hard to achieve without some calibration of C process related parameters (because there is a super huge C storage in peat soil, which tends to decompose either way to fast or way too slowly in the model). Yet, based on the manuscript, only a couple of hydrology-related parameters where calibrated. Please, clearly tell if any of C process-related parameters were calibrated or not. And if not, please advertise and discuss how the measured C balance was so well followed by the model!”

ANSWER:

We also appreciate the reviewer's comment on the model calibration regarding the carbon balance variables. Indeed, we have tried to avoid extensive calibration of many parameters due, as our goal was to test common model assumptions related to Norway Spruce dynamics and peatland decomposition that are embedded within the ForSAFE-Peat structure. To address this comment, we will expand on the parametrization of the model and clarify our calibration process:

• The calibration details are currently found in lines 172-175. We will add further explanation here to clarify why we only calibrated the rate at which water leaves the bottom soil layer and the fraction of wood that respires.
• We calibrate against ground water level observations and tree ring-derived biomass data available for the years 2008-2010. The rationale for this approach is to test the model outputs regarding carbon fluxes if forest stand dynamics and abiotic conditions (such as groundwater level (GWL) and soil temperature) are similar to those observed in the field. We will add this rationale and further explain how these calibrated parameters (fraction of wood that respires and hydraulic conductivity at the bottom soil) influence the results after lines 172-175.
• Additionally, we will include a note in the results section highlighting that the model successfully captures the transition from carbon sink to carbon source, consistent with field observations.
• In the discussion, we will also mention that under a reasonable water regime, the decomposition rate constant used in ForSAFE—similar to that in other models like Orchidee and LPJ—yields realistic peat losses. Meanwhile, the PnET default parametrization related to carbon assimilation, respiration and litterfall together with the calibrated fraction of live wood lead to reasonable tree biomass accumulation, resulting in an overall agreement with values reported in similar systems in the literature.

These additions will provide clarity on our modeling approach and explain why the measured carbon balance was well matched by the model without extensive parameter calibration of soil organic matter decomposition and soil gas diffusion.  

THIRD MAIN COMMENT:

“No sensitivity analysis is conducted nor sensitivity of results to model assumptions and parameter selection discussed. I am not saying that you should do some full-scale sensitivity analysis (might not be reasonable), but at least some discussion on how robust your results may be is needed! So please discuss sensitivity/robustness of the results.”

ANSWER:

We agree with the reviewer that a sensitivity analysis could make results more robust. In response, we propose the following:

• We will conduct a sensitivity analysis on the initial conditions of nitrogen content, as our focus is on conditions with high nitrogen content, but field measurements indicate spatial variability in the carbon to nitrogen ratio (C:N), ranging between 18 and 24. This variability motivates us to explore how different initial nitrogen contents impact model outcomes.
• To illustrate the response of the model to different water regimes, we will also run simulations under various future precipitation scenarios, adjusting precipitation levels by increasing and decreasing them by 20%. This will help demonstrate how changes in water availability affect the long-term carbon dynamics in the peatland. As in the first sensitivity analysis, also in this second one we are motivated by the high uncertainty in projected precipitation and potentially high impact of different precipitation regimes on the carbon cycle.
• Attached as a PDF is the figure with the runs associated to the sensitivity analysis.
• After the submission, we improved the default parametrization to address an issue discussed in line 446, that the model generated excessive wood litter. To rectify this, we updated wood decay rates based on literature values and increased the fraction of wood that respires. By using the same calibration parameter (fraction of wood that respires), we were able to simulate a consistent carbon balance for the tree, resulting in similar stand dynamics while achieving a more realistic distribution among carbon outflows. The overall results remain consistent with the original intent of this work, but with reduced woody litter generation. If encouraged to revise the manuscript, we propose to replace the default parameterization with this improved version.

These additions will provide clarity on our modeling approach, increase robustness of the results via sensitivity analyses, and explain why the measured carbon balance was well matched by the model without extensive parameter calibration.

FOURTH MAIN COMMENT:

“C loss from peat soil is a most central element of the drained peatland problem. You describe well the water table depth (the strongly controls C loss from peat), but the actual ditching (ditch depth, ditch spacing, both strongly impacting drainage both at the site and in the modelling world) is very poorly described. Also, how the ditch maintenance is carried out in the model? (indirectly you seem to say that ditch cleaning to 60 cm depth at stand regeneration, but please say clearly). Also how ditches develop, in addition to soil subsidence? They collapse, get filled with vegetation, etc. leading to faster decrease in ditch depth than what is caused by only soil subsidence). See e.g. these two studies on the subject: https://doi.org/10.46490/BF453; https://doi.org/10.14214/sf.10494, and better describe the ditch network, please!”

ANSWER:

We appreciate the reviewer's comment regarding the ditch network, as C loss from peat soil is a central aspect of drained peatlands. To address this, we will make several modifications:

• We will improve the description of the site in line 119 as follows and adding a figure with a picture of the site: "The site covers an area of 0.2 km², with a gentle slope characteristic of its history as a former fen valley. It was drained in the late 19th century for agriculture, then repurposed for forestry in 1951. The ditch network forms a grid-like pattern, with the main ditch running north to south for 0.8 km, draining into Lake Skottenesjön. Smaller parallel ditches are spaced at varying distances, as shown in Figure X."
• We will add a more thorough description of how the model considers drainage in the section "2.1 Model description" by modifying line 102 as follows: "In the model, ditch function is simulated by setting an initial drainage depth. Layers above this depth experience lateral outflow when water content exceeds field capacity, with outflow regulated by the layer’s hydraulic conductivity and width, as described in Zanchi et al. (2021b). The drainage depth adjusts dynamically with changes in the soil profile; when the soil profile height is reduced due to net losses of soil organic matter, the ditch depth is also reduced by the same magnitude."
• We will explain better the specifics of the simulation for the ditching by modifying the paragraph starting currently in line 149 within the section "2.1 Model description" as follows: "We set the initial ditch depth at 0.6 m, based on ditch depth estimations from previous work conducted at the site (Nystrom et al., 2016; He et al., 2016). We aimed at simulating common ditch network maintenance (DNM) practices. In reality, the ditch was not maintained after 2019 clear cut due to a rewetting experiment that began in 2022. Therefore, NEE observations for the year 2020 and 2021 happened in a post clear cut and no DNM context. To harmonize what happened in the field with our aim of representing conventional, we reset the ditch depth to 0.6 in our simulation starting in 2022. In the model formulation, lateral drainage is influenced by changes in ditch depth, which reflect alterations in soil profile depth and variations in hydraulic conductivity due to changes in the bulk density of layers susceptible to lateral drainage. In reality, ditch depth is affected not only by subsidence but also by infilling due to sedimentation, vegetation growth, and bank erosion (Hökkä et al., 2020). However, these processes are not yet accounted for in the model, which simplifies the representation of ditch dynamics."

FIFTH MAIN COMMENT:

“For the two C balance metrics, NCB has an unambigious interpretation and physical meaning (change in C storage). ICS on the other hand, although being mathematically solid (you can always integrate) and having a good goal (taking into account the residence time of C in the atmosphere) does not have an exact physical interpretation. This might not be a problem, if you didn't have better options. But you (excellent modellers) can very easily calculate radiative forcing, which in a physically sound way integrates C balance and residence time of C-gases in the atmosphere (your C-gases are virtually all CO2) into climate impact. And then you can look both at the time series of instantaneous radiative forcing and take mean of radiative forcing over study period. You have otherwise such an exemplary solid manuscript, that I think messing up with ICS that is "almost something" kind of spoils it.”

ANSWER:

We appreciate the reviewer's insightful comment regarding the use of the Integrated Carbon Stocks (ICS) metric. We recognize that the ICS may not be as intuitive as the Net Carbon Balance (NCB), and we acknowledge that explaining its importance clearly has been challenging.

We agree that calculating radiative forcing would provide a better understanding of climate impacts, and we are actively considering atmospheric-related metrics for a subsequent study where we will assess different land management alternatives. However, for this contribution, calculating radiative forcing would require an assessment of fluxes from an atmospheric perspective, which would involve altering the fluxes considered due to carbon leaching across the three system boundaries as well as harvested carbon in the ecosystem boundary. Additionally, it would necessitate changes in flux conventions, which could create confusion—for instance, a carbon uptake viewed from the ecosystem's perspective is a positive flux, whereas from the atmosphere's perspective, it would be negative.

We believe that NCB (gC/m²) is more closely related to radiative forcing (W/m²), while ICS (gCyr/m²) is akin to what some call cumulative radiative forcing (J/m²) (Murphy & Ravishankara, 2018; https://www.pnas.org/doi/epdf/10.1073/pnas.1813951115). In essence, if we consider NCB from the perspective of atmospheric carbon, it represents the change in atmospheric carbon content over time, accounting for atmospheric carbon dioxide decay functions and other gas concentrations. This approach would allow us to calculate radiative forcing, effectively translating changes in carbon content into changes in radiative forcing (which has units of power). Integrating radiative forcing over time yields units of energy (power × time)

A particular value of NCB metric can be achieved through various pathways, and ICS provides a straightforward means to compare those pathways from a climate change mitigation perspective. The importance of the time aspect in assessing the climate benefits of carbon sequestration in ecosystems has been highlighted previously (Sedjo & Sohngen, 2012; https://doi.org/10.1146/annurev-resource-083110-115941). The use of a mass × time metric has also been proposed for accounting for carbon permanence in carbon accounting problems (Fearnside et al., 2000; https://doi.org/10.1023/A:1009625122628). As expressed in Muñoz et al. (2024), studies like Sierra et al. (2021; https://bg.copernicus.org/articles/18/1029/2021/) have shown that ICS can effectively account for the time carbon spends stored in ecosystems, providing a more comprehensive means of analyzing and comparing trajectories of carbon accumulation, as demonstrated in the figure from Muñoz et al. (2024) Figure is found in https://onlinelibrary.wiley.com/doi/10.1111/gcb.17229 and reproduced under CC by 4.0 license. The figure was not modified.

Answer to anonymous referee # 3:

We appreciate the reviewer’s comments on the paper, we found them useful and will use them to improve the manuscript. We are also thankful for the small detailed comments on typos and similar, those ones are at some point hard to catch, so thanks.

Herewith we address the comments and propose ways to adapt the manuscript where possible.

FIRST MAIN COMMENT:

“However, there is no information or justification of applied model approach in the introduction. There are many models, which are used for this purpose, therefore authors should justify the selection of the approach applied in this study. It would be better to place the ForSAFE-peat model in the context of current modeling studies. The posed questions in the end of introduction can be answered depending on model features and complexity. If model is sufficiently complex and well parametrized, the answer on first question always will be yes.”

ANSWER:

We appreciate the reviewer’s comment highlighting the need for a justification of the modeling approach in the introduction. We agree that the paper would benefit from placing ForSAFE-Peat in the context of current modeling studies and providing a rationale for its selection. To address this, we have added a new paragraph in the introduction (line 72), which contextualizes ForSAFE-Peat among existing models and clarifies its features and relevance to this study. Additionally, we acknowledge the reviewer’s concern about the first posed question, given that a sufficiently complex and well-parameterized model might inherently yield a positive answer. However, we emphasize that our study does not perform extensive calibration. Therefore, this question remains important for evaluating the ability of the model to capture system dynamics under the given parameterization..  

Added paragraph in line 72:

Here we introduce the dynamic ecosystem model ForSAFE-Peat and use it to analyse long-term carbon dynamics in a forested drained peatland. Several models have been developed to represent carbon dynamics of coniferous forest and peat soils, ForSAFE-Peat integrates many common assumptions often use in these models and is our objective to reflect on their effects on the simulated carbon dynamics. ForSAFE-Peat simulates plant dynamics as a big leaf model where photosynthesis is a function of foliar nitrogen content using the same structure of the PnET model (Aber & Federer, 1992). This representation has been widely use to study managed coniferous forest in northern latitudes (Belyazid et al., 2011; Belyazid & Zanchi, 2019; de Bruijn et al., 2014; Gustafson et al., 2020). ForSAFE-Peat simulates soils as a group of layers that can expand or contract due to soil organic matter content changes similarly than in peat development models like HPM (Frolking et al., 2010). Soil organic matter is represented by several compartments, where the decomposition flux of compartments that represent litter fill a pool that represents peat, resembling approaches like the one in Yasso07 (Didion et al., 2014). This allows a simple representation of litter quality and peat. Decomposition is described with linear kinetics where peat decomposition rate constant is the same that has been used to evaluate future carbon dynamics of northern peatlands by land-surface models such as ORCHIDEE (Qiu et al., 2018) and LPJ-GUESS (Chaudhary et al., 2022). ForSAFE-Peat follows traditional assumptions, making it a suitable tool for exploring carbon dynamics in peatland systems and critically examining commonly used methods for their representation

In this study, we used the ForSAFE-Peat model to conduct a long-term simulation spanning two full forest rotations in a well-studied drained and afforested peatland in southwest Sweden, utilizing primarily pre-calibrated parameters. Model outputs were analysed to represent various system boundaries, and different metrics were applied to evaluate carbon exchanges across these boundaries. Consequently, we explore the following two questions in this contribution.

SECOND MAIN COMMENT:

“Concerning the applied forest model: it was tested only for forest site. In fact, more general approach for simulating C dynamics in drained (and not drained and restored) peatland would be preferable, where other land use options also considered like grassland and even arable land. The applicability of used model in broader sense should be at least discussed.”

ANSWER:

We appreciate the reviewer’s insightful comment regarding the broader applicability of the model to simulate carbon dynamics under various land use scenarios, including grassland and arable land, and the potential for its use in restored peatland systems. We agree that this reflection was missing and have addressed it in the revised manuscript. We acknowledge the need to discuss the general applicability of the ForSAFE-Peat model to broader ecological contexts. We are actively adapting the model to simulate other plant functional types, allowing it to represent different ecological states. Importantly, the soil dynamics associated with a restored peatland are expected to emerge naturally from the current model formulation. To reflect this, we have added a discussion at line 535

Added paragraph in line 72:

In general, the current model formulation of carbon dynamics, which is based on common representations embedded in other models, provides a reasonable platform for analysing peatland systems despite certain limitations. While this contribution focuses on a drained forested site, the model structure is flexible and applicable to other conditions, such as waterlogged soils (not drained) and natural vegetation, including grasses and mosses, provided appropriate parameterization of the vegetation submodel is implemented.

Specific comments:

“Reference list for the literature sources mentioned in the supplement with model description need to be provided.”

A: Indeed, it has been added.

“Figure 3a – it is not clear, which data are shown – observations from different locations, averaged observations or observations from one of the locations. Please extend the figure caption accordingly.”

A: Figure 3 and figure 4 captions have been updated to clarify that observation are from different locations within the site of interest. 

“Discussion section 4.1. The authors mention in L407 that soil emission factors can fall into several climatic and nutrient categories, but do not report which category the site in question falls into. It might be useful to mention that all published data taken for comparison with model predictions belong to the specific category. In this relatively long section, the authors look at various model results and compare them with various published data sets. The main message is somewhat blurred in this type of presentation. Perhaps a summary of all cited data in the table with a parallel presentation of the model prediction made in the current paper would be useful.”

A: We appreciate the reviewer’s comment and agree that the discussion could be clarified further. To address this, we have added a conclusive sentence at the end of the relevant paragraph in line 437:

Added text (line 437):

"Based on these conditions and considering the categories commonly used for emission factors in this land category, the simulated site would fall between the categories of nutrient-rich sites in the boreal zone and nutrient-rich sites in the temperate zone."

Regarding the suggestion to include a table summarizing all cited data alongside the model predictions, we opted not to include a table for the following reasons:

• The study reports different values for the same variable (e.g., GPP) at different points during the forest rotation, reflecting temporal variability that aligns better with the field-based data we are comparing against.
• For example, some GPP values are presented for specific periods when LAI (Leaf Area Index) was comparable to field observations, while others are presented for two years post-clear-cutting to match the context of specific reported field measurements.
• Creating a table would require including substantial supplementary information to explain the context and temporal variability behind each data point, which we believe is more effectively communicated through text.

“Discussion section 4.2, L481-490. The authors admit that linear decomposition rate constants are not adequate to reflect the real processes, but they also defend the approach used, saying that the model representativeness of the measurements is satisfactory.  Where is the truth, and what would the authors ultimately recommend? I assume that the use of multiple pools, as used in the current SOM decomposition concept, gives the model the necessary flexibility even for linear decomposition, but, unfortunately, there are no data to prove this, e.g. for cellulose, lignin and peat decomposition.”

A: We appreciate the reviewer’s insightful comment regarding the limitations of linear decomposition rate constants and the representativeness of the current approach. To address this, we have clarified and expanded the discussion to highlight the applicability and limitations of the model, as well as to provide recommendations for future improvements.

Modified text (from line 490):

“…where peat is defined as a conceptual compartment with unspecified chemistry, that decomposes linearly according to a rate constants modified by environmental conditions. Even though this description provides reasonable carbon dynamics at yearly and decadal time scale, limitations within this representation might explain why daily NEP fluxes were not well represented by the model."

Added conclusive statement:

"Generally, a model that explicitly represents the interactions between organic carbon substrates and microbial communities is desirable for exploring research gaps associated with priming effects and increased precipitation variability. The application of decomposition models based on linear kinetics to managed peatlands could benefit from representing SOM using measurable pools, especially if field-based decomposition data for different measurable SOM pools, coupled with field-based carbon balance data, become more readily available."

“L539: It is not entirely clear, what alternatives other than continuous forest cover were considered. Can your modelling approach be used to infer the applicability of alternative land uses, how should it be expanded? And the following statement… How can the full GHG budget (CH4 and N2O) can be considered? Is its consideration critical for the current study? If so, then the conclusions are somewhat vague...

A: We appreciate the reviewer’s comment and take this opportunity to clarify our approach and conclusions. Our study demonstrates that clear-cutting significantly affects the carbon metrics used, making the avoidance of clear-cutting a reasonable conclusion if the objective is to improve these metrics. While other management options exist, we highlight continuous forest cover (CCF) as a strategy worth exploring. To make this reasoning more explicit, we revised the paragraph in line 567 as follows:

Revised text (line 567):

"A very negative but improving ICS may be a representative pattern for these systems; therefore, avoiding clear-cutting is crucial to prevent deterioration of this metric. It is important to assess the effects of different management strategies, especially those that do not rely on clear-cutting, such as continuous forest cover (Laudon & Maher Hasselquist, 2023)

Regarding the full GHG budget, our intention is to discuss its limitations and applicability. For the current study, the ICS is useful as it provides a metric to assess the impact of the simulated site, with negative values being undesirable from a climate change mitigation perspective. This makes it suitable for evaluating a single trajectory of carbon exchange and comparing system boundaries. However, for comparative analyses between land-use scenarios, its applicability is limited. To clarify this point, we revised the paragraph as follows:

Revised text:

"While the ICS provides valuable information for evaluating the climatic impact of a specific trajectory of carbon exchange, as demonstrated in the present study, it does not account for the varying warming effects associated with different types of carbon compounds exchanged (e.g., methane emissions under waterlogged conditions). Therefore, to assess alternative land use scenarios for forested drained peatlands, such as rewetting, a metric that incorporates both temporal dynamics and the warming effects of all GHGs, such as cumulative radiative forcing (Murphy & Ravishankara, 2018), would be ideal."

“L 24: …consistent with the…”

A: Fixed

“L 30-31: A bit strange selection of three countries. What is about others? How big this coverage in respect to all drained peatlands in Northern Europe?”

A: The sentences has been modified to:

“Forestry on drained peatlands is a widespread land management practice in the northern hemisphere, covering approximately 15 million hectares, and it has significant implications on GHG budgets (Leifeld et al., 2019). This practice is particularly common in Fennoscandia, where it spans around 5.7 million hectares in Finland and 1.5 million hectares in Sweden (Vasander et al., 2003).”

“L 58: relevance (?)”

A: Fixed.

“L65: Mamkin et al is missing in reference list.”

A: Fixed.

“L75: I would suggest to formulate the first research question a bit differently, say more specific. For example, Do calibrated ForSAFE-Peat model capable to describe field-based observations.... ?”

A: Fixed, see next comment

“L77: Second research question also can be more specific. It is not clear what exactly you mean with system boundaries.”

A: To introduce what we mean by system boundary we modified the paragraph in the introduction as follows:

“The importance of the tree biomass components is clear from net ecosystem production (NEP) measurements performed with the eddy covariance technique, which indicate a persistent carbon sink in afforested drained peatlands despite high soil carbon losses (Korkiakoski et al., 2019; Meyer et al., 2013; Tong et al., 2024). It has been recognized that in cases of persistent and large soil carbon losses, compensation through forest uptake is limited because the tree component has a maximum carbon storage capacity lower than the carbon stocks of a typical peat soil. The magnitude and extent of this compensation are likely sensitive to how harvested wood products (HWP) are accounted for. When considering HWP, post harvesting periods are of special relevant suggesting that to understand the trade-off between tree biomass carbon and soil carbon, it is necessary to analyse carbon dynamics over more than one forest rotation. This shows how differences in system boundary definition, meaning considering the carbon balance within the soil, ecosystem, or the ecosystem plus the fate of HWP may lead to contrasting results.”

Furthermore, the questions were modified to:

1. Do output from calibrated ForSAFE-Peat model resemble field-based observations related to carbon dynamics in a northern drained afforested peatland?
2. Do output from calibrated ForSAFE-Peat model indicate different patterns of carbon exchange across different system boundaries for a northern afforested drained peatland?

“L124: Maybe add information about RCP 6.0 – e.g. moderate temperature increase”

A: The following sentence was added after the line mentioned by the reviewer:

“RCP 6.0 represents a medium stabilization pathway, where greenhouse gas emissions peak around 2080 and decline thereafter, reflecting a future with moderate climate change mitigation efforts”

“References: Sierra (2024). Please do not cite unpublished work. Check if it is accepted and delete if not at the moment of your publication.”

A: The reference has been eliminated and change for Muñoz et al 2024.

“L300: EDC – is not a common abbreviation, try to avoid it and at least explain at the moment of first appearance.”

A: The acronym is now defined before first appearance

“L302: It is not fully clear what you mean with the decreasing availability. Peat has always the same properties? Maybe explain this shortly here and not in the supplementary model description.”

A: We appreciate the reviewer’s comment and take this opportunity to clarify. Linear kinetics make the decomposition flux dependent on the amount of substrate available. In the second rotation, the first two soil layers experienced substantial decomposition, reducing the peat mass and, consequently, the amount of peat available for further decomposition. However, this reduction was compensated by higher decomposition rates, primarily driven by increasing soil temperature.

To address this point in the manuscript, we revised the sentence as follows:

Revised sentence:

 “The decreasing availability of peat in the first three soil layers over time, resulting from the reduction in peat mass due to decomposition, did not lead to lower decomposition fluxes because increasing soil temperature led to higher decomposition rates.

“L349: losses”

A: Fixed

“L483: It is well established…”

A: Fixed

 “Table 1: Full names for the NCB – net carbon balance and ICS- integral carbon storage should be given in the title. I would use the same dimension for ICS in all three system boundaries, namely 10^6. It helps to see immediately the difference.”

A: We agree, the table has been modified accordingly and now it the ICS for the system boundary Ëcosystem+HWP” reads -0.7×10⁶

“Figure 6: delete ‘plots’ in caption. I suggest to plot stacked areas for the combined figs a and b and for combined figs c and d. HWP in the plot description can be fully spelled.”

A: We agree is a good idea, the figure has been modified accordingly and is attached to this reply

“L285: Not clear what is said: CO2 emission is the most significant and then is stated that leached DOC contributed 10% and CO2 -4%, i.e. less than DOC.”

A: We meant leached CO2 dissolved in water, now is corrected in the sentence

“L511: what is DNM?”

A: Fixed, now in line 171 que introduce the acronym: Ditch network maintenance (DNM)

Answer to anonymous referee # 4:

We appreciate the reviewer’s insightful comments and the detailed attention given to the supplement. These observations have helped us identify areas for improvement that will enhance the clarity and quality of the manuscript. Below, we address the main comments and outline the proposed changes to the manuscript.

FIRST MAIN COMMENT:

“The authors use a pre-existing and tested model for simulating tree stand development and carbon dynamics, and incorporate a new description for organic soils, describing their carbon dynamics. The model is extensively described in the supplement. However, it was not immediately clear to me which parts of the described model are already published, and what is new here. This makes it hard to assess the work presented. I think the supplement would need a more clear distinction between these components.”

ANSWER:

The rationale behind the supplement was twofold: first, to provide a detailed explanation of how carbon dynamics are represented in the model, as this is the primary focus of the paper; and second, to highlight aspects of the model that have been improved for this contribution. However, we recognize that the distinction between previously published components and new developments was not made sufficiently clear. To address this, we will improve the first paragraph of the detailed model description in the supplement, explicitly clarifying which features are pre-existing and which are novel to this contribution. Additionally, we will incorporate clear statements during the mathematical descriptions to indicate where new features are introduced

“A general description of the previous version of the ForSAFE model can be found in Wallman et al. (2005). The model simulates plant dynamics based on the PnET-CN model (Aber et al., 1997), soil chemistry based on the SAFE model (Alveteg et al., 1998; Warfvinge et al., 1993), water dynamics based on the HBV/PULSE model (Andersson, 1988; Bergström, 1991) and soil decomposition dynamics based on the DECOMP model (Wallman et al., 2006; Walse et al., 1998).The model was further developed to include daily dynamics among other processes such as lateral water movement by Yu et al. (2018) and Zanchi et al. (2021).

In this contribution we refined carbon allocation and respiratory processes in the plant. Given different turnover and growing rates between foliage and roots, allocation to roots is constantly recalculated to follow foliage demand. Maintenance respiration is decoupled from photosynthesis and is now a function of plant tissue biomass and temperature. To better capture soil dynamic in peats, this version introduces a dynamic soil volume approach and considers anerobic decomposition of soil organic matter. Furthermore, soil nitrogen mineralization and mineral nitrogen transformations has been changed to recreate better soil microbial processes. Soil organic matter decomposition is coupled with a fixed microbial carbon to nitrogen ratio to calculate mineralization and immobilization, while mineral nitrogen is transformed through nitrification and denitrification with explicitly modelled microbial groups. Soil temperature dynamics have been improved by solving the heat equation tailored for peat soils.

While this description primarily focuses on carbon dynamics, other subcomponents are briefly explained, with an emphasis on changes made in this contribution. Carbon is represented as a set of compartments that denote different states of carbon within the system, such as carbon within foliage biomass, labile carbon allocated to roots, cellulose-like carbon in the soil, and so forth. A summarizing schematic of organic carbon pools within the model can be seen in figure S1.”

To further enhance clarity, we will explicitly highlight new developments during the mathematical descriptions. For instance:

Line 102 will be modified to: “For this contribution we change the description of carbon allocation to roots (equation 18). Unlike the other compartments, additional carbon for roots ( ) can be allocated at any time of the year if the carbon contents in root tissue ( ) and in the root labile compartment ) are not enough to fulfill the prescribed root to foliage ratio ( .”

Line 137 will be modified to: “Calculation of maintenance respiration was change for this contribution by decoupling respiration rates from photosynthesis rates. In this version we estimate total maintenance respiration  with tissue-specific ( ) linear rates scaled by tissue mass and temperature ( . Temperature response is the same for all tissues, following Aber et al. (1997). In the case of wood and branches, maintenance respiration occurs exclusively in the respiring wood fraction ( ).”

SECOND MAIN COMMENT:

“Also, more references would be needed: now many parameters were given, but with no source. Or at least clarify the source. For example, does Equation (59) in the supplement come from Manzoni & Porporato (2009), referenced in the previous sentence? Same applies to (at least, non-exhaustive list) lines 402 -- 404, 422 -- 423, 440 and the last three lines of Table S8.”

ANSWER:

To address this concern, we propose to be more explicit in indicating which equations are new to the model and their sources. This aligns with our response to the first comment, as both relate to improving clarity in the supplement.

For instance, we have revised the text around line 299 to specify the origin of Equation (59) and its adaptation for the model: “Nitrogen mineralization (denoted by 𝑁𝑀) from soil organic matter decomposition produces ammonium and is associated to the C:N ratio of the decomposition flux and the microbial biomass C:N ratio (Manzoni & Porporato, 2009). We have adapted the previously mentioned principle to both aerobic and anaerobic decomposition flux with equation 59.”

Regarding parameter sources, we recognize that linking parameters more directly to the equations in which they are used will enhance clarity. Below are proposed modifications as examples

Paragraph starting in line 397: “Variables used for plant carbon dynamics mentioned in the model description are grouped in Table S1 while parameters and their sources are grouped in Table S2. Most parameters come from articles presenting or applying the PnET model however we include new parameters associated to the new developments previously mentioned. Maintenance respiration rate constant (MR_TF, MR_TB, MR_TW and MR_TR) were obtained from Metzler et al. (2024), these parameters are necessary in equation 29. Also, a maximum root growth rate constant is introduced for Equation 18, with a value 0.05 d-1 derived from the upper range reported for Norway Spruce in a nutrient manipulation study (Sell et al., 2022). As mentioned in section 1.1, the parameter respiring wood fraction (RWF) modulates maintenance respiration of the woody tissue. The value was obtained by we manual calibrated against proxy variables observations related to biomass size (tree ring data for the period 2007-2009 and GWL data for the period 2008-2013).”

Paragraph starting in line 402: “Variables used for soil carbon dynamics are grouped in Table S5 while parameters and their sources are grouped in Table S6. Most parameters were derived from the original DECOMP model (Wallman et al., 2005; Walse et al., 1998) and recent modifications to ForSAFE (Yu et al., 2018). For this contribution the previous recalcitrant SOM compartment was renamed the peat compartment. Aerobic decomposition rate constant was taken from Clymo et al. (1998). This version of the model incorporates anaerobic decomposition based on linear kinetic as shown in equation 33. Therefore, we added to anaerobic decomposition rate constants ( ) for different soil organic matter compartments. For example, the anaerobic decomposition rate of the EDC compartment (PKan_SE) has been assumed to be 10% of its aerobic rate, while the anaerobic rate of the cellulose compartment (PKan_SC) is assumed to be 1%. based on ranges for anaerobic decomposition of polysaccharides in incubation studies (Benner et al., 1984). Limited lignin decomposition is assumed to occur anaerobically; the rate constant (PKan_SL) is assumed to be 0.01% informed by Reuter et al. (2024). For peat, the anaerobic decomposition rate constant (PKan_SL) is derived from catotelm decomposition rate constants (Clymo et al., 1998). Although anaerobic decomposition rate constants are less relevant under drained conditions, they are critical when the model is applied to waterlogged conditions.”

We have clarified the rationale for parameters related to plant-mediated gas transport (line 440): “Variables used for soil gaseous carbon dynamics are grouped in Table S7 while parameters and their sources are grouped in Table S8. Parameters associated with plant-mediated gas transport are not applicable to Norway spruce due to the absence of aerenchyma tissue in its root system. This explains the assigned values for parameters A_P, TRτ and RCS in table S8.

THIRD MAIN COMMENT:

“Related to that, the performance of the peat module may need more validation. Also, I found some pecularities about the peat fluxes: for example, seems that peat is formed from the other SOC pools. Is this realistic? ”

ANSWER:

To further explore and validate the model, we propose conducting a sensitivity analysis by altering nutrient content and water availability. This analysis will help us better understand the model's behavior under varying conditions. The results of this analysis will be included as an appendix. Additionally, we will provide a supplementary appendix showcasing the physical changes in the peat profile simulated by the model (e.g., bulk density and thickness changes). These simulated results will be compared with published data to offer an additional layer of validation for the peat module, which we believe will strengthen confidence in the model’s performance.

Regarding the question about peat fluxes, we argue that the current model description reflects a reasonable assumption grounded in field observations. Peat formation can be conceptualized as a humification process, where organic material is transformed through microbial activity and physical conditions. For instance, the Von Post scale, widely used to evaluate peat profiles, essentially measures the degree of humification. In our model, the first three soil organic matter (SOM) pools receive inputs exclusively from litter and represent a gradient of litter quality. The decomposition process, mediated by microbial activity and physical stress, results in the formation of an amorphous SOM compartment, which we term "peat." This approach aligns with other models, such as Yasso07 (Didion et al., 2014), where four SOM compartments associated with litter inputs decompose to form humus. In the YassoPeatland adaptation, this humus compartment is associated with the catotelm (Li et al., 2024). Similarly, the ESOM module in the SUSI peatland simulator follows this conceptual framework, where decomposition succession of litter-related compartments leads to peat formation (see reference: https://www.sciencedirect.com/science/article/pii/S0048969724053233).

Representing peat as a single conceptual compartment is a common practice in ecosystem models. For instance, the SUSI peatland simulator and models such as ORCHIDEE discretize peat into saturated and unsaturated layers. While we recognize the limitations of conceptual compartments, as noted in our discussion section, peat in real life is loosely defined. Chemically, it consists of varying proportions of microbial-derived compounds, cellulose-like compounds, lignin-like compounds, and others, with these proportions differing by peat type. Moving to a more analytically based compartmental structure might improve the representation of peat, but it would also reduce the model's practicality for many applications. Data on peat composition are far less common than data on peat thickness and carbon content, making a conceptual approach more viable for the moment.

To avoid confusion arising from the diverse chemical composition of peat, we propose renaming the SOM compartment unrelated to litter as "humus," as used in models such as Yasso07. This change would align with established conventions and clarify that, at the beginning of simulations, all organic material is assumed to be humus, similar to the assumption in He et al. (2016). This adjustment, combined with the sensitivity analysis and validation appendix, should address the reviewer’s concerns while maintaining the model's usability.

FOURTH MAIN COMMENT:

“Line 80 in main text: ForSAFE has a mineral soil weathering component. This plays a key role in nutrient availability and forest productivity, but the dynamics are radically different on organic soils. How is this accounted for?”

ANSWER:

The weathering component in ForSAFE is based on the initial conditions of the soil's mineral fraction, which is further subdivided into specific mineral types, as described in the SAFE model. This component is particularly relevant for nutrient availability as it provides a source of base cations. However, the differences in weathering dynamics between mineral and organic soils emerge naturally in the model due to the site-specific initial conditions.

In the peat-dominated site considered in this study, the mineral soil fraction is very small, comprising around 13% of the total soil composition, SOM is around 87% of the soil. While this fraction is incorporated into the model, its limited presence means that weathering has a negligible effect on nutrient availability due to a very small flux. Instead, nutrient dynamics in this site are predominantly governed by organic matter mineralization processes, which are more representative of nutrient cycling in organic soils. This distinction ensures that the model accurately reflects the key drivers of nutrient availability and forest productivity under the specific conditions of the study site.

By explicitly representing both the mineral and organic soil components and their respective contributions to nutrient availability, the model captures the unique dynamics of organic soils while maintaining the flexibility to simulate sites with different soil compositions.

MINORS COMMENTS:

Abstract: first line of abstract mentions greenhouse gases. For organic soils, also CH4 and N2O may play an important role: maybe emphasise in the abstract that you only focused on carbon dioxide?

A: We do include CH4, the site is just a small CH4 sink. To clarify this in line 16-18 has been modified as follows: “Simulated carbon fluxes were analysed and compared under different system boundaries (soil, ecosystem, and ecosystem plus the fate of harvested wood products named ecosystem+HWP) using the net carbon balance (NCB) and the integrated carbon storage (ICS) metrics. Carbon loses were indicated with negative values while carbon gains were indicated with positive values.”

Abstract: maybe give brief explanation on the metrics used, especially ICS? Also, clarify whether the negative or positive values indicate carbon loss or gain, conventions on this vary between metrics.

A: We will explore possibilities to add an explanation of ICS in the abstract.

Line 37: C in gC often not in subscript

A: I agree that is not often in subscript but that has never sat well in my head. When information is added to a unit (in this paper and in many others) often subscripts are used. For example, when is cubic meter of water you often find m³_water or m³_w I have not found things as m³W so in order to be consistent I rather go with g_C.

Lines 63 -- 64: add ref

A: We added Lehtonen et al. (2023) (https://www.nature.com/articles/s41598-023-42315-7)

Line 84: maybe intensively monitored instead of heavily?

A: Indeed, sounds better. We have changed it to “intensively”.

Line 103: is this a new development?

A: Not really, lateral flow was added by Zanchi et al. (2021). We have clarified this in a new paragraph, see answer to first main comment.

Line 106: mineral weathering: how is this described if the peat layer is thick?

A: We have tried to answer to this concern in the answer to the fourth main comment.

Lines 113 -- 114: add reference for decay constants

A: Here we just try to describe the model structure, we think that the supplementary material is detailed enough regarding the sources and values associated to the decay constant and other elements of the model

Line 117: is there a better reference for the station?

A: We changed to Klemedtsson et al. (2015) (https://ui.adsabs.harvard.edu/abs/2015EGUGA..17.7461K/abstract)

Line 118: high nutrient content: nitrogen, or also others? Due to agricultural history, was it fertilised earlier?

A: We don’t have information for every nutrient but CN ratios and NP ratios show high contents of nitrogen and phosphorus. Likely the agricultural history played an important role, but the deeper peat (1.5m to 2m) show high content of nitrogen and phosphorus suggesting nutrient content being associated with peat type.

Line 122: give also distances to the stations, coordinates hard to interpret

A: We have modified the paragraph to: “The model used daily mean meteorological data (1951 to 2023) from the Swedish Meteorological and Hydrological Institute (SMHI) Uddevalla (58°36’N, 11°93’E) and Vänersborg (58°35’N, 12°35’E) stations, both located approximately 12 km from the site. Future climate data (2023 to 2088) were obtained from projections for forest sites under the CLEO research program (Munthe et al., 2016). Climate projections were downscaled from regional projections based on ECHAM and HADLEY climate model under RCP 6.0 as in (Zanchi, et al., 2021a).”

Line 125: deposition of nitrogen?

A: Yes, among other elements for which deposition data is also available such as Cl and Na.

Line 130, Fig. 1 caption: is the difference between minima and maxima plotted, or minima and maxima explicitly? If not explicitly, why not?

A: We agree is unclear, we have a changed the caption to: “Figure 1. (a) Mean annual temperature (black line) and the range between the annual maximum and minimum temperatures (grey area)”

Line 140: is the assumption of only Norway spruce realistic?

A: We believe the assumption is reasonable, there was not any other species of tree in the site and the cycle of forest stand dominates carbon fluxes. However, we bring the limitations associated to not representing understory vegetation in the section (4.1) of discussion associated to representativeness.

Line 140 -- 141: little snippets of information like this make it hard to follow what is old, what is new. Maybe restructure to have all model related stuff in one place

A: This section aims to provide information of the main characteristics of the scenario simulated in the model and relate those characteristics to information we have about the site and similar sites. For this case we reformulated to avoid discussing detailed parametrization, we reformulate the paragraph to:

“The forest stand at the site was strongly dominated by Norway spruce (Picea abies), and the simulation used parameterizations specifically developed for this species in previous studies (Aber et al., 1996, 1997; Zanchi et al., 2021b). The modeled forest management replicated historical events at the site: spruce planting in 1951, a 72% tree biomass thinning in 1979, a 10% biomass loss in 2010 due to storm damage, and a 96% biomass removal in 2019 as part of a clear-cutting operation. Harvesting plays a crucial role in regulating carbon dynamics in such systems. The large thinning event, which removed 72% of the biomass approximately 28 years after planting, represents a non-conventional management practice (Metzler et al., 2024). This intensive management strategy was incorporated into the simulations to accurately reflect the historical management of the real site. The second modeled rotation (2020–2088) followed the same biomass removal timing patterns as the first rotation.”

Line 142: is three metres the depth of the peat?

A: Is an average because the site was a fen valley, therefore depth changes spatially being deeper at the center.  

Line 144: same properties for all layers: is this realistic? I would imagine different peat decomposition stages at different depths. Although I understand that this may be hard to verify

A: Indeed, giving the uncertainty about the conditions in 1950’s and spatial (horizontal and vertical) variability exhibited by the site we consider that giving all layers the same bulk density, CN ratio and OM content is a reasonable simplification. We know we can’t perfectly represent the site (digital twin) so the intention is to create a synthetic site that represents the main characteristics.

Line 146: ref for the OM content

A: We added a reference, now it reads: “while organic matter content was set to 87% based on Meyer et al. (2013)” 

Line 179: Is NCB same as net ecosystem exchange or productivity? What fluxes are included, also lateral?

A: We explain this in line 193.

Line 201, Fig. 2 caption: arrow a 1: this refers to leaching to runoff? Currently the arrow points upwards, and looks like offgassing from ditches

A: I agree but we try to clarify that in the text, we decided like this to simplify what is inflow (going down) and what is outflow (going up). In some cases, like leaching or harvested wood (6 in figure 2b), is not so logical, however we think it helps to quickly define inflows and outflows regardless of the true spatial characteristic of the flow (harvesting and leaching being lateral losses). 

Line 202, Fig. 2 caption: arrow a 3: does the soil itself take up atmospheric carbon dioxide? Litterfall is separately described. Also, SI says there are no processes that consume gaseous CO2 (line 243)

A: Carbon gas exchange is gradient controlled, so in theory could flow both ways. In these particular conditions (drainage) the site emits CO2 and consumes CH4.

Line 202, arrow a4: would this not be outflow? Or does this mean release of CO2 into the soil headspace/matrix? I would guess that is rapidly lost through arrow 2

A: We are concerned about flows of carbon so CO2 is not the only carbon flux. In this case this is carbon in litterfall being added to the soil

Fig. 2 caption: some arrows stay the same between a, b and c, such as arrows 1 and 2. Some change (like arrows 3, 4 and 5). I would suggest keeping the meaning of each number the same, and having a few more numbers.

A: I think you make good points to improve Fig 2. We will follow your advice regarding adding more numbers.

Lines 587 -- 589, code availability. Would it be possible to make the model openly available? It is possible to get a persistent doi for a model release e.g. through Zenodo github integration. Do the contact details here refer to DE? Are there any guarantees that these work in ten years, and the model is still available? Even if the details would be permanent, I would prefer open availability of the model. I understand that this view might be not shared by everyone, but I would urge the authors to consider this possibility to promote open science. 

A: We agree with the reviewer’s concern, the reason of the current statement on code availability is associated to current work being performed in the code to include certain aspect we need to include for a future contribution of the effect of restoration. Therefore, we would like to have a final version of the code with the new improvements before releasing it. 

Line 592, data availability. Great that the data is already available at this stage! However, please address the comment from another reviewer regarding the data coverage. 

A: We will indeed fix the file associated to the data used for this contribution.

Comment on the figures: please consider increasing the font size in the figures.

A: We will do that

Fig. S1: The definition of the abbreviations is not clearly given in the text, it took me some looking to find them in Table S1 much lower down. 

A: We have changed the figure caption to provide the missing information: “Figure S1. Scheme of carbon compartments within the model. The green boxes represent plant-related compartments, the brown boxes represent soil organic matter compartments, the grey boxes represent deadwood left after harvesting and harvested wood products (HWP), and the pink boxes represent CO2 and CH4 in the soil. The black arrows indicate carbon fluxes, with arrows not connected to any compartment representing fluxes leaving the system.  The subscripts denote specific carbon compartments: LC (labile central), LF (labile foliage), LR (labile root), LB (labile branch), LW (labile wood), TF (tissue foliage), TR (tissue root), TB (tissue branch), TW (tissue wood), HW (hardwood from harvest), HP (paper from harvest), HF (fuel from harvesting), TD (tissue deadwood), SE (soil easily decomposable compounds), SC (soil cellulose), SL (soil lignin), SP (soil peat) and SD (soil dissolved).”

Multiple points in manuscript and SI, for example line 157 in SI: linear decay of harvested wood. I assume first-order exponential decay is meant?

A: Thank you for noticing it, we will correct accordingly.

SI, Line 295: modified how exactly? Nitrogen dynamics on organic soil are very different from those on mineral soil

A: See answer to the first main comment. We improve the first paragraphs of section 1 of the supplementary material to provide more clarity to the changes performed to the model for this contribution.

SI, section 2.2: I initially was confused by the Harvest removal and Harvest intensity parameters, maybe these could be clarified in the text?

A: We have added this sentence to the section 2.2: “The Harvest Intensity parameter ( ) specifies the proportion of the forest stand that is cut during harvesting, whereas the Harvest Removal parameter ( ) defines the fraction of the harvested material that is removed from the site.”