Modelling boreal forest's mineral soil and peat C dynamics with the Yasso07 model coupled with the Ricker moisture modifier

As soil microbial respiration is the major component of land CO₂ emissions, differences in the functional dependence of respiration on soil moisture among Earth system models (ESMs) contributes significantly to the uncertainties in their projections.

Using soil organic C (SOC) stocks and CO₂ data from a boreal forest–mire ecotone in Finland and Bayesian data assimilation, we revised the original precipitation-based monotonic saturation dependency of the Yasso07 soil carbon model using the non-monotonic Ricker function based on soil volumetric water content. We fit the revised functional dependency of moisture to the observed microbial respiration and SOC stocks and compared its performance against the original Yasso07 model and the version used in the JSBACH land surface model with a reduction constant for decomposition rates in wetlands.

The Yasso07 soil C model coupled with the calibrated unimodal Ricker moisture function with an optimum in well-drained soils accurately reconstructed observed SOC stocks and soil CO₂ emissions and clearly outperformed previous model versions on paludified organo-mineral soils in forested peatlands and water-saturated organic soils in mires. The best estimate of the posterior moisture response of decomposition used both measurements of SOC stocks and CO₂ data from the full range of moisture conditions (from dry and xeric to wet and water-saturated soils). We observed unbiased residuals of SOC and CO₂ data modelled with the moisture optimum in well-drained soils, suggesting that this modified function accounts more precisely for the long-term SOC change dependency according to ecosystem properties as well as the contribution of short-term CO₂ responses including extreme events.

The optimum moisture for decomposition in boreal forests was found in well-drained soils instead of the mid-range between dry and water-saturated conditions as is commonly assumed among soil C and ESMs. Although the unimodal moisture modifier with an optimum in well-drained soils implicitly incorporates robust biogeochemical mechanisms of SOC accumulation and CO₂ emissions, it needs further evaluation with large-scale data to determine if its use in land surface models will decrease the uncertainty in projections.

1 Introduction

Soil moisture and soil C stocks in boreal forests are higher in forested peatlands on frequently paludified organo-mineral soils and in peatlands on water-saturated organic soils than in well-drained forests on mineral soils (Weishampel et al., 2009; Ťupek et al., 2008; Bhatti et al., 2006; Hartshorn et al., 2003). Almost a quarter of the total terrestrial C (440 Pg C) stored in boreal moist and dry soils has accumulated since the last glaciation (Scharlemann et al., 2014) and is expected to create large C losses under warming climates (Hararuk et al., 2015). Moist organic soils are crucial for modelling dynamics of the global C cycle as they store 5 times more carbon than dry mineral soils (Leifeld and Menichetti, 2018; Turetsky et al., 2015; Scharlemann et al., 2014). However, soil organic carbon (SOC) stocks modelled by Earth system models (ESMs) show large uncertainty due to structural model differences (Hashimoto et al., 2017; Hararuk et al., 2014, 2015; Todd-Brown et al., 2013) and differences in environmental drivers and their functional dependencies used by soil C models (Thum et al., 2020; Ťupek et al., 2019; Falloon et al., 2011).

Despite soil moisture being a dominant driver of variation in C dynamics (Humphrey et al., 2021), ESMs lack consensus on the response of decomposition to soil moisture and temperature (Yan et al., 2018; Sierra et al., 2015; Falloon et al., 2011). The functional forms of the temperature and moisture modifiers of default decomposition rates among models disagree in their representation of extreme cold–dry and hot–wet conditions (Sierra et al., 2015). For example, the moisture decomposition dependency in the Yasso07 soil C model (Tuomi et al., 2011, 2009) is based on annual precipitation, has a form of monotonic saturation, and is uninformed about soil characteristics. By a monotonic saturation function, we mean a function which is entirely non-decreasing, initially increasing rapidly and later slowly approaching a maximum. The use of the saturation function is limited to well-drained soils because under wet or poorly drained forest soils such a model results in underestimation of the C stocks (Dalsgaard et al., 2016, Ťupek et al., 2016). The soil module of the CENTURY model (Adair et al., 2008; Parton, 1996; Metherell et al., 1993) uses precipitation and basic soil data (bulk density, clay, and silt contents) to calculate soil moisture, which similarly to Yasso07 assumes saturation of decomposition rates. Other functional dependencies of moisture such as DayCent, DEMETER, STANDCARB, CANDY, Gompertz, Myers, Moyano, and Skopp assume various functional forms, for example, non-monotonic Gaussian increase with an optimum and reduction in decomposition (DayCent, STANDCARB, Moyano), linear increase until reaching an optimum and linear reduction (Skoop), monotonic linear and non-linear increase (DEMETER and Myers, respectively), or monotonic saturation functions (CENTURY, CANDY) (Sierra et al., 2015, and references therein). The wide variation in commonly used moisture functions may result from the variety of data from different soil types and climates used to constrain these moisture indices.

If environmental response functions were calibrated for mineral soils only, then these functions may not adequately represent responses in the moisture range characteristic of organic soils. For example, default response functions of soil C models cannot represent anoxic inhibition of decomposition rates in paludified peatland forest soils. However, the inhibition of decomposition can be accounted for even in monotonic functions, e.g. by a reduction parameter such as “anerb” in CENTURY (Methereall et al., 1993). Due to variable water levels and their determination of soil oxic/anoxic conditions and SOC accumulation in peatlands, peat SOC stocks are typically modelled with models specifically developed for peatlands (Bona et al., 2020; Kleinen et al., 2012; St-Hilaire et al., 2010; Frolking et al., 2010, 2001; Clymo, 1978, 1992). However, for global applications on peatlands, the general soil models in ESMs can be modified for peat soil by adjusting parameters such as the hydraulic conductivity, as seen in models like JULES (Chadburn et al., 2022) and ORCHIDEE (Qiu et al., 2018), or by reducing decomposition rates for wetlands, as in LPJ (Wania et al., 2010) and JSBACH (Goll et al., 2015). The land surface model JSBACH coupled with the Yasso soil C model adopts a heuristic 65 % reduction in decomposition for wetlands (Kleinen et al., 2021; Goll et al., 2015). Using the CENTURY model at the site level, Raich et al. (2000) opted for improvement in modelled SOC of wetlands by modifying the environmental function with the anoxic inhibition for sites with insufficient drainage. This approach improved CENTURY compared to the default Yasso07 in poorly drained forested peatlands in Sweden, though the SOC stocks of both models were still underestimated (on average by 10 and 13 kg C m⁻², respectively) (Ťupek et al., 2016). A similar magnitude of SOC underestimation of Yasso07 model with default dependency on precipitation was also observed for poorly drained forest soils (e.g. Gleysols and organic soils) in Norway (Dalsgaard et al., 2016).

Figure 1 Infrared areal image showing the location of nine studied forest–mire types forming a transect of approximately 450 m on a northern hillslope in Finland (61°47^′ N, 24°19^′ E). A series of hemispherical images of forest stands on the top of the aerial image show the increasing gradient in the canopy openness from upland forests (left) to mires (right). The series of soil profiles shows the increasing gradient of the organic layer depth. The images in the series arranged from left to right mimic the site type location on the slope from the hill to depression. Sites range from upland (1) xeric, (2) sub-xeric, (3) mesic, and (4) herb-rich forest types (CT – Calluna, VT – vitis-idaea, MT – myrtillus, OMT – Oxalis–myrtillus) through paludified forest–mire transitions (5–7) (OMT+ – Oxalis–myrtillus paludified, KgK – myrtillus spruce forest paludified, KR – spruce pine swamp) to sparsely forested mires and peatlands in depressions (8–9) (VSR1 and VSR2 – tall-sedge pine fen).

