Reply of Comment on gmd-2021-280

Referee comment on "Assessing methane emissions for northern peatlands in ORCHIDEEPEAT revision 7020" by Elodie Salmon et al., Geosci. Model Dev. Discuss.,

https://doi.org/10.5194/gmd-2021-280-RC1, 2021

We thank the reviewer for his/her constructive comments. Our manuscript has been significantly improved through addressing the reviewer’s comments.

Below, comments have been numbered and a response is exposed for each of them. Line numbers correspond to the Preprint version. Modifications are highlighted in bold in the present document.

Anonymous Referee #1

Comment 1.1: This paper presents a land surface model with an explicit representation of northern peatlands, ORCHIDEE-PCH4. The model simulations are compared to data from 14 wetlands. This paper focuses on methane emissions and refers to previously published work for peat carbon accumulation and carbon balance. The authors use the root mean square difference between simulated and observed methane emissions to optimize 7 model parameters. They first perform the optimization separately for each of the 14 sites, then perform a multi-site optimization.

General comments:

The paper is well written. The introduction, the model description and the site description are very clear. The optimization method is very hard to understand (although I am not a specialist). The results are sometimes hard to follow with very small figures. The discussion and conclusion are clear.

Response: We thank reviewers for positive comments and are glad that our work and conclusions are well understood. By addressing reviewer comments and suggestions below we hope that the optimized method section will be easier to understand.

Comment 1.2: My main comment is related to the representation of some processes in the model. From what is shown by the authors it seems like for most sites, methane emissions are pretty much independent of the water table depth. It is particularly obvious for the sites with

multiple years of data. Fr-Lag for instance has a simulated high water table the first summer followed by two summers and autumns with very low water tables. Methane emissions are low the first summer and increase the two following ones (contrary to the observed fluxes). This behavior can also be seen at DE-Sfn, Fi-Lom, Pl-Kpt, and to a lesser extent at DE-Hmm and Dk-Nuf. This contradicts most of the existing literature on observations showing a strong correlation between water table depth and methane emissions (the higher the water table, the higher the emissions).

Response: As we pointed out at the beginning of the discussion section L524-527 “Sensitivity analyses were previously performed to assess methane emission models responsiveness to parameters values (Meng et al., 2012; Riley et al., 2011; Spahni et al., 2011a; Wania et al., 2009; Zhu et al., 2014). These studies (Van Huissteden et al., 2009; Riley et al., 2011) suggested that temperature dependency of methanogenesis is the most influential parameter affecting methane production whereas methane emissions are mostly sensitive to oxidation and plant transport.”  Indeed, by definition (equation 3, L175-177) “the rate of methanogenesis (k_i in s^-1) depends on soil temperature and moisture according to the same function as for the heterotrophic respiration (Qiu et al., 2019)”. And our results displayed in the manuscrit, Figure 2 to 5 and in supplementary document, Figure S3, show that methane emissions are correlated with optimum when both simulated soil temperature and moisture conditions are the highest.

Besides, in our model, as explained L493-495 “the simulated water table position is a prognostic variable defined by the cumulative amount of soil water content over the soil column (Fig. S2 and Fig. S3).”  We also demonstrate for two sites, at US-Los and DE-Spw, in the supplementary document Figure S4 and S5 and discussed in the manuscript L495-496 that above the simulated water table position, soil moisture is still higher than 80% which is sufficient for methanogenesis to occur in the model. This explains that the correlation between the simulated water table position and simulated methane emissions is not as strong as the one observed in the field.

Comment 1.3:  In terms of processes, for a same site a higher water table is related to a higher soil moisture content and a lower oxygen concentration over the whole peat column, which favors methanogens (and also limits methanotrophs). It is particularly important to correctly represent the link between peat water content and methane emissions for a model designed to be used in climate change studies. I would request the authors to at least clearly discuss this issue in the paper and modify their conclusions accordingly.

Response: As explained in the response to comment 1.2, we demonstrate in the manuscript (L175-177, L493-495, L524-527) and with additional figures in supplementary materials (Fig. S3) that our methane model is correlated with simulated soil temperature and moisture conditions. Indeed, Fig S3 shows a strong correlation at each site of methanogenesis maximum with both soil temperature and moisture maximum. We also discussed in the manuscript (L493-495) with additional figures in the supplementary (Fig. S2, S4 and S5) that our prognostic water table position defined from simulated soil moisture content is not well correlated to the observation water table positions measured on sites.  We added a few sentences in the conclusion to highlight these results L707: “Our results show that as in previous methane emissions models (Meng et al., 2012; Riley et al., 2011; Spahni et al., 2011a; Wania et al., 2009; Zhu et al., 2014), simulated methanogenesis is strongly correlated to simulated soil temperature and moisture content whereas methane emissions are more strongly correlated to plant mediated fluxes and to soils methane oxidation proportion. We have to point out that in our model, a weak correlation has been established between the observed water table positions and the prognostic water table positions established from simulated soil moisture content. A correlation between soil moisture content and water table position in the field is needed to improve representation of the water table position in models.”

Comment 1.4: The second comment is related to the Dk-Nuf site. I happen to be familiar with this dataset and I noticed some imprecisions in the text (see specific comments), so I would ask each co-author responsible for a site to carefully proofread the manuscript.

Response: We addressed the reviewer's concerns by addressing specific comments below (comment 1.5).

Specific comments:

Comment 1.5: L 243: Dk-Nuf : The methane emissions are measured by automatic chambers on this site. There is a flux tower but it only measured CO2 fluxes (besides turbulent energy and radiative fluxes). Also, the water table depth is not available at this site. Table 1: measurements for Dk-Nuf don’t cover 2006-2009 but 2008 – 2014 (actually, the dataset extends to 2019)