We hypothesized that the SOC stocks and CO₂ emissions of mineral and organic (peat) soils can be modelled accurately by revising the original precipitation-based environmental modifier of a parsimonious model like Yasso07 with a function accounting for the reduction in decomposition based on the long-term near-surface moisture. Near-surface moisture is strongly correlated with the groundwater level depth in peatlands (Dimitrov et al., 2022), and the moisture values between mineral soil forests and peatlands are comparable on the same scale, which makes soil volumetric water content (SWC) a suitable variable for representing landscape moisture variation. Boreal forest SWC can be either measured in situ or derived in high resolution using hydrological models (e.g. Leppä et al., 2020; Launiainen et al., 2019) and at larger scale by remote sensing and machine learning (Han et al., 2023). We aimed to develop the original Yasso07 model with parameters as in Tuomi et al. (2011) by adding a revised unimodal moisture-based environmental function. We then optimized this function using Bayesian data assimilation of measurements from a boreal forest–mire hillslope catena of mineral, organo-mineral, and organic soils, and we tested whether we could correctly reconstruct observed SOC stocks and CO₂ emissions.

2 Methods

2.1 Study sites

Nine forest and mire site types in this study were situated along the hillslope from the Vatiharju esker to the Lakkasuo mire in southern Finland (61°47^′ N, 24°19^′ E) (Fig. 1) and formed a forest–mire ecotone, a gradient in soil moisture and nutrient status, vegetation composition, biomass production, and SOC stocks (Dimitrov et al., 2014). The sites were situated along a 450 m transect on a 3.3 % slope facing NE with a relative relief of 15 m. The site typology described below was based on the vegetation composition reflecting site wetness, fertility, and location on the slope according to Finnish forest and mire classification systems (Cajander, 1949; Laine et al., 2004).

On the crest of the esker was a well-drained xeric Scots pine forest (CT – Calluna type), which changed down the slope to sub-xeric mixed Scots pine–Norway spruce forest (VT – vitis-idaea type) and mid-slope to mesic and herb-rich Norway-spruce-dominated forest (MT – myrtillus type, OMT – Oxalis–myrtillus type), together referred to as mineral soil upland forests. On the toe of the slope were forest–mire transitions on gleyic organo-mineral soils or mixed spruce pine birch forests (OMT+ – Oxalis–myrtillus paludified type, KgK – myrtillus spruce forest paludified type, KR – spruce pine swamp type). On the level were water-saturated sparsely forested mires on Histosols (VSR1 and VSR2 – tall-sedge pine fen types).

The understorey or forest floor vegetation along the ecotone changed from being dominated by Calluna and Vaccinium vitis-idaea dwarf shrubs and typical forest mosses on the uppermost sites (CT, VT) to Vaccinium myrtillus dwarf shrubs with herbs on the mid-slope (MT, OMT); Vaccinium myrtillus dwarf shrubs with herbs and Sphagnum in the transitions (OMT+, KgK, KR); and Vaccinium oxycoccos and Betula nana dwarf shrubs with Menyanthes trifoliata, Carex, and Sphagnum species on the level (VSR1, VSR2) (Fig. 1). More detailed tree stand, soil, and climate characteristics for these sites were reported by Ťupek et al. (2015).

2.2 Auxiliary measurements

Soil temperature, water content, and CO₂ emissions ($\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) were measured simultaneously during the years 2004, 2005, and 2006. The measurement campaigns were conducted weekly, and we measured each plot once and all plots in 1 or 2 d between 07:00 and 18:00 eastern European time (EET) during the vegetative season of 2004 (July–November), 2005 (May–November), and 2006 (May–September) and monthly during the non-vegetative seasons (December–April). The summer seasons of the years 2004–2006 showed exceptionally different monthly weather patterns. Data from Finnish meteorological station, located in Juupajoki 3 km north-east of the ecotone, showed that the summer season in 2004 was rainy and colder in comparison to long-term typically mild weather, in 2005 the weather was typical, and in summer 2006 the weather was sunny and warm. The exceptional drought in 2006 caused by the lack of rain and increased temperatures in June and late July–early August (Gao et al., 2017) caused visible drying of the moss layer along all the sites of the ecotone. The 2006 summer drought ended with showers in mid-August and with more frequent rain in autumn when the soil moisture recovered to a normal level.

2.2.1 Soil temperature and moisture conditions

The soil temperature was measured at depths of 5 cm (T₅, °C) with a portable thermometer, and the soil volumetric water content was measured at a depth of 10 cm (SWC₁₀, %, m³ m⁻³) in all sites with a portable ThetaProbe (Delta-T Devices Ltd.) calibrated for each site type. The SWC calibration accounted for the bulk density–porosity of forest-type-specific soils (Ťupek et al., 2008, 2015). Because the forest–mire variation in the soil organic layer bulk density was relatively small at 0.34 ± 0.07 g cm⁻³ (porosity 74 ± 5 %) (Ťupek et al., 2015), the SWC values of the top 10 cm were on the same order of magnitude between the forest–mire site types. For missing field campaigns during months with snow cover (November 2004, February–April 2005, December 2005–April 2006), we interpolated the measured monthly mean T₅ and SWC₁₀ time series with a spline function. The SWC₁₀ values among the forest and mire site types ranged between 0 and 1 (or 0 % and 100 %) (Fig. 3), whereas in comparison to water level depth, the values range from 8 cm in tall-sedge mire to 881 cm in pine forest on the top of the esker (Ťupek et al., 2008).

2.2.2 Soil CO₂ emissions

Measurements of forest soil heterotrophic respiration (R_h, $\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$, positive sign) were taken using opaque cylindrical chambers (30 cm diameter, 21.2 L) placed on metallic collars (30 cm in diameter), which were installed permanently at 30 cm soil depth. The collar locations (3 for each site type, 12 for mineral soil forests, 9 for transitions, and 6 for mires, totalling 27) were selected to represent the spatial variation of each site type and the spatial variation along the forest–mire ecotone (e.g. dominant forest floor vegetation, microtopography, soil drainage, and nutrient status).

The aboveground forest floor vegetation inside each collar was clipped at the time of collar installation, and any plant regrowth of, for example, mosses was clipped approximately half an hour before the flux measurements. At the time of the collar installation, the roots of the understorey vegetation and trees were cut with a saw along the collars' perimeter. The metallic collars installed to 30 cm soil depth prevented the regrowth of the roots. Due to the vast majority of tree and understorey roots in boreal forest occurring in the humus layer, the 30 cm depth was considered sufficient to cut the roots and thus remove the signal of the root autotrophic respiration from the net CO₂ emissions. In transitions and mires, the depth of peat could be more than 30 cm (in a range from 0.15 m in OMT+ to 1.2 m in VSR2; Ťupek et al., 2008), but the prevailing high groundwater levels (in a range from 33 cm in OMT+ to 7 cm in VSR1, Ťupek et al., 2008) limit the root growth into the upper/sub-surface layer.

The soil CO₂ emissions were measured every 4.8 s during 80 s intervals with a portable infrared CO₂ analyser (EGM4, SRC-1, PP Systems Inc.). We calculated CO₂ flux rates from the development of CO₂ concentration over time inside the chamber.

2.2.3 Soil organic carbon stocks

The soil data from the 2006 sampling up to 30 cm depth (Ťupek et al., 2015) were combined with additional soil sampling cores of up to 100 cm depth in October 2015 (three per site) (Fig. S1 in the Supplement). The bulk density and C and N concentrations for new samples were determined as in Ťupek et al. (2015).

The SOC content (g cm⁻³) of separate soil layers was interpolated for the whole profile with the fitted spline functions and summed for each depth and each forest and mire site (Fig. S1). The SOC content was similar in the uppermost humus layers of all forest and mire types (below 0.3 g cm⁻³ in a layer 0–10 cm), but in the sub-surface level (10–30 cm), it clearly doubled from uplands to transitions and mires (from below 0.2 to above 0.4 g cm⁻³) (Fig. S1). In the soil layers below 30 cm, the SOC content showed differences in degrees of magnitude (around 0, 0.01, and 0.1 g cm⁻³ for forests, transitions, and mires, respectively) (Fig. S1). The SOC stock (kg C m⁻²) was a result of SOC content multiplied by a bulk density.

2.2.4 Biomass of tree stand and understorey vegetation

The breast height diameter and height of all Scots pine (Pinus sylvestris), Norway spruce (Picea abies), and silver birch (Betula pendula) trees for each forest site type were measured in 2006. The biomass components for each species (leaves, branches, stems, coarse roots) were estimated with biomass conversion functions (Repola, 2008, 2009) and for fine roots with functions by Lehtonen et al. (2016b). Forest floor plants from three 0.07 m² sample plots located near soil respiration measuring collars were harvested for each forest and mire site type in June–July 2004 (Ťupek et al., 2008). Plants were separated into herbs, mosses, and shrubs and dried and weighed for each category and each sample plot. The stand density and the tree biomass increased from xeric (CT) and mesic upland forest sites (VT and MT) towards the herb-rich forest site (OMT) and transitions (OMT+, KgK, and KR) and decreased to very sparse canopy in peatlands and mires (VSR sites) (Fig. 1b). The understorey aboveground biomass correlated negatively with the density of the canopy cover and thus positively with the light intercepted on the forest floor (Ťupek et al., 2008).

2.3 Data analysis

2.3.1 Yasso07 SOC and CO₂ modelling

Equilibrium SOC stocks of up to 1 m depth, SOC changes, and soil CO₂ emissions (R_h) for the forest and mire types were modelled using the Yasso07 soil carbon model (Tuomi et al., 2009, 2011) with specific litter input and weather data in accordance with the method of the Finnish greenhouse gas inventory (Statistics Finland, 2023). The temperature and precipitation data for the weather input were from the nearest Finnish Meteorological Institute (FMI) weather station, located 3 km away from our study sites. We first ran the Yasso07 model using the original formulation of the environmental function with precipitation and air temperature data, and then we ran the Yasso07 model fitted with the environmental modifier function based on the SWC₁₀ and T₅ of the forest and mire site types using the Bayesian data assimilation technique.

Yasso07 is a semi-empirical process-based soil carbon model where soil C is divided based on the organic matter solubility into five pools (C_A, C_W, C_E, C_N, and C_H) from which three are fast (acid, A; water, W; and ethanol, E, all soluble), one is slow (non-soluble, N), and one is almost stable (humus, H) (Tuomi et al., 2011). The rates of C decomposition in each pool and C transfers between the pools are affected by climate. The model can be expressed mathematically as a set of differential equations where decomposition of the entire structural matrix of C pools ${\mathrm{C}}_{\mathrm{A}},\mathrm{\dots },{\mathrm{C}}_{\mathrm{H}}$ defined by default mass flow parameters ${\mathit{\alpha }}_{\text{A,W}},\mathrm{\dots },{\mathit{\alpha }}_{\mathrm{H}}$ and decomposition coefficients $k_{\mathrm{A}},\mathrm{\dots },k_{\mathrm{H}}$ (A_YS) is scaled by the time-step-dependent scalar of the environmental rate modifier ξ(t) (Eq. 1).

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}c\left(t\right)}{\mathrm{d}t}}=\left(\begin{array}{c}{\mathit{i}}_{\mathrm{A}}\\ {\mathit{i}}_{\mathrm{W}}\\ {\mathit{i}}_{\mathrm{E}}\\ {\mathit{i}}_{\mathrm{N}}\\ {\mathit{i}}_{\mathrm{H}}\end{array}\right)\left(t\right)\\ & +\mathit{\xi }\left(t\right)\left(\begin{array}{ccccc}-k_{\mathrm{A}}& {\mathit{\alpha }}_{\text{A,W}}{k}_{\mathrm{W}}& {\mathit{\alpha }}_{\text{A,E}}{k}_{\mathrm{E}}& {\mathit{\alpha }}_{\text{A,N}}{k}_{\mathrm{N}}& \mathrm{0}\\ {\mathit{\alpha }}_{\text{W,A}}{k}_{\mathrm{A}}& -k_{\mathrm{W}}& {\mathit{\alpha }}_{\text{W,E}}{k}_{\mathrm{E}}& {\mathit{\alpha }}_{\text{W,N}}{k}_{\mathrm{N}}& \mathrm{0}\\ {\mathit{\alpha }}_{\text{E,A}}{k}_{\mathrm{A}}& {\mathit{\alpha }}_{\text{E,W}}{k}_{\mathrm{W}}& -k_{\mathrm{E}}& {\mathit{\alpha }}_{\text{E,N}}{k}_{\mathrm{N}}& \mathrm{0}\\ {\mathit{\alpha }}_{\text{N,A}}{k}_{\mathrm{A}}& {\mathit{\alpha }}_{\text{N,W}}{k}_{\mathrm{W}}& {\mathit{\alpha }}_{\text{N,E}}{k}_{\mathrm{E}}& -k_{\mathrm{N}}& \mathrm{0}\\ {\mathit{\alpha }}_{\mathrm{H}}{k}_{\mathrm{A}}& {\mathit{\alpha }}_{\mathrm{H}}{k}_{\mathrm{W}}& {\mathit{\alpha }}_{\mathrm{H}}{k}_{\mathrm{E}}& {\mathit{\alpha }}_{\mathrm{H}}{k}_{\mathrm{N}}& -k_{\mathrm{H}}\end{array}\right)\\ & \times \left(\begin{array}{c}{c}_{\mathrm{A}}\\ c_{\mathrm{W}}\\ c_{\mathrm{E}}\\ c_{\mathrm{N}}\\ c_{\mathrm{H}}\end{array}\right)\left(t\right),\end{array}\end{array}$$

where t is time, i defines a vector of initial carbon pools ${\mathit{i}}_{\mathrm{A}},\mathrm{\dots },{\mathit{i}}_{\mathrm{H}}$, and subscripts to α indicate mass transfer pools (e.g. α_A,W defines mass transfer from pool W to pool A). The total soil respiration or CO₂ efflux (R_h) is a product of a column vector by a row vector C(t), where the elements of the column vector are the fractions that were not transferred among the pools (Sierra et al., 2012).

The model was originally calibrated for running on annual time steps (Tuomi et al., 2009), but it can run on monthly steps with monthly decomposition rates ($\mathrm{1}/\mathrm{12}$ of annual k_A, …, k_H) and monthly litter and climate data (Ťupek et al., 2019). Then ξ(t_m) is defined by a combined function of monthly air temperature (T_m) and $\mathrm{1}/\mathrm{12}$ of annual precipitation ($P_{\mathrm{a}}/\mathrm{12}$) (Eq. 2).

$$\begin{array}{}\text{(2)}& {\mathit{\xi }}_{\mathrm{T}}\left(t_{\mathrm{m}}\right)=e^{\left({\mathit{\beta }}_{\mathrm{1}}{T}_{\mathrm{m}}+{\mathit{\beta }}_{\mathrm{2}}{T}_{\mathrm{m}}^{\mathrm{2}}\right)}\left(\mathrm{1}-e^{\mathit{\gamma }\frac{P_{\mathrm{a}}}{\mathrm{12}}}\right),\end{array}$$

where β₁, β₂, and γ are parameters of the environmental function and t_m is the monthly time step. To test our hypothesis of running the model for a catena of soils with gradually increasing moisture content (from xeric to mesic, paludified, and saturated), we re-defined the ξ(t_m) function for use with soil temperature based on a Q₁₀ exponential function applied to T₅ (used by Davidson et al., 2012, as an alternative to Arrhenius kinetics) and with moisture data using an adjusted hump-shaped Ricker function (Bolker, 2008) for response to SWC₁₀, which limits the decomposition rate outside the optimum soil water content (ξ_AR, Eq. 3).

$$\begin{array}{}\text{(3)}& {\mathit{\xi }}_{\text{AR}}\left(t_{\mathrm{m}}\right)=Q_{\mathrm{10}}^{\left(\frac{T_{\mathrm{5}}-\mathrm{10}}{\mathrm{10}}\right)}a{\mathrm{SWC}}_{\mathrm{10}}{e}^{\left(-a{e}^{-\mathrm{1}}{\mathrm{SWC}}_{\mathrm{10}}\right)},\end{array}$$