Response: We thank the reviewer for noticing this mistake, we carefully checked the observation period of table 1 and modified the one of FI-Lom to 2006-2009 and DK-Nuf to 2008-2013 that were wrong.

Table 1. Sites ecological characteristics summary. Sites identification includes the country initials and the short three letters name of each site, locations of the sites are provided by the country, latitude (Lat) and longitude (Lon) values. Hydrological characteristics are distinguished by the type of ecosystem, fen, bog, tundra and marsh. Y and N indicate presence and absence of snow cover in winter, permafrost soil, forest above the peat. Temporary drawdown of the water table level is specified by presence and absence indicators Y or N. Grey color highlight groups of peatlands organized by amount of methane emissions in ranges 0-10, 10-150, 150-400, 400-600 mg m^-2 d^-1.

Comment 1.6: L 266 : I don’t understand the 0.5 degree grid cell. What is actually run at this resolution?

The authors say at line 270 that they impose site level meteorological forcings. They also seem to indicate that the spin-up to reach close to observed peat carbon content and depth was done by using the site specific meteorological data. So is it the texture that is at 0.5 degree or the hydrology? If it is the hydrology to calculate the peatland fraction then what is used to force this calculation? A gridded meteorological forcing or the site specific one ? If a gridded meteorological forcing was used then it should be mentionned. This is a bit confusing.

Response: to improve the description of simulation setup, we modified paragraphs L266: “Each peatland site is a sub-grid area embedded in the 0.5°x 0.5° grid cells whose extent is determined by a fraction of grid area as defined in Table 2. These sub-grid areas enable the representation of ecosystems variability in which a specific scheme simulates soil hydrology, vegetation characteristics and soil carbon cycling for northern peatlands. The fraction of peatlands per grid cell was defined by modifying the prescribed values employed by Qiu et al., (2018) in order to collect enough water to fill the peatland by runoff from the other soil fractions and elevate the water table level for northern peatlands. We employed vegetation phenotype properties and peatland fractions described in (Qiu et al., 2019) and peatlands hydrology and carbon model as described in Qiu et al., (2019). Site simulations were then constrained at the grid cell scale with a half hourly time series of meteorological conditions e.g. air temperature, wind speed, wind direction, longwave incoming radiation, shortwave incoming radiation, specific humidity[ES1] , atmospheric pressure, and precipitation. These time series are flux tower measurements that were gap filled by 6-hourly CRU-NCEP 0.5◦ global climate forcing dataset (Qiu et al., 2018). ”

Comment 1.7: Table 2: I am surprised by the value of observed carbon stock at DK-NuF. The only study known to me (Morel et al, Earth Syst. Sci. Data, 2020) gives 36.3 kg/m2. This is much lower than the 54.6 given in the Table.

Response: We added in table 2 citations in which we collected the maximum peat depth and the soil carbon stock. For DK-Nuf we used and refer to measurements from Bradley-Cook and Virginia, (2016)

Table 2. Simulations conditions and framework to constrain peatlands soil carbon stock. Grey color reports the groups of sites with equivalent levels of methane emissions (Table 1).

Comment 1.8: Section 2.3: this is very hard to follow. I am absolutely not a specialist but I wondered if the whole time series of observation was used, and at what time step? (hourly, daily,monthly, yearly?) and why.

Response:  Optimization simulations were performed over site-specific observation periods as defined in Table 1. We added these precisions in the sentences L305 “Two types of simulations are performed over the site-specific observation period defined in Table 1: single site (SS) experiment for which parameters are optimized for each site and a multi-site (MS) experiment that aims at refining one set of parameters considering all sites together.”

The timestep of the model is half-hourly but the model output is the mean monthly value. We choose to have the same timescale for the outputs throughout the sites following measurements timescale that we all converted to monthly timescale because it was the timescale of the chamber measurements. We added measurement timescale L284-288: “These sites are a subset of the 30 peatlands sites collected for the calibration of ORCHIDEE-PEAT (Qiu et al., 2018) for which, in addition of eddy-covariance data and physical variables (water table, snow depth, soil temperature), methane emissions were measured by eddy-covariance at daily time scale at US-Los, hourly timescale at DK-Nuf and otherwise at half-hourly timescale or chamber measurements at monthly timescale for FR-Lag and RU-Che. All methane emissions data were monthly average.”

Comment 1.9: L326-327: this is very strange. If zroot is increased to 0.75m then, if I am not mistaken, froot=0 (I am assuming zsoil=0.75m since this is the peat depth), so that fpmt=0. If zsoil is not equal to 0.75m, increasing zroot decreases fpmt (in absolute value). Is this wanted?

Similarly, why increase the rate of methanotrophy to get higher methane net emissions?

Response: Indeed, as in Walter and Heimann (2000), froot is decreasing with depth meaning that proportion of root decreases with depth. Nevertheless, since the plant mediated transport scheme is limited to layers containing roots by increasing zroot we increase the number of layers that will be involved in methane transport. Concerning the rate of methanotrophy, when the rate is higher, a lower proportion of methane is oxidized which leads to higher content of methane in the soil that will eventually be emitted. We modified the sentence L324: “Three parameter ranges were modified for DK-Nuf, the minimum value of q_MG was lowered to 7.0, z_root maximum is increased to the maximum peat depth at 0.75m in order to consider plant mediated transport in all the peat layers, the maximum value of T_veg was increased to 40.0 and the maximum rate of methanotrophy k_MT was enlarged up to 8 d^-1to decrease the methane oxidation and to obtain in the simulation methane emissions higher than 150 mg CH₄ m^-2 d^-1. ”

Comment 1.10: Another question: what did the authors do with missing data for methane emissions? There are very few winter measurements of methane emissions at Dk-Nuf. Did the authors gap-fill the data and then optimized their parameters on this? It should be clearly stated.