where the Q₁₀ parameter represents the increase in the temperature function over a 10 °C difference in T₅ and a controls both ascending and descending slopes of the moisture function when the peak is set to 1. In the Ricker function with a and b parameters and an independent variable vector x (axe^−bx), the height of the peak can be inferred as $(a/b)e^{-\mathrm{1}}$ and the x value of the peak location as $\mathrm{1}/b$ (Bolker, 2008). Thus, in our formulation, by setting the peak in the Ricker function to 1, we could substitute the b parameter with ae⁻¹ and the SWC₁₀ optimum (the SWC₁₀ when decomposition is optimum) was inferred as $\mathrm{1}/a{e}^{-\mathrm{1}}$.

Table 1 The posterior distribution of parameters of the Yasso07 soil carbon model (parameters the same as in Table S1) coupled with the environmental function ξ_AR (Eq. 3, parameters Q₁₀ and a) optimized with observations of SOC stocks (p(θ|SOC)), SOC stocks and CO₂ emissions (p(θ|SOC−CO₂)), and CO₂ emissions (p(θ|CO₂)) from forest–mire ecotone sites using Bayesian data assimilation (Hartig et al., 2012). PSRF stands for the Gelman–Rubin potential scale reduction factor and MAP for a maximum a posteriori probability. The SWC₁₀ optimum (the SWC₁₀ when the SOC decomposition in the boreal forest–mire ecotone was optimum) was inferred as $\mathrm{1}/a{e}^{-\mathrm{1}}$ and ranged between 14 % and 27 % (for a_MAPSOC, $a_{\text{MAP}{\mathrm{SOCCO}}_{\mathrm{2}}}$, and $a_{\text{MAP}{\mathrm{CO}}_{\mathrm{2}}}$).

The Yasso07 model versions in this study run accordingly as follows:

1. The Yasso07.ξ_TW version is Yasso07 coupled with the original ξ_T (Eq. 2) and with the original global parameter set (Tuomi et al., 2011) but with two k-rate parameter sets, (i) the original $k_{\mathrm{A}},\mathrm{\dots },k_{\mathrm{H}}$ rates for application on mineral soils applied to mineral and organo-mineral soil forests (CT, VT, MT, OMT, OMT+, KgK, KR) and (ii) with an inhibitor reducing k rates by 35 % (0.65 $k_{\mathrm{A}},\mathrm{\dots },\mathrm{0.65}{k}_{\mathrm{H}}$) for application on wetlands (Goll et al., 2015; Kleinen et al., 2021) applied to mire sites (VSR1 and VSR2).

2. The Yasso07.ξ_W wetland version is the same as Yasso07.ξ_TW but with a fine-tuned k-rate inhibitor to fit the SOC of mire sites (VSR1 and VSR2).

3. The Yasso07.ξ_AR version for the soil moisture gradient from mineral to peat soils (e.g. as in the boreal forest–mire ecotone) is the Yasso07 model coupled with ξ_AR (Eq. 3), with the original global parameter set of the structural matrix and optimized parameters of ξ_AR.

The initial equilibrium SOC stock (C_o) for each forest and mire type for the pre-trenching period was simulated analytically (Xia et al., 2012) (Eq. 4).

$$\begin{array}{}\text{(4)}& C_o=-\mathit{\xi }{A}_{\text{YS}}^{-\mathrm{1}}\stackrel{\mathrm{‾}}{u},\end{array}$$

where ξ is the environmental modifier, A_YS is a structural matrix formulation of the Yasso07 model's differential equations, and $\stackrel{\mathrm{‾}}{u}$ is the litter input (mean annual litter sum of foliage, branches, stem, stump, roots, and understorey).

The Yasso07 model source code used here was built in R software (R Core Team, 2023) on the platform of the SoilR package (Sierra et al., 2012) according to the mathematical description and parameters of Tuomi et al. (2011). The model outputs are monthly SOC stocks and soil CO₂ emissions. The model was run with data inputs of above- and belowground litterfall (accounting for its chemical composition) and climate data (described in more detail below). Monthly model outputs of heterotrophic soil respiration were compared to mean monthly R_h measurements.

2.3.2 Climate and litter C input data for the Yasso07 model

Yasso07.ξ_TW was run with monthly air temperature and precipitation from the nearby Juupajoki weather station of the Finnish Meteorological Institute. Yasso07.ξ_AR was run using site-type-specific continuous monthly T₅ and SWC₁₀ time series.

Figure 2 Observations specific to the forest or mire site type of soil organic carbon (SOC) stocks (kg C m⁻², summed to 1 m) and the mean volumetric soil water content (SWC) at 10 cm depth (%) (with error bars showing 1 standard deviation of all measured values) (a) in comparison to their distributions of heterotrophic soil CO₂ emissions/respiration measurements (R_h, $\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) and R_h expressed as the emitted C fraction per site-specific SOC stock ($\mathrm{\mu }\mathrm{g}\mathrm{C}{\mathrm{g}}^{-\mathrm{1}}\mathrm{SOC}{\mathrm{h}}^{-\mathrm{1}}$) (b). The CT, VT, MT, and OMT types represent upland forests; OMT+, KgK, and KR forest–mire transitions; and VSR1 and VSR2 mires. The boxplot horizontal lines show the 25th and 75th interval with the median in between, and the 5th and 95th confidence interval is shown by whiskers.

The litter C input of the forest and mire types (Figs. S2 and S3) used by Yasso07 was estimated as in Lehtonen et al. (2016a) based on turnover rates of tree stand biomass components (including fine and coarse roots, stump, branches, and foliage) and understorey vegetation. The litter C input was separated into Yasso07 A, W, E, and N pools according to the component and species-specific (or species-group-specific) A, W, E, and N ratios taken from the literature (Berg et al., 1991a, 1b, 1993; Gholz et al., 2000; Trofymow, 1998; Vávřová et al., 2009; Straková et al., 2010) (Table S2 in the Supplement). The annual litter was distributed either according to a monthly resolution by accounting for seasonal trends in foliage, fine roots, and the understorey (Ťupek et al., 2019; Zhiyanski, 2014) or evenly. The litter input before trenching was assumed to represent the long-term average of the equilibrium state forest (Figs. S2a and S3). During trenching, the severed fine and coarse roots made up the major component of the total litter (Fig. S2b) and resulted in a clear peak in the monthly litterfall time series (Fig. S3). After trenching, the monthly litter levels decreased as the sum of components excluded the roots (Figs. S2c and S3).

2.3.3 Bayesian SOC and CO₂ data assimilation

The Bayesian Markov chain Monte Carlo (MCMC) data assimilation has proved useful in improving soil organic carbon estimates (e.g. Xu et al., 2006; Hararuk et al., 2014). The Bayesian posterior uncertainty provides updated information on parameter values based on pre-existing information on the parameters and the data through the likelihood function (Speich et al., 2021). The Q₁₀ and parameters of ξ_AR (Eq. 3) coupled with the Yasso07 model were optimized on the level of the forest–mire ecotone using the Bayesian data assimilation technique (Luo et al., 2011; Hartig et al., 2012; Speich et al., 2021) with observed SOC stocks and monthly R_h data of forest and mire types with prior information on the best parameter values and the defined parameter range in Table S1. During the optimization, the Yasso07.ξ_AR model was run first only with observed SOC stocks, second with both SOC stocks and CO₂ data combined, and third only with CO₂ data obtaining a probability distribution of model parameters of ξ_AR (the posterior uncertainty p(θ|y) conditional on the observations (y) and prior knowledge on the parameter values p(θ)). The sum of the probability density for the target parameter set (θ) between the model predictions and observations was maximized for the best agreement using the likelihood defined by a modified Laplace probability density function p(y|θ) (the probability of observing the data y with the model parameters set θ), where we allowed the width of the distribution to be affected by the observed SOC and CO₂ values (Eq. 5).

$$\begin{array}{}\text{(5)}& p\left(y\mathrm{|}\mathit{\theta }\right)={\prod }_{j=\mathrm{1}}^{N_j}{\prod }_{i=\mathrm{1}}^{N_i}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}({\mathit{ϵ}}_{{\mathrm{a}}_j}+{\mathit{ϵ}}_{{\mathrm{b}}_j}{x}_{j,i})}}{e}^{\frac{-\left|\left(x_{j,i-}{\mathit{\mu }}_{j,i}\right)\right|}{{\mathit{ϵ}}_{{\mathrm{a}}_j}+{\mathit{ϵ}}_{{\mathrm{b}}_j}{x}_{j,i}}},\end{array}$$

where μ_j,i is the observed jth variable (e.g. SOC, CO₂, or SOC and CO₂) of ith observations, x_i is the modelled prediction, N is the total number of observations, and ϵ_a and ϵ_b are parameters affecting the width of the distribution. In the combined SOC–CO₂ likelihood, the likelihood function p(y|θ) was then the multiplication of the distributions of SOC and CO₂ at all observation times. We evaluated the variation in the estimated parameters by separating data for fitting the models and testing them with a 9-fold cross-validation technique.

Table 2 The SOC and CO₂ performance statistics of Yasso07 (Y07) model versions against the measured data in the boreal forest–mire ecotone, where MAE is the mean absolute error, MBE is the mean bias error, RMSE is the root-mean-square error, $R_{\text{adj}}^{\mathrm{2}}$ is the adjusted coefficient of determination, and AIC is the Akaike information criterion. The units of MAE, MBE, and RMSE are in kg C m⁻² and kg CO₂ per square metre per month for SOC and CO₂, respectively.

The model parameters of ξ_AR and p(y|θ) were sampled from an assumed uniform distribution within their prior ranges (Table S1). Posterior probability distributions of parameters (Table 2) were derived using the differential evolution “DEzs” MCMC sampler used by the runMCMC function from the BayesianTools package in R (Hartig et al., 2019) and by computing three chains in parallel. The convergence of MCMC runs was evaluated using the Gelman–Rubin multivariate potential scale reduction factor (PSRF) (Brooks and Gelman, 1998). The MCMC simulation was considered converged if PSRF was below 1.03 for all parameters (1.01, 1.01, and 1.016 for p(θ|SOC), p(θ|CO₂), and p(θ|SOCCO₂), respectively). Trace plots of MCMC runs for target parameters showed effective sampling and unimodal parameter density with clearly defined peaks. The differences in parameter uncertainties (difference between 97.5 % and 2.5 % quantiles of the 95 % confidence interval) were not significant (p = 0.99) when evaluated with a Welch two-sample t test between two posterior distributions p(θ|SOC) and p(θ|SOCCO₂) (Table 1).

2.3.4 Performance evaluation of Yasso07.ξ_TW and Yasso07.ξ_AR

The performance of Yasso07 model versions (i) Yasso07.ξ_TW and (ii) Yasso07.ξ_AR with ξ_AR parameter set θ from three posterior distributions, p(θ|SOC), p(θ|CO₂) and p(θ|SOCCO₂), was evaluated with the modelled SOC and CO₂ outputs against the observed data in the forest mire ecotone with the coefficient of determination (R²), the mean absolute error (MAE), mean bias error (MBE), the root-mean-square error (RMSE), and the Akaike information criterion (AIC), considering the number of model parameters in the error calculation as in Abramoff et al. (2022) and with the fitted linear trends of normalized SOC and CO₂ model residuals with observations against T₅ and SWC₁₀ data.

3 Results

3.1 Distributions of SOC stocks and R_h in relation to SWC

The SOC stock measurements (to a depth of 1 m) in a forest–mire ecotone were distributed in range between 20 kg C m⁻² in well-drained soils of upland forests and 125 kg C m⁻² in poorly drained soils in peatlands and mires (Fig. 2). The SOC stock values strongly correlated with the long-term moisture levels. The median R_h values ranged between 0.4 and 0.6 $\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ for upland forests, 0.4 and 0.5 $\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ for forest–mire transitions, and 0.3 and 0.4 $\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ for mires (Fig. 2). The forest and mire site type differences in median R_h values expressed per square metre were small and poorly correlated with the mean soil moisture levels.

Figure 3 The 3-year time series (2004 – wet; 2005 – typical; 2006 – dry) of instantaneous measurements of (a) soil heterotrophic respiration (R_h, $\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$, positive sign), (b) soil temperature at 5 cm depth (°C), and (c) soil moisture at 10 cm depth (%) of nine forest and mire sites (four upland forests (CT, MT, and OMT), three forest–mire transitions (OMT+, KgK, and KR), and two mires (VSR1 and VSR2)). The sites are arranged from left to right according to their position on the slope (see Fig. 1).

3.2 Distribution of R_h in climate space of soil T and SWC

The site-specific time series of hourly R_h measured in the forest–mire ecotone during the years 2004, 2005, and 2006 followed a typical seasonal pattern of temperature and was distributed in range between 0.08 and 1.6 $\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ depending on the corresponding soil temperature and moisture conditions (Fig. 3). The R_h values were generally larger during wet years than during a typical year and were lowest during dry years (Fig. 3).

The T₅ and SWC₁₀ values showed a typical seasonal variation (in a range between around 0 and 20 °C, driest in the summer and wettest in late autumn–spring) (Figs. 2b and 3c). T₅ showed a similar magnitude among the forest and mire sites, whereas the SWC₁₀ increased from the driest (upland forest) to intermediate (forest–mire transition) sites and from upland to lowland for the wettest (mire) sites located on the slope (Fig. 3). The volumetric SWC₁₀ values (%) were comparable and on the same order of magnitude between the forest and mire site types because the forest–mire variation in the soil organic layer bulk density was relatively small at 0.34 ± 0.07 g cm⁻³ (porosity 74 ± 5 %). The forest–mire ecotone soil moisture at 10 cm depth ranged from 5 % to 91 %. The minimum, maximum, and mean SWC at 10 cm depth between forests, transitions, and mires clearly differed, showing a gradient of increasing moisture from forests to mires (Figs. 2a and 3). Due to highly variable weather (wet, typical, and dry years), all ecosystems experienced periods of extremely low and high SWC₁₀ values. The SWC₁₀ ranged between 5 % and 25 % in upland forest, between 17 % and 70 % in transitions, and between 49 % and 91 % in mires (Fig. 3). The variation in soil temperature at 5 cm depth along the ecotone was similar among the forest and mire types and ranged between −3 and 22 °C (Fig. 3). The R_h values during the dry year of 2006 were, in comparison to previous years, clearly reduced, mostly in upland forest and forest–mire transitions (Fig. 3). The spatial soil moisture gradient of forest and mire types was not clearly reflected in the distribution of R_h values when expressed in $\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ (Fig. 2) (unless expressed as a C fraction of SOC). The short-term SWC variation impacted the typical seasonal levels of R_h values, mainly during the extreme events (the rainy summer period in wet years or drought summer period in dry years) (Fig. 3).

Figure 4 The optimized environmental modifier of default decomposition rates ξ_AR (Eq. 3) (coupled with the Yasso07 model) drawn with mean posterior values of parameters and their confidence intervals (dashed lines) (Table 1) for separate responses to (a) soil temperature at 5 cm, ξ_AR = f(T₅) when f(SWC₁₀) = 1 and (b) soil water content at 10 cm, ξ_AR = f(SWC₁₀) when f(T₅) = 1. The functions were fit based on only CO₂, SOC and CO₂, or only SOC data.

Figure 5 The colours and contour lines showing the optimized environmental modifier of default decomposition rates ξ_AR (Eq. 3) (coupled with Yasso07 model) drawn with mean posterior values of parameters (Table 1) for combined responses to soil temperature at 5 cm, ξ_AR = f(T₅), and to soil water content at 10 cm, ξ_AR = f(SWC₁₀), based on only SOC (a), SOC and CO₂ (b), or only CO₂ (c) data. In the panels of combined ξ_AR white circles show pairs of corresponding monthly means of T₅ and SWC₁₀, and the black circles show the annual T₅ and SWC₁₀ means for 9 forest and mire site types.

3.3 ξ_AR optimized with Yasso07.ξ_AR

The optimization of ξ_AR (Eq. 3) coupled with Yasso07 showed that in the catena of mineral and organic soils of the boreal forest–mire, the optimum moisture content for decomposition and CO₂ emissions was in well-drained mineral soil forests (SWC_opt medians between 14 % and 27 %, Table 1, Fig. 4). The MCMC fit with CO₂ data produced larger SWC_opt and a larger tail in the Ricker function (compared to the SOC or SOCCO₂ fit). The decomposition rate outside the moisture optimum reduced decomposition similarly for the two data sources (SOC and SOCCO₂) used for calibration. However, the temperature and moisture functions were different when only CO₂ was used for calibration.