Response: To clearly explain data availability and timescale, we added fews sentences L243 : “For the optimization simulations, at DE-Sfn, DE-Hmm, FI-Lom, PL-Kpt, PL-Wet, and US-Wpt, year-round data were available and zero values were filled for the first and the last month of years at the beginning and the end of the observation period. Otherwise, winter months were filled with zero and during spring, summer and fall months missing data were filled gapped using a linear regression.”

Comment 1.11: L327: to be coherent with the rest I believe Table 4 should have values in parenthesis for qMG at PL-Wet site (4 is outside the range given in table 3)

Response: PL-Wet out site range has been added under parenthesis to table 4:

Table 3. List of parameters driving the methane production, oxidation and transport scheme in ORCHIDEE-PCH4.

Comment 1.12: Figure 2-5 a): why not give the observed water table depths when available

Response: The observed water table depths when available are already provided in the supplementary document Figure S2 and overlapped with simulated water table depth. More details about correlation between observed and simulated water table depth are available in the response to comment 1.2 and 1.3.

Comment 1.13: L 459: how does permafrost explain a deeper simulated water table position ? explain.

Response: The simulated water table position is defined as the accumulation of water content height in each soil layer. When a soil layer, or part of it, is frozen water infiltration is reduced involving reduction of soil moisture in the deepest soil layers and a lower water content height and consequently a lower water table position. We modified the sentence L458: “Indeed, both sites are underlaid with permafrost which limit water infiltration to the deepest soil layers and can explain these deeper simulated water table positions.”

Comment 1.14: L 562: I couldn’t agree more with the authors comment on the need for more data on vascular plants in peatland

Response: This sentence is: “While a significant number of studies provide insight on gas exchanges through vascular plants, densities of vascular plants with aerenchyma in peatlands is poorly characterized.”

Comment 1.15: L 590-591: I am not sure the authors really showed that these 2 sites were limited in methane substrate. It is likely the case in the model, but is it the case in reality ? Because it seems that this model result might be related to the partitioning between active, slow

and passive C pools.

Response: We agree with the reviewer that for these 2 sites methane substrate is limited in our model. We modified the sentence to specify this L590: “Only two values have been defined above 10 at US-Wpt and DK-Nuf which are two sites that are limited in methane substrates in the model which explains these high values of T_veg. ”

Comment 1.16: L 727: I couldn’t agree more with this last sentence.

Response: This sentence is: “This demonstrates the complexity of interactions of the methane cycle with environmental conditions considered at various scales and the need for more detailed on-site studies.”

Technical comments:

Comment 1.17: L103: Qiu et al, missing year

Response: We added citation year L103 “Recent developments of ORCHIDEE land surface model lead to simulate soil hydrology, permafrost thermodynamic and carbon cycle in the northern latitudes (Guimberteau et al., 2018) and in the northern peatland specifically (Qiu et al., 2018), including peat carbon decomposition controlled by soil water content and temperature as well as CO₂ production and consumption processes (Largeron et al., 2018; Qiu et al., 2018)”

Comment 1.18: L 285 : there is something strange with the end of the sentence “… and driven”

Response: Words ‘ and driven’ have been removed, the sentence read now L285: “Then during the site-specific measurement periods (Table 1), methane variables are calibrated against observed monthly average methane fluxes times series.”

Comment 1.19: L316 : reached instead of reach

Response: The sentence read now :  “Successive runs serve to ensure that the minimum reached is not a local minimum.”

Comment 1.20: L318 : that emits instead of that emitting

Response: L318 ‘One of these four sites, DE-Spw, is among the sites that emits the fewest amount of methane (up to 7 mg m^-2 d^-1) and features a larger stock of carbon of 84 kg C / m² than at US-Los that features 27 kgC / m² and emits up to 4 mg m^-2 d^-1. “

Comment 1.21: L357 : “discharged” is not the right verb here

Response: The verb has been removed, sentence L 357 read now “At sites that emitted between 10 and 150 mgCH₄ m^-2 d^-1, RMSD values fluctuate between 4 and 26 and when methane fluxes were between 150 and 400 mgCH₄ m^-2 d^-1, RMSD is of 38 - 80.”

Comment 1.22: L365 : “significantly lower” : that is quite an understatement

Response: L365: “The temporal and the average magnitude are equivalent than in measurements except for the US-Wpt and RU-Che for which simulated emissions are much lower than observed emissions.”

Comment 1.23: L 385 : “simulated diffusion of atmospheric methane” instead of “diffusion of simulated atmospheric methane”

Response: L385:”This explains the negative methane flux (Figure 2c) produced in winter by the model via simulated diffusion of atmospheric methane in the snow cover. “

Comment 1.24: Figure 2, 3, 4 and 5: b and c : the units of the flux should be mg/m2/d on the Y axis e : on the Y axis on the right

Response: Units of Figures 2 to 5 line b,c and e have been modified to mg/m2/d for the Y axis.

Comment 1.25: L 482: Table 6 instead of 5

Response: We modified the table number, the sentence now reads L482: “Multi-site optimized parameters values acquired by using average values of parameters defined at each site and the initial ranges (Table3) are shown in Table 6. ”

Reply of Comment on gmd-2021-280

Referee comment on "Assessing methane emissions for northern peatlands in ORCHIDEEPEAT revision 7020" by Elodie Salmon et al., Geosci. Model Dev. Discuss.,

https://doi.org/10.5194/gmd-2021-280-RC1, 2021

We thank the reviewer for his/her constructive comments. Our manuscript has been significantly improved through addressing the reviewer’s comments.

Below, comments have been numbered and a response is exposed for each of them. Line numbers correspond to the Preprint version. Modifications are highlighted in bold in the present document.