The two SOC- and SOCCO₂-based Q₁₀ functions showed a similar increase with T₅ until 10 °C (Fig. 4a). The combined non-linear temperature and moisture response in the whole climate data range showed larger change in ξ for mineral soil forests compared to forest mire transitions and peatlands (Fig. 5). The ξ_AR values in Fig. 5a and b for SOC- and SOCCO₂-based ξ_AR are similar, showing that both SOC parameterization and SOCCO₂ parameterization are almost the same, whereas the ξ_AR in Fig. 5c is different. The CO₂-based temperature sensitivity of the median ξ_AR Q₁₀ value of 4.5 resulted in a more pronounced increase in decomposition rates, especially in climate space with high water saturation (Figs. 4 and 5c). The SOC- and SOCCO₂-based ξ_AR Q₁₀ values for the forest–mire ecotone (medians 2.3 and 2.5, Table 1) were lower.

3.4 Performance of Yasso07.ξ_TW, Yasso07.ξ_W, and Yasso07.ξ_AR

The model performance evaluation showed that the soil water and temperature modifier ξ_AR coupled with the Yasso07 model (Yasso07.ξ_AR) outperformed the original Yasso07 environmental function even after a 65 % reduction in decomposition rates for wetlands was applied (Yasso07.ξ_TW) (Table 2, Fig. 5). Although the Yasso07.ξ_TW model version accurately predicted SOC stocks of mineral soil forests (CT, …, OMT), it heavily underestimated the SOC stocks of organo-mineral forested peatlands and mires (OMT+, …, VSR2); thus it showed the most biased model performance metrics (highest RMSE, MBE, MAE, and AIC and lowest $R_{\text{adj}}^{\mathrm{2}}$) among the model versions compared (Table 2). A reduction in decomposition rates of 65 % for mires in Yasso07.ξ_TW was not sufficient to simulate their SOC stocks as the simulated SOC of mires was only about 10 % of measured values (Fig. 5a–c). The SOC simulations for VSR mires with Yasso07.ξ_W would have required as much as a 96 % reduction in the decomposition rates. The optimized Yasso07.ξ_AR model version accurately simulated SOC stocks throughout the forest–mire ecotone.

The Yasso07.ξ_AR,p(θ|SOC) version outperformed Yasso07.${\mathit{\xi }}_{\text{AR},p\left(\mathit{\theta }\mathrm{|}{\mathrm{SOCCO}}_{\mathrm{2}}\right)}$ when evaluated against SOC data, and both models were similar when evaluated against CO₂ data (Table 2, Fig. 5). Yasso07.${\mathit{\xi }}_{\text{AR},p\left(\mathit{\theta }\mathrm{|}{\mathrm{CO}}_{\mathrm{2}}\right)}$ outperformed the functions based on SOC or SOCCO₂ against measured CO₂ data but failed when evaluated against measured SOC values. The Ricker function improved the representation of decomposition for drier soils and the representation of optimal SWC for decomposition.

Figure 6 Scatterplots between observed SOC (kg C m⁻²) and R_h ($\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) from the forest–mire ecotone against modelled values with the two versions of Yasso07 model: (i) Y07.ξ_TW – Yasso07 coupled with the default environmental modifier (ξ_T, Eq. 3) based on air T and precipitation with a global parameter set (Tuomi et al., 2011) (applied for CT, …, KR mineral and organo-mineral soil forest sites) and with a reduction in decomposition rates by 65 % for wetlands (Goll et al., 2015; Kleinen et al., 2021) (applied to VSR1 and VSR2 mires sites) – and (ii) Y07.ξ_AR – Yasso07 coupled with an environmental modifier (ξ_AR, Eq. 3) based on SOC, SOC–CO₂, or CO₂ data (a, d) compared with the 1:1 line (dashed red line). The normalized model residuals (residuals divided by observations) are plotted against T₅ and SWC₁₀ with the trend lines of the linear fits and with their confidence intervals (dashed lines) (b, c, e, and f).

The SWC optimum was derived from the fitted ascending slope parameter, and its values were between 14 % and 27 % (depending on the data used for fitting; 14 % for SOC and CO₂ and 27 % for CO₂). The MCMC fit with CO₂ data produced larger SWC_opt and a larger tail in the Ricker function (compared to the SOC or SOCCO₂ fit). However, the CO₂-only fit also underestimated SOC stocks of forested peatlands (Fig. 5a). The normalized SOC residuals of the two Yasso07.ξ_AR models based on SOC or SOC and CO₂ did not show any large T₅ or SWC₁₀ trends (Fig. 5a–c), although measured R_h fluxes during the cold season may be slightly overestimated.

The soil CO₂ emissions simulated with the original Yasso07.ξ_TW agreed unexpectedly well with observed R_h values (Table 2, Fig. 6d–f), outperforming the Yasso07.ξ_AR,p(θ|SOC) version (Table 2, Fig. 6e). On the other hand, the R_h simulated with Yasso07.ξ_AR,p(θ|SOC) performed similarly to Yasso07.ξ_TW in terms of RMSE (the same RMSE of 0.16 $\mathrm{g}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ for both models). The performance statistics for Yasso07.ξ_AR based on SOC and SOCCO₂ were better compared with those of Yasso07.ξ_TW for SOC and about the same for R_h (Table 2). The performance statistics for Yasso07.ξ_AR based on CO₂ were best among the models for CO₂ but failed for SOC in a similar fashion to Yasso07.ξ_TW. The normalized modelled residuals showed that both Yasso07 model versions (Yasso07.ξ_TW and Yasso07.ξ_AR,p(θ|SOC)) showed small R_h biases in extreme, very cold and very warm, temperatures (Fig. 6e). The normalized CO₂ residuals of the Yasso07.${\mathit{\xi }}_{\text{AR},p\left(\mathit{\theta }\mathrm{|}{\mathrm{CO}}_{\mathrm{2}}\right)}$ showed no bias among the functions (0 MBE, Table 2). The normalized CO₂ residuals evaluated against SWC₁₀ did not show any bias for any of the models (Fig. 6f).

4 Discussion

The Yasso07 model (Tuomi et al., 2011), coupled with a revised and optimized empirical Q₁₀ soil temperature and Ricker moisture function ξ_AR (Eq. 3, Fig. 4), successfully reconstructed observed variation in SOC stocks and soil heterotrophic CO₂ emissions with increasing soil wetness in mineral, organo-mineral, and peat soils in boreal forest (Fig. 6). The original Yasso07 monotonic precipitation function is effective due to easily available data on upper boundary conditions, but it is also flawed in the case of a shallow water table when the lower boundary is equally important in defining the water content on the soil. Therefore, the usage of soil water content as a variable is structurally superior, and this can be proved by inductive reasoning, e.g. from the test model runs. Separating the effect of structure from calibration would require more test runs with data from a larger number of study sites.

Our application of Yasso07 models on the hillslope accounted for the continuity in moisture conditions, which was reflected in the modelled gradient of mineral and peat soil carbon stocks. The Yasso07 model initially developed for mineral soils was improved for application in peatlands by accounting for the soil temperature and volumetric moisture, as these are better predictors of heterotrophic respiration than air temperature and precipitation (Jian et al., 2022). Although the empirical Ricker function ξ_AR used here was heuristic, its form implicitly accounted for prevailing intrinsic micro-scale processes on the hillslope controlling R_h and SOC accumulation, e.g. plant and microbial communities and the long-term and short-term limitation of oxygen and substrate with moisture (Davidson et al., 2012; Moyano et al., 2012; Ghezzehei et al., 2019).

ξ_AR being able to simulate gradually increasing SOC stocks from mineral to organic soils makes it a preferable rate modifier for the Yasso07 model, instead of simply adjusting decomposition with a reduction constant for wetlands (e.g. Goll et al., 2015; Kleinen et al., 2021), which underestimated the SOC stocks of peatlands (Yasso07.ξ_TW in Fig. 5). In this study, a constant 96 % reduction (0.04 × k rates) was proposed for the existing Yasso07.ξ_TW model for more accurate SOC modelling in mires, a value comparable to rates of anaerobic decomposition (Schuur et al., 2015). The 96 % reduction is comparable to JULES, which accounts for oxygen inhibition with a gradual reduction in decomposition from the maximum rate of 1 at the moisture optimum (30 %–75 % SWC) to a reduced rate of 0.2 in water-saturated peat soils (Chadburn et al., 2022).

4.1 The moisture response

The use of gradually changing near-surface soil moisture avoids biases in land surface modelling related to ignoring high SOC stocks of organo-mineral soils of forested peatlands (Dalsgaard et al., 2016, Ťupek et al., 2016), e.g. forest–mire transitions (Figs. 1 and 2). The accurate modelling of decomposition rates with diffusion-based moisture functions accounting for microbial processes requires the correct representation of the drivers of heterotrophic respiration (Yan et al., 2018; Moyano et al., 2012, 2013; Manzoni et al., 2012; Ghezzehei et al., 2019). However, there is uncertainty in functional moisture – soil respiration dependencies are high (Sierra et al., 2015; Falloon et al., 2011) and dependencies vary with the soil properties; e.g. the SWC optimum increases for soils with higher organic C content (from 30 % to 75 % SWC; Moyano et al., 2012, 2013). The ξ_AR function's SWC_opt found in dry and well-drained conditions and the reduction in default decomposition rates (k) with increasing soil wetness contrasted responses from short-term laboratory incubation soil respiration studies (weeks, months) that showed an increase in decomposition from dry conditions until reduction in very wet conditions (Sierra et al., 2017; Moyano et al., 2012, 2013; Kelly et al., 2000; Skopp et al., 1990; Yan et al., 2018). The ξ_AR optimized with SOC and CO₂ data showed that the optimum/maximum decomposition rate in the forest–mire ecotone in dry, well-drained conditions was around 14 % of mean long-term near-surface SWC (around 20 % water-filled pore space (WFPS), corresponding to sub-xeric and mesic forest site types) (SWC_opt parameters inferred from a parameter in Table 1, Fig. 4b), whereas the moisture optimum of studies based only on respiration from laboratory soil incubations was around 40 %–60 % (Fairbairn et al., 2023; Moyano et al., 2013; Kelly et al., 2000; Skopp et al., 1990; Yan et al., 2018).

The moisture optimum derived from the field site soil respiration datasets from a larger moisture range was found for 50 % WFPS and corresponded to around 31 % SWC assuming a mean porosity of 62 % (Hashimoto et al., 2011). Our SWC_opt between 14 % and 27 % SWC (Table 1) was comparable to the optimum derived from the field site data, which was lower compared to laboratory incubations. The SWC_opt discrepancy of the ξ_AR function highlights the difference between (1) the responses from the field-based or long-term soil respiration measurements reflecting moisture responses of older, stabilized, and slowly decomposing SOC and (2) the short-term incubation-based soil respiration studies which predominantly capture the decomposition of the newly available, labile, and rapidly decomposing SOC pool (González-Domínguez et al., 2019; Huang and Hall, 2017). Over longer periods of incubation, high Q₁₀ can be observed (Zhou et al., 2014). The enhanced C mineralization can occur during periods of elevated moisture under Fe reduction when microbes can access previously protected labile C (Huang and Hall, 2017). The incubations are short term (from a few days to a few months) and are useful to identify short-term processes. Moreover, they are performed on disturbed soils (sometimes even sieved), and therefore the soil structure is not representative of the field.

The ecosystem-scale application of moisture reduction functions obtained in the laboratory can be hindered by several factors. There are a number of feedback mechanisms which modify the response obtained on a limited size soil sample. Among them is a change in microbial community composition, the texture- and structure-dependent effect of pore-scale connectedness of soil solutions, and competition between plants and microorganism for resources under different environmental stress conditions. Under a changing climate, these feedback mechanisms may lead to the system behaving unpredictably in relation to extrapolation. Therefore, the validation of the models at the site level with a series of various in situ stress levels is necessary for reliable future predictions.

In their impact on decomposition, the ξ_AR functions (calibrated with SOC, SOCCO₂, and CO₂ data) incorporated into the Yasso07 soil C model were comparable (e.g. all found the moisture optimum in well-drained soils of the forest–mire ecotone). Although the soil temperature and moisture functions showed relatively small differences in Q₁₀ between the model fits, the a parameter of the moisture functions of the CO₂-based fit was larger than those from the SOC and SOCCO₂ fits (Table 1). In terms of the Yasso07 model constants, if temperature and moisture conditions are favourable for organic matter stabilization, then ξ_AR is reduced (Figs. 4 and 5), which in turn reduces decomposition rates of fast and slow C pools, reduces their CO₂ emissions, and increases C storage. The forest–mire sites' heterotrophic respiration per unit of area did not show a clear difference between well-drained and water-saturated soils, whereas the C mineralization per unit of SOC was clearly reduced in soils with mean long-term field soil moisture (Fig. 2). A reduction in decomposition rates in the environmental gradient from low to higher field moisture indicates a possible difference in the soil C stabilization mechanisms under low and high water content (Das et al., 2019). Ghezzehei et al. (2019) suggested that empirical moisture sensitivity curves should be calibrated individually for each soil type. However, our study shows that the common modifier function, based on the SWC of the topsoil humus layer, which has comparable properties across the soil types, could provide insights into a more generalizable moisture sensitivity function. The mechanistic diffusion-based moisture functions (e.g. by Ghezzehei et al., 2019; Yan et al., 2018) could, in follow-up studies, be compared against deterministic moisture functions to evaluate their applicability and interpretation.

The ξ_AR function with its reduction with increasing wetness from dry soils was based on a large range of forest and mire soil C stocks (between 11 and 134 kg C m⁻²), reflecting a spatial long-term moisture gradient between forests and mires (Fig. 2) and short-term moisture and CO₂ dynamics over years with contrasting climate (Fig. 3). The soil respiration data from 3 years covered exceptionally contrasting wet and dry summers and likely captured a full range dependency on the soil moisture induced by short-term weather variation in a spatial long-term forest–mire gradient for soil moisture, soil C pools, vegetation litter input dynamics (Figs. S2 and S3), and microbial composition. The short-term deviations in respiration indicative of wetting–drying cycles (Barnard et al., 2020; Patel et al., 2021) could be seen by the respiration increases in wet summers or during and after a period of drought (Fig. 3). Thus, the ${\mathit{\xi }}_{\text{AR},p\left(\mathit{\theta }\mathrm{|}\mathrm{SOC}\text{-}{\mathrm{CO}}_{\mathrm{2}}\right)}$ curve calibrated with highly variable SOC and CO₂ data from a forest–mire ecotone represented a mean robust moisture decomposition dependency smoothing short-term weather-dependent fluctuations with the spatial variation in organic matter decomposition across ecological gradients. This function could meet the land surface modelling criteria for not only spatial accuracy on small scales but also cost efficiency for running or forecasting the C dynamics at large scales (Luo and Schuur, 2020).

Including SOC data or a combination of SOC and CO₂ data in model fitting resulted in lower SWC_opt, and the model fitting based only on CO₂ showed larger SWC_opt and a larger tail (descending slope) of the Ricker moisture function. Thus, in comparison to other studies, whose dependencies were limited to relatively short-term responses of only soil heterotrophic CO₂ respiration from mainly mineral soils in laboratory conditions, the differences in SWC_opt observed in our studies could be expected from differences in the data source used in model calibration. Unlike the data from controlled laboratory conditions, we used data from field measurements (mineral soil and peat SOC stocks, litter input, soil CO₂ respiration, T₅, and SWC₁₀ measured under extreme weather variability for 3 years). In optimizing model performance with a multi-variable dataset, Keenan et al. (2013) found that a combination of data with fast and slow turnover (e.g. soil respiration and soil carbon stocks) leads to the largest improvement in model performance. Yasso07.ξ_AR based only on slow turnover (SOC) was as good as that constraining with SOC and CO₂, as both approaches accurately observed soil CO₂ emissions and SOC stocks across the site types of the forest–mire ecotone with no clear bias in residuals (Fig. 6). Thus, in a catena of mineral and peat soils of the forest–mire ecotone and in the combined measured SOC and CO₂ data assimilation in ξ_AR (9 and 2369, respectively), the relatively small number of SOC stocks (nine forest and mire types) largely determined the SWC response form, reflecting both a spatial moisture gradient and its temporal variation. Whether the deterministic modifier rate was estimated correctly or not for the drained peatlands as well should be tested in follow-up studies, as our data did not include drained peatlands. The Ricker functional dependency has performed well for the drier region, but the performance in soils with high water status could still be improved. This could be deduced from the better statistical performance of the CO₂-only fit with CO₂ data (compared to the SOC or SOCCO₂ fit), which produced a larger tail of the Ricker function. The CO₂-only fit also underestimated SOC stocks of forested peatlands.