Anonymous Referee #2

Comment 2.1: This is my first time reviewing this manuscript. The authors evaluated single-point applications of one land model ORCHIDEE in simulating methane emissions at 14 northern peatland sites. The PEAT model version they employed has been incorporated with several previous methane algorithms.

General comments:

The paper is well written and organized. The topic on peatland methane simulation is also very interesting because it remains a challenge to perform the methane processes well for most land surface models. The ORCHIDEE-PEAT model can be a potential helpful addition for the wetland modeling studies. However, in the current manuscript there are still some problems regarding the methods and results the authors need to clarify.

Response: We carefully revised the manuscript by addressing each of the comments and suggestions.

Comment 2.2: 1. The study sites are northern peatlands and the authors have mentioned the importance of snow layer on calculating the diffusion of ch4 flux in the study, but I never found any data or plots of snow observations/simulations to support their claims, such as L385-386, L403-407, L432-433, L466 in the manuscript.

Response: The methane model has been optimized employing observed methane emissions. Our goal has never been to evaluate or optimize every process (methanogenesis, methanotrophy and gas transport) since there is not enough observation data available in the literature to achieve this. Nevertheless, it is of significance to highlight the contribution of each process in order to fully explain simulated results of methane emissions. More specifically these “claims” described results from Fig 2 to 5, line e, that display simulated methane content in the model snow layers. The snow model has been fully evaluated in Guimberteau et al., (2018) which has been reported L118.

Comment 2.3: 2. For wetland simulations, the anaerobic environment is crucial for methane production/oxidation. It is mainly controlled by water conditions and microtopography. I found the bad performance of methane emissions always with poor water table simulations in the study. for example the US-Wpt and RU-Che sites in fig5, Fig2s, L455-466. How did you consider the impacts or limitations from the biases of hydrological simulations?

Response: ORCHIDEE is a global scale land surface model (LSM), while it is commonly employed to perform site simulations for calibration purposes or study specific processes, it is mainly dedicated for global scale studies. At global scale, all LSM methane models are optimized against observed methane emissions which implicitly enabled compensation effects between processes which are at the origin of methane emissions. In order to assess model uncertainties that will occur during global scale simulations, we choose to consider simulation conditions as if the model was employed for global scale simulations and to enable compensation effects between processes involved in methane emissions that include biases of the hydrology model. For instance, water table positions at US-Wpt in our simulations are fluctuating between 0 and -40m which is also the case at FI-Lom, DK-Nuf and DE-Sfn.  While simulated methane emissions are in good agreement with observations for FI-Lom, DK-Nuf and DE-Sfn, methane emissions are largely underestimating at US-Wpt. Therefore, we have not found any obvious relationship between a lower water table position and an underestimation of methane emissions. In addition, simulated water table positions are prognostic variables estimated from the simulated soil moisture content and the correlation between observed and simulated water table position is provided in supplementary Figure S2. We also demonstrate in supplementary document Figure S3 that at each site maximum methane emissions are correlated with the optimum of both soil temperature and moisture. Finally, we show at US-Los and at DE-Spw, in supplementary document Figure S4 and S5 that above the simulated water table position soil moisture content is higher than 0.8 which is sufficient for methanogenesis to occur. Our results show that conditions or processes promoting the largest methane fluxes such as US-Wpt and RU-Che and the lowest one such as at US-Los and DE-Spw are the least understood.  A better understanding of these processes/conditions will serve to better constrain models.

Comment 2.4: 3. Some key information on parameter definitions, parameter values and units were missed in the method section. It is a bit hard for readers to understand the methane modeling structure. See specific comments below.

Response: In order to improve the method section and ease its reading we have addressed the reviewer's specific comments. Point by point responses are described below (see comments 2.10 to 2.25).

Comment 2.5: 4. The authors should include more details on the parametrization method. (1) Please provide the details of how you determined the final accumulation year for ending the simulations at L279. Because there are two factors of soil peat depth and carbon content, which factor is prior? What are the referential resources of soil peat depth and carbon stock?

Response: We added some details on how we defined soil carbon profile in the model, L278 “The peat model (Qiu et al., 2019) enables a vertical buildup of peat by simulating a downward movement of C when the discretized organic layers reach a threshold defined from a regression relationship between the carbon fraction and measured bulk density. This scheme in ORCHIDEE-PCH4 serves to constrain the vertical distribution of the soil carbon stock to the observed maximum peat depth. Simulations with ORCHIDEE-PCH4 driven by repeated site-specific meteorological conditions were performed for various periods of time to reach the observed soil carbon content and maximum peat depth (Table 2). “

We also added references in Table 2 and full references details have been added to the references list at the end of the manuscript:

Table 2. Simulations conditions and framework to constrain peatlands soil carbon stock. Grey color reports the groups of sites with equivalent levels of methane emissions (Table 1).

Comment 2.6: (2) For you parameter sensitivity analysis at L305-310, which period did you run for the analysis?

Response: We are not sure to understand this comment, optimization experiments serve to define “the set of parameters that reduce the gap between observation and simulation data.”(L304). Therefore, for each site we run the simulation over the entire observation period defined in Table 1 for which we have temporal monitoring of methane emissions. Moreover, this is not a sensitivity analysis per se, this is more a comparison between two optimization strategies. The first one is based on single sites and we assumed it will perform better for a site but each site will at the end have their own parameter sets. Since ORCHIDEE is a global model, we needed a single set of parameters that we obtained using the multi sites parameterization. Thus, the single site optimization is the best the model can do whereas the multi-site is the best trade off we can find.

Comment 2.7: 5. I have also a couple of concerns on the parameter settings of single-site and multi-site simulations. (1) There are several conflict descriptions on parameter ranges between single-site and multi-site simulations the authors should clarify. See specific comments below.