Although the moisture representation of the ξ_AR environmental function was accurate at the forest–mire ecotone level, at the forest site level, the contrasting respiration responses to moisture (i.e. either respiration reduction during soil drying or increased CO₂ emissions with rewetting for dry soils or the opposite for wet soils; Barnard et al., 2020; Patel et al., 2021) were likely not captured sufficiently. Soil C modelling might be further improved using a moisture response that accounts separately for microbial respiration with increased water availability and for oxygen limitation in soil reaching water saturation (Sierra et al., 2015). However, as the aim of the environmental modifier used in this study was to apply the above concepts in a cost-efficient way using an empirical function with easily interpretable parameters (Q₁₀ and a which informs us about the SWC optimum) (Davidson et al., 2012), the mathematical representation of the moisture function with increased complexity still needs to be evaluated in further studies testing different functional forms with larger regional data availability. Ghezzehei et al. (2019) suggested that empirical moisture sensitivity curves should be calibrated individually for each soil type. The mechanistic diffusion-based moisture functions (e.g. by Moyano et al., 2013; Yan et al., 2018; Ghezzehei et al., 2019) could, in follow-up studies, be tested against deterministic moisture functions (e.g. as in Davidson et al., 2012) to evaluate their applicability and interpretation. However, our study shows that the common modifier function, based on the SWC of the topsoil humus layer, which has comparable properties across the soil types, could provide insights into a more generalizable moisture sensitivity function.

In this study, we constrained the soil carbon model using both SOC (stock) and CO₂ (flow) data. Few studies have constrained the soil carbon model to both SOC and CO₂ data. Our study demonstrates the importance of extensive constraints on the soil carbon model to obtain a reliable model output. The SOC constraint improved the model performance; at the same time, the intensive SOC and CO₂ constraint did not result in the improvement of model performance, which implies the need for further model development and testing. One potential improvement in modelling could be to account for the different responses to the environment (e.g. soil moisture) among different pools, like the temperature dependency separated between the soil layers and soil C fractions in more recent versions of the Yasso model, e.g. Yasso15 and Yasso20 (Viskari et al., 2020, 2022). The Yasso07 model adapts one common response function among different pools for simplicity; however, the fresh plant litter moisture limitation of decomposition may be expected to differ from the moisture limitation on older stabilized C in the humus horizon and mineral-associated C. Another factor could be the vertical process. The SOC is vertically distributed in the soil, and soil C fractions differ among soil depths. Accounting for the depth of the soil layer with the largest proportion of net CO₂ emissions (Davidson et al., 2006, 2012), which is expected to vary with fluctuating water levels in forested peatlands, may further improve the soil respiration estimates for organic soils. On a process level, the key to understanding the difference in the moisture reduction function at different soil depths may lie in the nature of the physical and biochemical availability of substrate to enzymes released by microbial decomposers (Sainte-Marie et al., 2021).

4.2 The temperature response

The original air-temperature-based modifier in Yasso07 was replaced by the Q₁₀-type temperature function driven by soil temperature (Davidson et al., 2012). This function was found to best represent the enzyme kinetics under unconstrained substrate and oxygen (Sierra et al., 2017). The optimized temperature function with SOC and combined SOC–CO₂ data produced slightly more biased R_h estimates than the modifier based on CO₂ data with Q₁₀ around 4.5 (Table 1, Fig. 6). The Q₁₀ values around 4.0 were comparable with the well-drained, moderately drained, and poorly drained forest soils for similar climates (Chen et al., 2020; Davidson et al., 1998; Karhu et al., 2010; Pumpanen et al., 2008). However, the optimization of the Q₁₀-type temperature response only with SOC data and with SOC and CO₂ (Table 1) showed lower $R_{\text{adj}}^{\mathrm{2}}$ values than for the Gaussian type of temperature dependency in Yasso07.ξ_TW (Fig. 5). Thus, due to the comparable predictive power when using soil and air temperature (Jian et al., 2022), the original Gaussian air temperature dependency could be more accurate than the Q₁₀ response for the optimization with soil temperature (Tuomi et al., 2008).

Discrepancies in modelled respiration during winter (Fig. 6) could also be caused by a scarcity of winter field CO₂ measurements, potentially resulting in larger random errors (e.g. due to difficulties of measuring relatively small respiration fluxes during soil freezing–thawing cycles, measurements on soil covered by a snow layer, and reduced precision of gas analysers during measurements in a lower temperature range). The density of CO₂, temperature, and SWC measurements can be seen in Fig. 3. The less frequent measurements during periods of near-zero soil temperature might have affected the fit of the temperature function. However, our main emphasis was on the moisture, which in near-zero temperature conditions plays only a minor role in controlling respiration.

5 Conclusions

The Yasso07 soil carbon model was developed and parameterized at a global scale for mineral soils; however, it has also been applied to land surface modelling coupled with the JSBACH model with a 65 % reduction in default decomposition for wetlands. In this study we emphasized the improving representation of the response of soil organic C stock change and respiration to soil moisture in the Yasso07 model for forest–mire ecosystems. At the forest site level, we evaluated the performance of the Yasso07 model with an original climate modifier based on air temperature and precipitation against the model coupled with a revised environmental modifier based on soil temperature and moisture. We found that the Yasso07 model coupled with revised climate dependencies performed similarly for mineral soils but outperformed the original configuration with the JSBACH modification for undrained peatland soils.

The optimization of moisture dependency conducted in this study accounted for both a spatial moisture gradient and its temporal variation. The moisture optimum for dry soils has not changed depending on whether the function was optimized using both slow (SOC) and fast (CO₂) turnover data (combined SOC and CO₂ data) or only slow (SOC) data.

The SOC stocks in peatland forests were an order of magnitude larger in comparison to forests on mineral soil. On a landscape level, these peatland SOC stocks had the largest influence on calibrating the moisture optimum, when they were included along with fluxes in optimization. This could be inferred from the same calibrated moisture optimum when using calibration with only SOC or SOCCO₂ as the data source, whereas for only CO₂-based calibration, the optimum differs. The function implicitly accounted for a relative contribution of C fluxes from short-term biogeochemical processes in long-term SOC accumulation. For accurate estimates of the boreal forest soil carbon pools with the Yasso07 model, the SOC accumulation related to inhibition of decomposition with increasing wetness was more pronounced than the one related to dryness.

This study illustrated the limitation of the default moisture functions used for peatland forest soil C modelling. Also, the non-monotonic Ricker function with a moisture optimum in well-drained mineral soils needs further evaluation with regional boreal forest data. The exact representation of the functional form of the soil moisture dependency is characteristic of conditions of our study, e.g. the distribution of organic and mineral soil forests in the data. Broader extrapolation of the conclusions, e.g. regarding climate change or forest management on drained peatlands, would require more model testing with spatially larger data and lower water levels in forests on organic soils. However, if the soil moisture optimum of litter decomposition in forests on well-drained mineral soils of boreal landscapes proves to be robust, then in the future warmer and drier climates, the boreal forest could be expected to enhance soil C emissions to the atmosphere due to water level drawdown of presently water-saturated peat soils with large C stocks. In contrast, rewetting of previously drained peatlands could be expected to reduce soil C emissions, turning SOC loss to long-term C sequestration.