Response: We  addressed the reviewer's concerns by addressing specific comments below (comments 2.28 to 2.30 and 2.40 and 2.41).

Comment 2.8: (2) When doing parameter optimization for the multi-site simulation, why didn't you set all parameter ranges to cover the already obtained single-site optimized value, such as L482, L491-499 qmg, zroot and Tveg. For example, the specific optimized value of qmg at PL-Wet is 4; according to the range of qmg from 9 to 10 in table 6, the final multi-site optimized value will miss its single-site optimized value. This may be one reason why the simulations with multi-site optimized values for PL-Wet and DK-Nuf were worse than their single-site simulations.

Response: We performed simulation experiences that the reviewer proposed using a so-called extended range that include all the parameters values obtained after the single site optimization. Results are shown in the supplementary document Table S3 to S5 and Figure S9 and S10. These results are also reported in the manuscript L500-503.

Comment 2.9: 6. How did you consider the limitation for peat depth simulations in section 4.2? In table 2 I found there are several sites whose peat depths do not consist of observations, such as DE-Zrk, DK-NuF, which should impact the distribution of soil carbon within 0-0.75m.

Response: The peat depth limitation is defined based on the diagrams of maximum methane production depth display in Figure 2. To explain this, we added that sentence L625 “In Figures 2e to 5e, that display the depth of maximum methane production, reveal that the deepest methane production depth is 0.75 meters in all the simulation results.“

We also checked the carbon content value integrated up to 1.5m which corresponds to the depth of the model layer below the 0.75 meter layer. Even though the numbers change for some sites, the conclusions drawn are the same.

Specific comments:

Comment 2.10: In equation 1-13 of section 2.1, do variables, such as fmga, fmgs, fmagp,.... fmt, also change with layer depth and time? I assume z is depth and t is time in the equation 1, despite you not giving their definitions.

Response : To ensure that the reader understand that these variable change in depth and time we modify the sentence L150: “where each term that varies in time (t) and with depth (z), expresses methane production […] “

Comment 2.11: Also lacks the related information on parameters. For example, are kmt of equation 5 at all layers the same value?

Response : Yes, k_mt and all the other parameters are the same value for all model layers unless it is mentioned otherwise.

Comment 2.12: L164-165 What's the units of [C], ki? 

Response: The unit for [C] has been added L164 :“each type of soil carbon pools ([C]_i, i = a, s, p; in gC m³ of soil), active, slow and passive “ and ki  L176: “where the rate of methanogenesis (k_iin s^-1) “

Comment 2.13: L165 For soil, it should be 'soil moisture', not 'soil humidity'.

Response : We changed it to soil moisture:“where the rate of methanogenesis (k_i) depends on soil temperature and moisture “

Terms ‘soil humidity’ were changed to ‘soil moisture’ in all the document, L128, L166, in Table3 (‘Connectivity of soil moisture‘), L381, L535, L550 and L606.

Comment 2.14: L169 Why didn't you mark O2, p and O2* using brackets as [CH4] and [C]? For example, is the O2 here different from [O2]soil in Fig1?

Response : O2, p and O2* are notations employed by Kvonostianov et al. (2008) but the review is right it is more consistent to harmonize the notation within the manuscript. 

We modified equation 3 to

And L169 to : “ [O₂]_p is the oxygen concentration in the soil per unit porous volume, [O₂]_anoxia is the soil oxygen”

Comment 2.15: L174 What's the unit of clay content? For these peatland sites, what are the values of clay you used?

Response : The clay content is a fraction whose values are defined between 0 and 1. Here it is defined at 0.2. We modified L174 “where clay is the clay fraction and has a value of 0.2, “

Comment 2.16: L177 Again, what's the difference among [O2]soil, O2 and O2*? And their units should be added.

Response : To improve the definition of O2,p and O2* We modified L169-170 “ [O₂]_p is the oxygen concentration in the soil [O₂]_soil(in gO₂ m^-3 of soil) per unit porous volume ( , 𝝅_soil is the soil porosity), [O₂]_anoxia is the soil oxygen concentration at which anoxic conditions are reached and enable methane production. This oxygen concentration threshold is assumed to be 2 g m⁻³(Duval and Goodwin, 2000).”

And [O2]soil is defined L175-176 : ”The amount of methane consumed by methanotrophy is limited by the soil oxygen concentration, [O₂]_soil, following a 1:2 CH4:O2 molar ratio”

Comment 2.17: L178 What is the exact unit for kmt, hour or day?

Response : We modify the unit L178 “where k_MT is the rate of methanotrophy which value range from 0.06 to 5 d^-1 (Morel et al., 2019) “ to be consistent with k_MT unit in the tables.

Comment 2.18: L184 Are you sure that kebu= '1 h-1'? This would mean the constant is unuseful during calculation.

Response : Here the kinetic variable is defined for the physico-chemical sense of the equation. However, since no observation value are available, it is not possible to constrain the threshold for ebullition by a kinetic rate therefore in the present study and in Khvorostyanov et al., (2008a, 2008b), the kinetic rate is defined to 1h-1.

Comment 2.19: L185 What is the unit of Psoil? Why did you set the 0.75m depth as limitation?

Response : The unit of Psoil is Pascal (Pa) we added to the manuscript L185 “pressure (Psoil in Pa) “

The limitation at 0.75m is defined L191-194 “We estimated that in our model below the layer corresponding to 0.75 m the hydrostatic pressure is always higher than the partial pressure of dissolved gases. Therefore, we considered below 0.75 m that methane ebullition threshold is constant and equal to the value defined at 0.75 m in order to avoid methane accumulation in the deeper layers. “

Comment 2.20: L192-193 Again, how did you estimate the layer of 0.75m? Please add more information or results.

Response : The full explanation is provided L 190 “It has been suggested that ebullition in soil occurs when the partial pressure of dissolved gases exceeds the hydrostatic pressure (Chanton and Whiting, 1995). We estimated that in our model below the layer corresponding to 0.75 m the hydrostatic pressure is always higher than the partial pressure of dissolved gases. Therefore, we considered below 0.75 m that methane ebullition threshold is constant and equal to the value defined at 0.75 m in order to avoid methane accumulation in the deeper layers. ”

This description is in the method section, we think it will be confusing to provide results in that section. For more details readers should refer to the work of Chanton and Whiting, (1995) that is cited.

Comment 2.21: L199 What's the unit of tortuosity n? Does it have a constant value?  

Response : L200 have been modified to “the tortuosity n that depicts the sinuous path of bubbles, is defined to be 2/3 (Hillel, 1982). “

Comment 2.22: L201-202 Do you mean wsize = 1cm here?

Response : The sentence, L201-202, has been modified to “Khvorostyanov et al., (2008a, 2008b) defined wsize = 1cm for a carbon rich loess deposit of the Yedoma. “

Comment 2.23 : L206-207 It is not clear. Do you mean that the oxidized methane in the root zone accounts for 39-98% of total methane oxidation in soil? Please provide the exact value used in your study.

Response : We modified L206 to “The proportion of methane oxidized (M_rox) in the root zone is emitted as CO₂ to the atmosphere. Walter and Heimann, (2000) estimated M_rox to range between 39 and 98% of methane located in the root zone. “ 

We explained in section 2.1, L139-142 that the parameters values are defined in section 2.2. In order to put forward the presentation of the method subsections we moved sentences L139-142 ahead after the section title, 2.Model description, and added a sentence to introduce subsection 2.1 ORCHIDEE-PCH4. This part reads now L111:”A general presentation of ORCHIDEE-PCH4 and of processes considered in the methane model are exposed in section 2.1. Implementation of methane production and oxidation and transport are specified respectively in sections 2.1.1 and 2.1.2 whereas parameters values established for the formatting site simulations conditions before observation periods are given in section 2.2. Then, section 2.3 describes the parameter optimization approaches. “

Comment 2.24: L210 Again, what's the value of Tveg used in your study?

Response : As explained in response to comment 2.23 above, parameters values are defined in section 2.2. This is specified in line 111-115.

Comment 2.25: L215 I didn't find the root biomass in the root distribution function 10, only the root depth. Please clarify this.

Response : We defined in a different way the f_root  L215 to “This function describes the vertical distribution of roots in the soil in which z_root is the rooting depth and z_soil the soil depth. “

Comment 2.26: L219 Is the z in this equation the same as that in other equations? How did you distinguish the depths of soil, water or snow coverage?

Response : In the model, soil carbon and the hydrology schemes are vertically discretized into 32 layers corresponding to a depth of 38 meters. Each layer is ascribed a soil depth and a layer thickness. The snow is divided in 3 layers on top of the soil layers. We added L217 “The gas diffusion scheme features the diffusion of CH₄ and O₂ in the 3 top layers of snow when snow cover is formed and in the 32 soil layers that correspond to 38 m depth. This scheme considered (1) the diffusion of oxygen from the top soil to the soil layer, (2) the diffusion of methane produced and remaining in the soil and (3) methane exchange between the soil and the atmosphere at z=0:”

Comment 2.27: L224 Please provide the sources for the values.

Response : We added the reference L224 (Khvorostyanov et al., 2008a)

Khvorostyanov, D. V., Krinner, G., Ciais, P., Heimann, M. and Zimov, S. A.: Vulnerability of permafrost carbon to global warming. Part I: Model description and role of heat generated by organic matter decomposition, Tellus, Ser. B Chem. Phys. Meteorol., 60 B(2), 250–264, doi:10.1111/j.1600-0889.2007.00333.x, 2008a.

Comment 2.28: L286-287 Did you run the carbon accumulation process again for the historical period of each site with these calibrated parameters?

Response : We added this sentence to be more accurate, L287 “Then site-specific simulation over the observed period is run again using the optimized parameters.”

Comment 2.29: L305-310 For the parameter sensitivity analysis. Which period did you run? Please clarify this.

Response :  We added this precision to the sentence L305: “Two types of simulations are performed over the site-specific observed period defined in table 1: single site (SS) experiment for which parameters are optimized for each site and a multi-site (MS) experiment that aims at refining one set of parameters considering all sites together.”

Comment 2.30: L324-327 Can you explain why the optimized qmg at PL-Wet is so low at 4? The value is beyond much from the observed ranges in other studies. Do you have any evidence to support this value? The questions are also about qmg, zroot and Tveg for DK-Nuf.

Response : Since, methane results from the decomposition of soil carbon in order to produce large amount of methane, large amount of soil carbon need to be degraded in the model. Therefore, L321 we started by explaining that “The other three sites for which some of the optimized parameters are out of the initial range, DK-Nuf, PL-Wet and US-Wpt, are among the sites that emit more than 150 mg CH₄ m^-2 d^-1. The carbon stock at DK-Nuf and PL-Wet are respectively 55 and 38 kg C / m² which is lower than at FI-Lom and PL-Kpt that accumulated more than 200 kg C / m².” This suggests that for both PL-Wet and DK-Nuf methanogenesis is limited or near to be limited by the soil carbon content in the model. Then we explain L327 “PL-Wet required also to modify range values of q_MG to 1.0-11.0 leading to the lowest optimized q_MG value of 4.0 which significantly reduced the RMSD from 227.4 to 80.5 (Fig. S1 and Table S1).” Indeed, a smaller minimum of the cost function have been found when extending the range values of the parameters leading to q_MG = 4.0. In the discussion section we point out that L532 : “This ratio [q_MG] was first established from experimental studies that determine the microbial production ratio CO₂ to CH₄ (Potter et al., 1996; Segers, 1998) for various water table positions. These ratio values were found to be between 0.58 and 10000." Then we provided model values employed in the literature for q_MG L 534 : “Because of this wide range of values, process - based models employed this CO₂ to CH₄ ratio as an adjustable parameter that is weighted by environmental factors such as soil humidity and temperature. Wania et al., (2009) performed a sensitivity analysis study of the LPJ-WHyMe model using 7 sites in which the multi-site optimization value of the CO₂/CH₄ ratio was defined at 10 while other models as CLM4Me use a value of 5.” Therefore, we disagree with the reviewer, a q_MG value of 4 for PL-Wet or 7 for DK-Nuf is not beyond much from the observed ranges nor from values employed in other models. This value only reflects, as we pointed out L321, explain here above and discuss in section 4.2, the limitation of methanogenesis by soil carbon content in the model.

For DK-Nuf, we added some details in the results description L324 “Three parameter ranges were modified for DK-Nuf, the minimum value of q_MG was lowered to 7.0, z_root maximum is increased to the maximum peat depth at 0.75m in order to consider plant mediated transport in all the peat layers, the maximum value of T_veg was increased to 40.0 and the maximum rate of methanotrophy k_MT was enlarged up to 8 d^-1to decrease the methane oxidation and to obtain in the simulation methane emissions higher than 150 mg CH₄ m^-2 d^-1.” Here again, these parameters values result from the minimization of the cost function for a site where methanogenesis is limited in soil carbon content in the model. In the discussion section 4.1 we explain L572 that in the present study “Optimized z_root values at sites ranges between 6 and 68 cm depth with the average depth defined at 26 cm which is also the value obtained using the multi-sites approach. These values are consistent with values utilized by Walter and Heimann (2001) that ranged between 0 and 74 cm.” For T_veg, we explain that it is arbitrary defined L586 “T_veg has been introduced by Walter et al., (1996) to describe the density of plants and their efficiency in methane transport for site estimation. It is an adjustable parameter that was scaled to be between 0 and 15 with lower values for ecosystems dominated by trees and shrubs and the highest values for ecosystems dominated by grasses and sedges.” The range of T_veg has been established by Walter and Heimann (2001) based on the study of 4 peatlands sites and a swamp site which may not be representative of all northern peatland sites.

Comment 2.31: L335 Why the range of qmg here is not 4.0-10.7 as before?

Response : We modified L335 to “Across sites, q_MG values extend between 4.0 and 10.7, optimized k_MT values vary between 1 and 5.25 d^-1. ”

Comment 2.32: L381/L550 'soil humidity' should be ''soil moisture'

Response : Terms ‘soil humidity’ were changed to ‘soil moisture’ in all the document L128, L166, L381, L535, L550 and L606 and Table3.

Comment 2.33: L385-386 Please provide the snow cover depth data to support the claim on diffusion transport.

Response : These sentences are referring to results in Figure2 line c and d which show line e, the simulated amount of methane fluxes via diffusion, ebullition and plant transport and line d, on the right hand side of the y axis, the simulated amount of methane contained in the simulated snow layers. We added references to Figure 2 in the sentence L385 “This explains the negative methane flux (Figure 2c) produced in winter by the model via simulated diffusion of atmospheric methane in the snow cover (Figure 2d).”

Comment 2.34: L403-407, L432-433, L466 Again, please provide the snow depth evidence.

Response : Similarly than for comment 2.33, we added references to figures in sentences L402 “As for US-Los and DE-Spw, at CA-Wp1, during the winter simulations show that in the top soil layers some methane is transferred by diffusion (Figure 3c) to the snow cover (Figure 3d). Then a small part of the simulated methane is temporarily stored in the snow (Figure 3d) and the other part is released to the atmosphere via diffusion (Figure 3c). More simulated snow accumulated at DE-Sfn, DE-Zrk, CA-Wp1 and US-Bog where up to 0.8 - 0.04 gCH₄ / m² are temporarily stored in the snow (Figure 3d). At FR-Lag and DE-Hmm, fewer methane, less than 0.005 gCH₄ / m², are contained in the simulated snow cover (Figure 3d).”

L432 ”In the winter the methane fluxes are stored in the simulated snow cover at FI-Lom (Figure 4d), therefore the simulated surface fluxes above the snow are driven by diffusion (Figure 4c). However, during summer simulated methane fluxes essentially originate from plant mediated transport. At DK-Nuf, PL-Kpt and PL-Wet, simulation results show that fewer methane, less than 0.4 gCH₄ m^-2 d^-1, accumulates in the simulated snow during winter (Figure 4d). ”

L464 “Though, at RU-Che, simulated methane production rate is higher around 100 mg CH₄ m^-2 d^-1 and occurs at 20 cm depth during summer and few centimeters below the surface during winter.”

Comment 2.35: L411, Can you explain why there are so large differences in the depth of maximum methane production among sites. The main controlling factors for methanogenesis are temperature and water conditions. Did you check the relationships/dynamics of soil temp along depth?

Response : The reviewer is correct, variability of the depth of maximum methane production results from the effect of soil temperature and moisture. This relationship derived directly from the implementation of peat decomposition function described in detail in Qiu et al. (2019). In the supplementary document, Figure S3 displays the temperatures and soil moisture at the depth of maximum methane production 

Comment 2.36: L440 "not sufficient to cause methane ebullition" How much ch4 concentration is enough for the ebullition under your study conditions? Please clarify this.

Response : The ebullition threshold is defined by equation 7 (L186). It varies in space and time and depends on the soil pressure and temperature and on the mixing ratio mxr_CH4. In addition to this threshold, the remaining amount above the methane threshold is affected by a probability function defined in equation 8 that depends on soil moisture and wsize. Therefore, it is difficult to define how much CH4 concentration will result into a flux by ebullition since it is a complex multivariable problem. However, Figure 4c shows the simulated amount of methane that is computed by the ebullition, the diffusion and the plant transport scheme. There are no methane fluxes by ebullition for FI-Lom, there are only methane fluxes by diffusion (in blue) and plant transport (in green).

We added references to the Figure 4 in the sentence L438 : “In contrast, at FI-Lom simulated soil methane concentrations are near 50 gCH₄ m^-2 during summer whereas the winter concentrations are near 80 gCH₄ m^-2 (Figure 4d) which is not sufficient to cause methane ebullition (Figure 4c). Indeed, the ebullition (equation 7 and 8) results from the balance of soil temperature, pressure and gas content, which explain the large diversity of methane fluxes response by ebullition at each site. ”

Comment 2.37: L534 Can you give some citations or explanations for the qmg range of 0.58~10000 because in table 3 you state the range is 9-11.

Response : We do not understand the reviewer comment 2.37 since the citation for the qmg range of 0.58~10000 are provided L528 (Potter et al., 1996; Segers, 1998), the value used by other models are reported L533 and explanation for the qmg range in table 3 is provided L537 “Khvorostyanov et al., (2008a) and Morel et al., (2019) used respectively q_MG value of 9 and 10 to simulate methane emissions from arctic peatlands therefore in the present study at first q_MG were optimized in the range of 9-11 then this range was enlarged only for sites that underestimate methane emissions. “

Comment 2.38: L554 Still lacks snow depth or accumulation evidence.

Response : Here we are discussing the model implementation and results to enforce these two aspects we modified sentences L551 to “Thus, the optimization of the oxidation rate results from the balance between model inputs and outputs that are respectively available methane and oxygen substrates and methane fluxes which explain this large variability in oxidation rate. In addition, in our model, snow is considered in the diffusion scheme which is in part controlling diffusivity of oxygen from the atmosphere to the ground in winter (e.g. Figure 2c).”

Comment 2.39: L641 What are the full names of 'SS' and 'MS'?

Response : To ease the reading of the manuscript we added in the section title L 314 “3.1 Single site optimization (SSO)”, L475: “3.2 Multi-site optimization (MSO)” and L641: “by employing single site (SS) and multi-site (MS) optimized parameters”

Comment 2.40: L645-646 I am curious that did the plant transport dominate total ch4 flux for the whole year, or it only occured in the summer?

Response : Plant transport process is defined by equation 9 in which “The leaf area index (LAI) influences the methane flux depending on the growing stage of the plants” (L216). This means that methane plant transport occurs during the plant's growth period. We modify the sentence L645: “Plant mediated transport (PMT) were the largest simulated fluxes, during the plant’s growth period.”

Comment 2.41: Table 3, For kMT,  it is not consistent with the values at L178; for Mrox,  it is not consistent with the values at L206-207; for mxrch4, the citation is not the same at L188;

Response : We revised the citation of Morel et al. 2018 to Morel et al. 2019 in table 3. The citations of Walter and Heimann, 2000; Riley et al., 2011 have been added for mxr_CH4. The range for M_rox defined by Walter and Heimann, (2000) is based on the study of 4 peatlands sites and a swamp site which may not be representative of all northern peatlands sites therefore we employed a larger range 0-1 rather than 0.39-0.98.

Comment 2.42: Table 6, Why is the span of mxrch4 not the same as table 3?

Response : We do not understand this comment, the range for mxrCH4 are the same in both table 3 and 6. Table 3 presents the parameters and the default values employed for the simulation set up as explained in section 2.2 and Table 6 presents the multi-site optimization results discussed in section 3.2.

Comment 2.43: Fig 5 with Fig2s, The performances of water conditions for US-Wpt and RU-Che look not so good. How do you consider this?

Response : This comment is the same as comment 2.3, please see our response to comment 2.3.

Technical corrections:

Comment 2.44: L149 Full names of fmga, fmgs and fmgp?

Response : L150 “methane production (fMG, MG: methanogenesis, a: active pool, s: slow pool, p: passive pool) “

Comment 2.45: L186-187 I found both mxrch4 and mwrch4, which one is correct?

Response : L187 “where mxrCH4 is the methane mixing ratio in the bubbles”

Comment 2.46: L189 What's the full name of 'RR'?

Response : RR is the ideal gas constant and is defined at L189. To avoid any confusion we modified it to its more common name  R in equation 7 and L189: “It is converted to gCH4 per unit porous volume by the ideal gas constant (R), “

Comment 2.47: L189 "MwCH4 and BCH4 the Bunsen... " should be "MwCH4 and BCH4 are the... "

Response : To ease the reading, we modified this sentence L188-189 to “It is converted to gCH4 per unit porous volume by the ideal gas constant (R), MwCH4 and the Bunsen methane solubility coefficients (BCH4). “

Comment 2.48: L293 Again, mwrch4 or mxrch4?

Response : L293 “transport (mxrCH4, wsize, Tveg, zroot) “

Comment 2.49: L685 Citations?

Response : L685 “global methane budget for natural wetlands located northern of 30°N of 12-61 Tg CH4 y-1 for bottom-up approach and 31-64 Tg CH4 y-1 for top-down approach (Saunois et al. 2020) “

Comment 2.50:

L185 'depend' to 'depends'

L482, Not correct.  'table 5' must be table 6.

Fig 1, 'Lai' should be 'LAI'

Response : The manuscript has been modified respectively at L185, L482 and in Fig 1, following specific review suggestions.