Impact of vegetation composition and seasonality on sensitivity of modelled CO2 exchange in temperate raised bogs

Encroachment of vascular plants (VP) in temperate raised bogs, as a consequence of altered hydrological conditions and nutrient input, is widely observed. Effects of such vegetation shift on water and carbon cycles are, however, largely unknown and identification of responsible plant physiological traits is challenging. Process-based modelling offers the opportunity of gaining insights into ecosystem functioning beyond observations, and to infer decisive trait shifts of plant functional groups. We adapted the Soil–Vegetation–Atmosphere Transfer model pyAPES to a temperate raised bog site by calibration against measured peat temperature, water table and surface CO2 fluxes. We identified the most important traits determining CO2 fluxes by conducting Morris sensitivity analysis (MSA) under changing conditions throughout the year and simulated VP encroachment. We further investigated transferability of results to other sites by extending MSA to parameter ranges derived from literature review. We found highly variable intra-annual plant traits importance determining ecosystem CO2 fluxes, but only a partial shift of importance of photosynthetic processes from moss to VP during encroachment. Ecosystem respiration was dominated by peat respiration. Overall, carboxylation rate, base respiration rate and temperature sensitivity (Q10) were most important for determining bog CO2 balance and parameter ranking was robust even under the extended MSA.

Bog peatlands are characterised by a permanent water saturation and low pH values, hampering organic matter from decomposition 1 and vascular plants (VP) from encroaching into peat moss (Sphagnum ssp.) dominated plant communities 2 .However, due to changing climate and artificial drainage for land use the pristine hydrological regime is often disturbed.For example, lower water levels foster aeration of the peat and hamper capillary water supply for Sphagnum mosses, both supporting an encroachment of graminoid and tree species 3 .This phenomenon can be observed in both natural and restored bog peatlands 4,5 .
The presence of VP can have contrasting effects for peatland water and greenhouse gas (GHG) cycling and can even improve conditions for their own expansion.Abundant root penetration of the peat could increase transpiration rate, leading to soil aeration and consequently to a further decrease of the groundwater level, creating a positive feedback loop for establishment of VP 6 .However, it has been shown in other studies, that birch encroachment does not necessarily have an impact on water table 7 , while increased graminoid biomass even diminished evaporation losses by attenuating wind speed at the moss surface 8 and thus protecting the moss from desiccation 9 .Further, there's evidence, that rather temperature increment seems to be the main driver of graminoid and tree dispersion compared to water regime 6,10 .Some findings suggest, that Sphagnum mosses desiccate permanently after extreme heat waves 11 while others indicate recovered growth and a decrease of tree cover if drought events are followed by a normal water regime 6,12 .Thus, although encroachment is observed also in natural raised bogs, restored bogs are probably more vulnerable due to disturbed hydrological regime and potentially insufficient rewetting measures and consequently, lower ability for recovering.
In a modelling study, Heijmans et al. 10 reported a negative relationship between VP expansion and peatland carbon accumulation.The presence of VP seems to prime microbial decomposition of organic matter, resulting in increased losses of carbon due to respiration or in dissolved form due to destabilization of organic matter 13,14 .However, short-term observations show that removal of VP resulted in decreased net CO 2 uptake by 50% during vegetation period in two investigated sites at different elevations 13 .Many research papers focus on how GHG emissions respond to changing water tables 15,16 or how vegetation composition is affecting biodiversity 17,18 .Due to the large amount of stored carbon in peatlands 1 , even small changes in climatic or hydrological conditions or vegetation composition can have significant impacts on the global carbon balance.When loss of peat carbon due to encroachment exceeds the enhanced sequestration by VP photosynthesis, the bog might turn from a sink into a source of carbon and thus current efforts of peatland protection and restoration for emission mitigation will be diminished.As shown above, implications of VP encroachment on water and carbon balances in bogs require more research, especially regarding predicted future climate conditions.There's also a need to understand physiological and physical key parameters of raised bogs responsible for alterations of the carbon cycle following restoration measures.
Conducting empirical studies to investigate these key parameters impacting carbon cycling would require years or even decades of measurements.Utilizing Soil-Vegetation-Atmosphere Transfer (SVAT) models is an effective tool to get insights into relevant physical and ecological processes and infer crucial parameters from flux simulations.A review of process-based models used for peatland carbon and water balance simulations 19 showed strong focus on northern peatlands and less on temperate ones.Of the 45 investigated models in that study, only three models specifically address eco-hydrology.Further, only 18 of them consider peatland specific properties such as microtopography, shallow water tables or peat accumulation and there's a large variation regarding temporal and spatial scales.Consequently, using an existing process-based peatland model for simulating CO 2 fluxes in temperate bogs that (i) represents eco-hydrological processes (ii) at local scale in temperate climate and (iii) in high temporal resolution, will most likely require some amount of structural model adaptation.
Here, we applied the multi-layer SVAT model pyAPES 20,21 , developed and tested under boreal conditions, to a temperate raised bog with shallow water table, peat properties and vegetation composition.We calibrated the model using soil temperature (T s ), water table dept (WTD) and CO 2 flux data 22 .pyAPES consists of sub-models for canopy and soil processes.Within the canopy sub-model, a separate moss layer is described that can be parametrized with specific traits measured in-situ or obtained from literature.This includes physical properties such as moss height (h m ) and bulk density (ρ m ), water retention characteristics, and photosynthetic traits 23,24 .Further, the model provides the opportunity to include vascular plants such as graminoids and tree species, characterised by a multi layered leaf area density profile and physical and physiological traits.As the model includes feedbacks between the multi-layer vegetation, microclimate and soil hydrology, pyAPES can partition fluxes originating from soil, moss and VP under changing environmental conditions.
The overall objective of this study is to identify the main eco-physiological traits of the two functional groups (Sphagnum moss and VP grass) and peat properties that control ecosystem level CO 2 uptake (gross primary production, GPP) and emissions (ecosystem respiration, R eco ), as well as their dynamics caused by environmental conditions.To achieve this goal, we address the following research questions: I. Which photosynthetic and respiratory parameters (traits) are most crucial for modelling CO 2 fluxes in a temperate raised bog?II.Is the observed relative importance of parameters/processes affected by seasonal variability in the environmental conditions (especially T s and WTD)?III.Does simulated VP encroachment (increasing leaf area index and decreasing Sphagnum dry mass) alter the importance of functional groups for CO 2 balances?
The results will contribute to better understand the interactions between plants and soil in raised bogs and how these interactions regulate the release and accumulation of carbon.Further, our investigation will give insights about pyAPES' performance in modelling CO 2 exchange and hydrology of raised bogs.

Sensitivity of annual fluxes
We found only minor changes in parameter ranking by extending the boundaries of Morris sensitivity analysis (MSA), either regarding GPP or R eco .However, extending boundaries for MSA enlarged the amplitude of all parameters.Figure 1 shows absolute mean (μ*) and standard deviation (σ) of EE.While dominating parameters that affected GPP annual sums were moss related (Fig. 1a), the highest impact on R eco annual sums was forced from soil parameters (Fig. 1b).
In case of the standardized ranking (± 30% boundaries), overall seven of total 22 parameters were identified to be important (μ* > 5 mol m −2 year −1 ) in affecting GPP annual sums.Parameters that are included in photosynthesis processes directly (Farquhar model in pyAPES) play a major role compared to moisture, temperature or radiation related traits.Maximum moss carboxylation velocity (V cmax,moss ) was clearly found to be the most important parameter, followed by its ratio to maximum transport rate of electrons (J/V c,moss ) and desiccation induced reduction in photosynthetic capacity (Δ desic ).The last four parameters that affected GPP most shared the same rank as their confidence intervals overlapped.No parameters of VP appeared significant for GPP of the natural bog, likely due to the low VP leaf-area.The highest rank of this plant functional group was eight, shared by V cmax,vas and J/V c,vas .Annual balances of R eco were most sensitive to peat respiration base rate and temperature sensitivity (r 10 and Q 10 ), followed by V cmax,moss and its ratio to dark respiration rate (R d /V c,moss ) and, finally, temperature dependency of dark respiration (tR d,moss ).Similar to GPP, no parameters of VP were considered important for R eco balances.
By expanding parameter boundaries of MSA from the previously used ± 30% to boundaries based in a literature review, the amplitude of all parameters increased.In case of GPP, two more parameters became important: the air-to-chloroplast conductance for CO 2 when external water has evaporated (g max ) and the temperature sensitivity of electron-transport limited photosynthesis rate (tJ moss ).On the other hand, the temperature sensitivity of RuBisCO limited photosynthesis rate (tV c,moss ) was no longer part of the group of most important parameters.Regarding R eco , only tR d was no longer found to be important, while Δ desic and VP parameter g 1 (stomatal slope) became important here.Overall, all parameters showing a high μ* (V cmax,moss and J/V c,moss for GPP; r 10 , Q 10 , V cmax,moss and R d /V c,moss for R eco ) did not shift more than one rank.Increased μ* and σ values were found, as well as wider confidence intervals (CI) and coefficients of variation (cv, σ/μ*) for all shown parameters (Fig. 1) compared to standardized boundaries for both, GPP and R eco fluxes.

Dependence of flux sensitivity on environmental conditions
With regard to research question II, we found an effect of WTD on both, GPP and R eco flux sensitivity for nearly all parameters (Fig. 2), and the parameter rankings for GPP were the most sensitive to changing environmental conditions.
GPP was much more susceptible to rate of decrease of photosynthetic capacity in dry moss (Δ desic ) at deeper WTD lower than -5 cm (Fig. 2a).In contrast, maximum symplast water content (θ symp ) was more relevant for GPP at water levels close to surface and showed decreasing importance with both declining and increasing WTD, but with a stronger intensity in the direction of deeper WTD.Lower sensitivity values for V cmax,moss only appeared, if water level and temperature were low at the same time.Increased sensitivity to J/V c,moss ratio was found with rising WTD compared to deeper WTD, given similar temperature levels.Temperature optimum of V cmax,moss (tV c,moss ) showed high and low sensitivities mixed over the whole range of WTD.However, at T s above 25 °C only very high ranks occurred.For R eco we found a strongly declining influence of r 10 and Q 10 with water table above soil surface (Fig. 2b).V cmax,moss and R d /V c,moss respectively showed slightly increasing ranks with shallower WTD while there was no clear pattern for the other parameters.
Temperature also affected flux sensitivities.Regarding tV c,moss there was a strong drop in its relevance for GPP with decreasing temperatures, but only when WTD was close to surface (i.e.no water limitations).J/V c,moss was declining more gradually with temperature for the whole range of WTD.Despite a small temperature dependency for V cmax,moss at deeper WTD, other important parameters showed no clear pattern.Most of the important parameters for R eco showed a relationship to T s .While the effect was minor, stronger impact of T s could be shown for tR d,moss with lower ranks by increasing temperatures.The only parameter that increases ranking with increasing T s , given lower WTD, was Q 10 .
To quantify these relationships, linear regression was used.For both GPP and R eco , the results showed strong impact of WTD and T s on average elementary effects μ* for all shown parameters (Table 1), illustrating clearly changing absolute sensitivities (μ*) with varying water table and soil temperature.Except tR d,moss , all regression models showed a better fit (R 2 ) to R eco values than to GPP, confirming less variation in ranking compared to ) at standardized (± 30% variation; blue dots) and literature boundaries (orange triangles).Inset plots show ranking of variables for both boundary conditions and thus display a change in ranking (ranking 1 = most important).Dashed/dotted grey lines representing isolines of the ratio of standard deviation (σ) and absolute mean (μ*) of the elementary effects.At a high ratio, the respective parameter has likely a non-linear effect and/or is involved in interactions with other parameters.
GPP.Nearly all parameters responded with increasing μ* to increasing WTD and T s , respectively.Only Δ desic was found to affect GPP fluxes less with higher water table.Further, the decreased impact of tR d,moss with rising soil temperature was confirmed.In all cases varying WTD by one unit affected fluxes much more than varying T s .Thus, in average CO 2 fluxes became more sensitive to all presented parameters (Fig. 2, Table 1) with increasing WTD and T s, except Δ desic and tR d .Relative standard errors were low for all coefficients, resulting in a narrow CI 95 and thus a robust estimate.

Impact of vegetation composition on flux sensitivities
The simulated shift in vegetation composition (research question III) had a clear impact on parameter rankings, more distinct for annual GPP than for R eco (Fig. 3).However, the V cmax,moss , r 10 and Q 10 remained the most important irrespective of decreasing m dry nor increasing vascular LAI.
The maximum carboxylation rate of vascular plants (V cmax,vas ) became increasingly more important for GPP with decreasing m dry , reaching even rank two when vascular LAI was at its highest.Same pattern, albeit at lower ranking, was observed for several other VP parameters (J/V c,vas , quantum yield α vas , stomatal conductance parameters g 1 and g 0 ).In contrast, moss parameters affecting GPP increased their importance with enhanced  www.nature.com/scientificreports/m dry and lower vascular LAI (R d /V c,moss , α moss , tV c,moss , tJ moss , tR d,moss ).In case of α and tV c , this led to a switch in the order of importance between moss and VP.
Ranking of parameters with respect to annual R eco was again more robust than that of GPP.Both soil parameters (r 10 and Q 10 ) were the dominating parameters followed by V cmax,moss and R d /V c,moss , both alternating at rank three and four.All of them were unaffected by changed m dry or vascular LAI.However, there were also some parameters that changed importance with simulated encroachment.The role of Δ desic and θ desic both decreased with VP encroachment, whereas the importance of V cmax,vas and R d /V c,vas increased.Of the VP parameter V cmax,vas was ranked highest with a value of seven.

Discussion
The calibration of the model to represent temperate raised bogs was conducted using observations from the Meerkolk site.Since pyAPES was originally developed for canopy covered boreal forest and peatland ecosystems 20,21 , adapting the model to temperate bogs required parameter calibration, mainly concerning peat hydrological properties.The most challenging part was the choice of shape parameters (α and n) of the unimodal van Genuchten-Mualem (vGM) water-retention model.McCarter and Price 25 showed that peat air-entry potential (α) can be highly variable in the uppermost peat layers consisting on poorly decomposed dead peat forming Sphagnum, ranging from 0.4 to 2.6.Values of n used in that study 25 were more robust (ranging between 1.2 and 1.5) for Sphagnum peat of different bulk densities.Price et al. 26 presented much higher values for weak decomposed Sphagnum mosses, which is used in the present study.High α and n values indicate high proportion of macropores and consequently an early entrance of air, well reflecting field observations, as the first soil horizon in Meerkolk was characterized as Von Post degree 1 (very weakly decomposed) and bulk density of living moss was comparably low.Further, higher values of α lead to higher K unsat , which may partly compensate the missing shrinkage of the modelled peat profile and thus better reproduce observed WTD.The applied combination of high α (0.8 27 ) and n (2.5 26 ) describing water retention of the first soil horizon up to − 18 cm in the present study was thus probably necessary to compensate a distinct swelling/shrinkage behaviour as discussed above.Using bimodal vGM curves 28 to describe dual porosity might thus be beneficial to better capture changes in porosity due to these processes in organic soils 31 .
The current parametrization of moss pF (α = 0.4; n = 2.8) is similar to those reported in Price et al. 26 (α = 0.5; n = 2.5) for living/undecomposed Sphagnum.Given the maximum gravimetric water content of moss layer w max that can be hold against gravity, saturated (θ sat ) and residual (θ res ) volumetric water content of moss are computed by using ρ m .ρ m varies strongly between Sphagnum species and even at the same plot between two years of measurement for the same species.For Sphagnum papillosum (dominating Meerkolk species), Bengtsson et al. 29 reported ρ m comparable to measurements at Meerkolk site (8.3 ± 0.97 kg m −3 ), but about 50% higher one year earlier (13.4 ± 1.63 kg m −3 ).Considering all reported species, ρ m in their study 29 ranged from 3 to 20 kg m −3 .Applying Eqs. ( 6) and ( 7) using this low ρ m , the resulting water retention curve was spanning a very narrow range of θ.This was consequently causing an implausible breakdown of capillary rise in summer 2018, since www.nature.com/scientificreports/unsaturated conductivity K unsat decreased too rapidly with θ, although water table was still high (about -10 cm).Since shrinkage due to desiccation is not implemented, it was necessary to constrain the minimum K unsat of moss to ~ 0.004 mm h −1 to enable capillary connection between the peat and the living moss layer throughout the modelling period.Overall, we highly recommend a parametrisation of both, ρ m and w max , by measured values, and if possible, measured water retention and conductivity characteristics.
Modelling soil temperature and ground water level was particularly challenging in the hot and dry summer conditions during European 2018 heatwave, while the WTD and soil temperature profile was well simulated during less extreme conditions in 2017 (Fig. 5).This was most likely related to shrinking and swelling of the peat body due to desiccation and rehydration that is currently not implemented in pyAPES.Nijp et al. 30 found, that in pristine bogs a change in absolute ground water table of -15 cm can change peat elevation up to -10 cm.Howie and Hebda 31 present similar or even higher values for very wet conditions and inundation, comparable to Meerkolk site with shallow WTD.Further, shrinkage reduces the pore size, which tends to lead to higher water retention and an increase in unsaturated conductivity (K unsat ) 32 , enhancing water supply from deeper soil layers.Including a dynamic soil profile might be more important in temperate compared to boreal latitudes fostering a more accurate description of the moss layer hydrology, even under increasingly extreme climate conditions.
Measurements of CO 2 fluxes of Oestmann et al. 22 were conducted from one hour before sunrise until late afternoon.It is common practice to measure GHG fluxes from close before sunrise until afternoon [33][34][35] , as maximum PAR and soil temperatures are already reached at this time of the day.Measurements over the whole day and even night time would potentially increase quality of the manually conducted model calibration.However, main focus of this paper are outcomes of the subsequent model sensitivity analysis, which does not take into account the observed fluxes explicitly.Thus, the impact of additional flux observations on results on model sensitivity is considered marginal.
It was necessary to scale Farquhar parameters with decreasing m dry for the encroachment simulation, since the model defines these parameters of the moss layer per unit ground area.Despite Scartazza et al. 36 , who found linear relationships of V cmax and J max with leaf mass of European beech (Fagus sylvatica) trees, evidence for the nature of this relationship to moss dry mass is missing.Thus, linear scaling was used.
Annual bog GPP was most sensitive to maximum carboxylation rate V cmax,moss , followed with greater distance by J/V c,moss , Δ desic , tV c,moss , θ symp , α moss and θ desic .Thus, under the studied environmental conditions, the bog GPP variability was primarily controlled by moss photosynthetic traits, water holding capacity, as well as the shape of photosynthetic moisture response.The relative importance (ranking) of V cmax,moss was barely affected, neither due to extended boundaries, WTD, T s nor increase of VP LAI.Only temperatures below 20 °C and WTD < − 7.5 cm were reducing its impact.This relation to WTD and T s could be observed for all important parameters, and is probably the joint effect of desiccation reduction and temperature dependant increase of photosynthesis.J max,moss (represented by its ratio J/V c ) is also determining photosynthesis rate of moss directly.Lower ranks of J max,moss were found, similar to V cmax,moss , only at low water tables and temperatures.Sensitivity to J max is thus also declining as maximum electron transport rate is never reached with deficit of both water and light.High temperature, on the other hand, is often associated with high irradiance which is increasing J and attenuating this limitation.Consequently, very low ranks of V cmax and J max , respectively could mainly be observed when both, WTD and T s were low.Albeit ranking of V cmax,moss did not change under simulated encroachment, i.e. increasing VP LAI, shading of VP started to decrease absolute sensitivity (μ*) of moss V cmax and J max (data not shown).This is not reflected in relative importance (ranking) since μ* of moss traits was still higher than those of VP.
Simulated VP encroachment (decreased moss dry mass and increased vascular LAI) both decreased μ* of V cmax,moss and J max,moss .While V cmax,moss still remained the most important parameter for GPP sums, J max,moss became less important than V cmax,vas at low m dry and high LAI.This is likely caused by shading (and thus decreasing electron transport due to decreased radiation) and indicates the onset of VP dominance.Measurements of grass LAI in restored, VP encroached bogs reached values up to 5.5 m 2 m −234 .This illustrates that LAI can be much higher than the maximum of 1.5 m 2 m −2 used here, which consequently also might reduce moss coverage.The results suggest an increasing importance of VP photosynthesis in encroached bog ecosystems while the share of Sphagnum mosses will decrease.Although mosses still remained as the main contributor of GPP (at least at VP LAI up to 1.5 m 2 m −2 ), accurately describing VP photosynthetic traits becomes more important at encroached or drained sites.
Effects of moss desiccation on GPP are prominent as all three moisture related parameters of moss photosynthesis are shown to be important, although WTD was relatively close to surface, even in summer 2018.These moisture related parameters will likely become even more important at sites with deeper WTD, like insufficiently restored sites or under strong VP encroachment.Especially Δ desic showed very high importance up to comparably shallow WTD of ~ − 0.07 m below soil surface and was the third most important parameter for annual GPP.Above this threshold, its impact declined very rapidly as the capillary rise from the water table was sufficiently strong to maintain Sphagnum moisture above θ desic , i.e. at non-limiting levels.The strength of decline in photosynthetic capacity (Δ desic ) seems to be of higher importance than the critical moisture value (θ desic ).This is emphasizing the strategy of Sphagnum to avoid rather than to tolerate desiccation, as several studies already suggested 37 .The decreasing importance of both, Δ desic and θ desic at higher LAI (even more distinct for R eco ) supports the finding of Heijmans 8 and Nichols and Brown 9 that grass canopy protects moss from wind and desiccation (reduced evaporation rate).Further, the role of these two parameters declined with higher m dry as the increased water holding capacity reduces the amplitude of moss moisture variations.The third parameter related to hydraulic traits (θ symp ) affects the shape of the water content response of air-chloroplast CO 2 conductance.It showed a comparable strong relationship between μ* and WTD (Table 1).At low water levels (and consequently low moss moisture) CO 2 diffusion is not inhibited by water films of external water.At high moss moisture θ symp becomes more influential, as it describes the moisture content above which the CO 2 conductance starts to decrease.Morris rank of θ symp , however, did not change linearly but showed an optimum instead, www.nature.com/scientificreports/i.e. decreased rank above and below water levels close to soil surface.Consequently, other parameters become more crucial here, e.g.maximum CO 2 conductance g max (related to efficiency of air-chloroplast CO 2 transport in non-limiting moisture) and quantum yield α (not presented due to μ* < 5 mol m −2 ).Overall, even under shallow WTD and near-natural Sphagnum cover as found at our test-site, moss desiccation behaviour can strongly impact ecosystem GPP and great care needs to be taken parameterizing its model representation, especially for future climate scenarios characterized by increasing magnitude and frequency of dry periods.An effect of shading from vascular plants is visible regarding the importance of moss photosynthetic temperature response (tV c,moss and tJ moss ).Due to shading (i.e.increasing overlying vascular LAI), moss temperature becomes less variable compared to direct radiation exposure.Thus, variability of photosynthetic rates will become less dependent on temperature, and importance of temperature response parameters will decrease (Fig. 3).As moss m dry per unit ground area is expected to decrease due to shading 38 under VP encroachment (as in our synthetic experiment), the moss will hold less moisture, leading to reduced heat capacity and resistance to temperature variations (i.e.decreased cooling effect).The risks for extreme moss temperatures and desiccation seem, however, to decrease as the sheltering from direct sunlight and wind exposure reduces evaporation rates.
Modelling annual R eco was strongly dependent on the soil parameters base respiration rate (r 10 ) and temperature sensitivity (Q 10 ).The behaviour of these parameters regarding environmental conditions was similar, and decreased in importance only at shallow water tables when the peat profile was saturated.As WTD is comparably close to surface in all cases, a reduction in microbial activity due to drought was not occurring, and sensitivity of R eco to these parameters was only marginally affected by environmental conditions.Ranking of r 10 and Q 10 was not altered by composition of VP, indicating a high contribution of soil respiration to total R eco .Since dark respiration rate of both moss and VP are defined proportional to V cmax , the R eco was sensitive to both parameters, regardless of the selection of boundaries.Thus, in contrast to GPP, moss respiration has a small and robust effect to R eco .Although there was a pattern of increasing importance of some VP parameters (e.g.V cmax and R d /V c ), they did not reach higher ranks as observed for GPP.
Overall, we demonstrated that sensitivity of CO 2 fluxes (GPP and R eco ) depends more on intra-annual changes of climate conditions compared to bog's plant traits or shifts in vegetation composition, i.e. there was a higher fluctuation of ranking for each parameter throughout the year than changed sensitivities of annual balances caused by wider boundaries or encroachment.This might be the result of the extreme conditions in 2018.However, despite this high variation throughout the year, parameter rankings were robust against extension of MSA boundaries.This suggests the transferability of the model's parameter importance to other temperate raised bog sites even at extreme years.We further illustrated initial dominance in photosynthesis contribution of VP traits (especially carboxylation rate) during simulated encroachment, shown by decreased rankings of moss parameters and at the same time increasing importance of VP parameters.In contrast, respiration is mainly driven by peat properties and water table regardless of vegetation composition.However, the current approach of computing R soil is not distinguishing between heterotrophic and autotrophic soil respiration.Since VP roots were present in saturated conditions, respiration fostered by aerenchyma would be valuable to include in the model.Especially at increasing LAI (and consequently increasing RAI) contribution of VP to R eco might be underestimated in the present results.
We conclude that pyAPES is capable for simulating CO 2 fluxes (Fig. 6), soil temperature and water table (Fig. 5) of temperate raised bogs.Under (near-) natural conditions, carboxylation rate of moss had dominating impact on bog GPP while peat respiration was the main factor affecting R eco .The impact on almost all plant traits was found to be highly non-linear, i.e. very high and low parameter values led to major changes in annual CO 2 fluxes.Thus, the extent of necessary parameter calibration highly depends on site conditions, but most important ones are robust in ranking.However, sensitivity of bog GPP and R eco to plant traits can change substantially depending on climate conditions.Especially WTD (and consequently traits associated with moss moisture) showed highly fluctuating ranks and were important for annual sums even though they had very low ranks during most of the modelling period.Therefore, if possible, we recommend measurements of moss moisture as well as hydrological properties of mosses as extensive as possible due to its crucial impact on bog photosynthesis.Consequently, it is necessary to consider moisture related traits even at comparable shallow WTD in raised bogs.During encroachment, VP related parameters became increasingly important.Most important upcoming traits were carboxylation and electron transport rate.This gives evidence that even at only slightly encroached bog sites photosynthetic traits of vascular start to become relevant for CO 2 balance.Here, moss V cmax still remained the most important trait but decreased in absolute importance (μ*) inversely to VP.We were able to delineate the most relevant plant physiological traits determining ecosystem carbon fluxes dependent on climate conditions and vegetation composition, which is useful for future modelling CO 2 fluxes from raised bog sites.

Site description
The model parameterization and sensitivity analysis were conducted using 'Meerkolk' bog, a nature conservation area located in Emsland, Northern Germany (52°38′N; 07°08′E) as a test area.With its peat thickness (3.5 m), vegetation cover (100% total coverage; 100% Sphagnum; 15% VP), water table depth WTD (− 0.02 m annual average) and peat properties, it owns characteristics of near-natural raised bogs in temperate climate conditions.The underlying peat was classified as Ombric Fibric Histosol 39 with an increasing degree of peat decomposition with depth from Von Post degrees 40 of H1 to H3 within the upper three horizons (53 cm depth) and H6 to H9 below.Mean annual precipitation is 791 mm and mean annual temperature is 9.8 °C according to a weather station of the German Weather service in Lingen (nearby Meerkolk site), documented from 1971 to 2000.During the modelling period (2017-07-01 to 2019-01-01) the sum of precipitation, aggregated to annual scale was 719 mm year −1 and average temperature was 11.1 °C.Monthly average temperatures and sums of precipitation are shown in Fig. 4. Most abundant VP species were Rhynchospora alba, Erica tetralix, Vaccinium oxycoccos and Eriophorum angustifolium.For a more detailed site description see Oestmann et al. 22 .
VP one-sided leaf area index (LAI m 2 m −2 ) was determined at the end of summer 2018 using LI-COR LAI-2200C Plant Canopy Analyzer.Measurements were replicated at three plots to capture spatial heterogeneity and conducted inside the frames of the flux measurement chambers (75 × 75 cm).At each plot eight above and eight below measurements were conducted, using a 45° angle view cap with view direction to the centre, as suggested in the LI-COR LAI 2200c manual for point measurements.Results were computed with the LI-COR 'FV2200' software (version 2.2.1).All measurements were averaged to yield a representative mean of annual maximum LAI of 0.47 m 2 m −2 .
Dominating bryophyte species was Sphagnum papillosum.Moss dry biomass m dry [kg m −2 ] was determined using a grid of 5 × 5 cm squares covering the whole plot resulting in overall 225 sections.Number of moss capitula were counted in every second section and scaled to square meter.Ten moss shoots (only green, i.e. photosynthetically active parts) per plot were sampled and length was measured, resulting in an overall mean of 3.5 cm.We assume that moss length is equal to moss height h m .Subsequently, samples were oven dried at 60 °C for at least 48 h to determine dry biomass m dry [kg m −2 ].Moss dry mass was then scaled-up using capitula counts of all sections: where m s [kg shoot −1 ] is the average dry mass of the sample, N represents the quantity of moss capitula, a defines the section area [m −2 ] and i and j are indices for plots and sections, respectively.m dry is the average of the three plots (0.27 kg m −2 ).Living moss bulk density (ρ m ) is the result of m dry divided by moss height (h m ) and equals 7.48 kg m −3 .
Further, we measured maximum water content of moss layer w max that can be retained against gravity.Ten sample cylinders with 3.5 cm height were cut from the moss carpet and saturated with water.Subsequently, samples were drained on suction plates applying a suction of 0 cm at the bottom of the cylinders until weight constancy.Mass was measured for drained and oven dried (105 °C) samples and w max was determined from the mean difference of both.h m , ρ m and w max were used for model parametrization.
All measurements are summarized in Table 2.

Calibration data
For Meerkolk site, continuous measurements of soil temperature (T s ) and water table depth (WTD) at 30 min resolution, and campaign-wise data on gross primary production (GPP) and ecosystem respiration (R eco ) were available.In the following, negative carbon fluxes indicate uptake, while positives indicate release of carbon.T s at 2 cm soil depth were measured by a weather station at a site close to Meerkolk with similar site conditions.WTD [m] was determined using groundwater pressure transducers; negative signs indicate WTD is below soil surface.In total 18 measurement campaigns were conducted to measure GPP and R eco at each of the three replication plots.While NEE was determined with transparent chambers, measurements of R eco were conducted using   opaque chambers.The measurements covering the whole amplitude of photosynthetic active radiation (PAR) and T s , starting an hour before sunrise (lowest values) and reaching maximum values in the afternoon.GPP was calculated by subtracting R eco measurements from the nearest in time NEE measurement.The measurements and data-processing is presented in detail in Oestmann et al. 22 .

Implementation of photosynthesis in pyAPES
The pyAPES is a multi-layer, multi-species soil-vegetation-atmosphere transfer model, which models ecosystem GPP, R eco , evapotranspiration, sensible heat and momentum fluxes using 1-D representation of the canopy-mosssoil continuum.Both moss and VP photosynthesis are implemented using Farquhar-von Caemmerer-Berry (FvCB) biochemical model 23 , the standard approach for modelling CO 2 assimilation in vascular plants 41 and mosses [42][43][44] .The solution of leaf (or moss layer) photosynthesis is coupled with solution of energy balance (leaf/ moss temperature and (stomatal) conductance), using iterative solutions between the leaf/moss processes and the ambient microclimatic gradients (CO 2 , air humidity and temperature) and short and long-wave radiation regime.For complete model description see Launiainen et al. 20 and for applications in boreal sites 21,45,46 .The model assumes, that net assimilation rate (A n ) is either limited by Rubisco (A c rate) or regeneration of Ribulose-1,5-bisphosphate (A j rate): where Rd is the dark respiration rate.For the subsequent sensitivity analysis, we used maximum carboxylation velocity (V cmax ), defining A c , as well as maximum transport rate of electrons (J max ) and quantum yield (α) which determine A j .The V cmax , J max and R d temperature response follows Medlyn et al. 24 , and we included optimum temperature of V cmax and J max and activation energy of R d into the sensitivity analysis.
Optimal stomatal conductance of VP (g s * ) is described in pyAPES following the Unified Stomatal Optimization (USO) approach of Medlyn et al. 47 : where g 0 is residual conductance, g 1 stomatal slope, C a the ambient CO 2 concentration, A the net assimilation rate [Eq.( 2)] and VPD the water vapour deficit (VPD).The g 0 and g 1 were also included in the MSA.
For mosses, CO 2 conductance decreases exponentially with water content (θ), based on an empirical function fitted to data of Williams and Flanagan 44 for Sphagnum and Pleurozium: With g max describing CO 2 conductance at maximum symplast water content (θ symp ), i.e. when all external water is evaporated and a 0 and a 1 are empirical parameters of the fitted function.g max is thus actually a shape parameter instead of maximum conductance.
The effect of desiccation on photosynthetic capacity (cap photo ) at low moss moisture content (a multiplier for moss V cmax ) is described as Thus, desiccation leads to decrease of moss assimilation capacity (V cmax , J max ) and dark respiration rate (R d ) if moss water content (θ) is lower than a threshold (θ desic ), while the strength of the decay is defined by the slope parameter Δ desic .Both reduction and recovery of photosynthesis is assumed instantaneous and reversible in the model.Both Δ desic and θ desic were included in MSA.The moss moisture content is modelled at 30 min timestep, and accounts for rainfall interception, evaporation as well as capillary rise from underlying soil (affected by WTD and peat hydraulic properties).

Model parametrization and adaptation to site conditions
Measurements of T s , WTD and carbon dioxide fluxes, i.e.GPP, R eco and NEE were used to calibrate the model for temperate bog.Wherever possible we used measured values for parametrization.If no measurements were available, we conducted a literature review.In contrast to VP, literature values of Sphagnum moss hydraulic or photosynthetic properties are rather scarce.When no value could be found for a certain parameter, we used values from VP studies instead.If this was still not possible, we inspected respective functions and set parameter values to best fit measurements and plausibility of related model outputs.
In a first step we calibrated the heat flow of pyAPES' soil sub model using measured T s in 2 cm soil depth.No adaptation of parameters was necessary here, since observations were overall well reproduced by pyAPES with a RMSE of 2.5 °C.However, in some periods (mainly summer 2018), modelled T s is oscillating stronger than measurements.
As several studies have pointed out the importance of moss moisture for photosynthetic capacity 37,48 , we decided to calibrate WTD subsequently.It acts as a proxy for moss moisture, from which no data was available.Simulated moss moisture and capillary rise were observed carefully to show realistic dynamics.Calibrating the model to water table depth was challenging, and parameters of the water retention curve (WRC) 49 of both living moss and peat were most crucial for this part.We tested numerous sets of retention parameters for Sphagnum peat reported in McCarter and Price 25 , but none of these could reproduce WTD dynamics during the hot and (2) dry summer of 2018.Finally, we applied a mixture of Liu & Lennartz 27 for the air entry value α, and Price et al. 26 for the shape parameter n to the upper three soil horizons.Hooghoudt equation 50 describes lateral flow in the model.At Meerkolk, lateral drainage to the surrounding, lower agricultural areas plays a significant role, especially in times of water excess.Thus, amount and depth of lateral outflow was set to well represent winter water tables.Finally, the model was able to reasonably well reproduce the dynamics of changing water tables, even in summer 2018, with a RMSE of 2.8 cm (Fig. 5).
Water retention parameters of living moss micropores are, similar to the peat soil horizons, described by Van Genuchten Mualem equation, following the analogy of Voortman et al. 51 .The authors summarize several studies, showing high variations for α and n, as well as differences in their combination to define the shape of water retention curve.In the model the water content at saturation (θ sat ) and residual water content (θ res ) are defined by maximum gravimetric water (w max ) content of moss and maximum symplast water content (θ symp ), respectively, which is implemented in pyAPES as a ratio of w max .They are converted in the model to volumetric water contents using moss bulk density (ρ m ).
Using measured ρ m , moss WRC was very narrow and default α and n parameters resulted in a too strong decrease in WTD in dry periods.We increased the n value, which improved the modelled WTD but caused a break in capillary rise in summer 2018.By enhancing minimum hydraulic conductivity this issue could be solved.Roughness height of moss was further set to 1/10 of h m .
As VP are not abundant in Meerkolk, maximum vascular LAI is 0.47 m 2 m −2 on average.This value only represents a maximum LAI, which is not sufficient for model parametrization.To display seasonal dynamics, measurements from a comparable site 34 were used to calibrate development of modelled LAI by adjusting the base temperature for degree-day calculation for phenological development simulation (Fig. S1).
In the default settings of pyAPES, Farquhar parameters (V cmax , J max and R d ) for moss community are presented per unit ground area, and light attenuation within the single-layer moss is not accounted for.Since one of our research questions aims at changing vegetation composition (and thus changing m dry ), it was necessary to include a scaling of Farquhar parameters with moss dry mass.Currently, published relationships of V cmax etc. and moss dry mass are scarce, however at least J max seems to decrease almost linearly up to a depth of ca. 4 cm in a Pleurozium schreberi moss canopy 46 and non-linear functions with respective shape parameters are missing in the literature, we decided for a simple linear scaling with m dry .
Maximum velocity of carboxylation (V cmax ) for VP was set to average values of Ekberg et al. 53 for Eriophorum angustifolium at 20 °C.Maximum electron transport rate (J max ) and basic dark respiration rate (R d ) are set as ratios of 1.97 48 and 0.03, respectively of V cmax , as for mosses.
Heterotrophic and autotrophic soil respiration (R soil ) in pyAPES is computed using the approach of Pumpanen et al. 54 and moisture response of Skopp et al. 55 .By default, R soil is scaled by relative root area density in the respective soil layer.For bogs with shallow water table, we assume bulk peat respiration takes place only in the unsaturated zone.R soil is thus computed as layer-thickness weighted average from all soil layers above WTD, assuming no respiration within water saturated layers: where R and dz are the soil respiration and thickness of layer i, z the soil depth and the subscript n denotes the index of the corresponding layer of current WTD.
Using the parametrization of Table S3, the model slightly overestimated fluxes.Due to higher overall GPP bias (− 0.7 g CO 2 m −2 day −1 ) compared to R eco (0.5 g CO 2 m −2 day −1 ), simulated NEE is tending to show more CO 2 uptake (− 0.6 g CO 2 m −2 day −1 ) than observed (more details in Table S2).Both, bias and RMSE varied strongly between campaigns and had lowest errors in winter months (at low fluxes) and highest in summer (Fig. 6).All adapted parameters of pyAPES are presented in Table S3.

Morris sensitivity analysis
To find most relevant parameters of the pyAPES' carbon module for bog GPP and Reco, sensitivity analysis was conducted using the approach of Morris (MSA) 56 .The method is based on evaluation of elementary effects (EE), i.e. evaluating the alteration of output variables by adjusting only one factor at a time (OAT).These local OAT approaches are less expensive in computation time, as they need less model runs in contrast to global approaches.However, different from other OAT methods, Morris approach uses randomly sampled parameter combinations.Thus, the whole input space of parameters is covered and Morris method can be seen as computationally less expensive approximation of global sensitivity analysis.( 6) To compute an EE of one parameter, a base parameter set must be created.For this purpose, lower and upper boundaries must be defined for each of the k parameters.Within these boundaries, a parameter can take a value out of p levels in an even distance from each other.This results in a k times p dimensional grid, representing the input space and called Ω. From Ω a random value is drawn for each parameter X i to get the base parameter set.In a next step, each X i is changed by Δ, without returning to its initial base value.This procedure requires k + 1 parameter sets (base set plus k changes, one for each X i ) and is called a trajectory.One trajectory yields exactly one EE for each input parameter.Repeating this procedure r times, one receives a distribution of EE, F i, for each input parameter with computational costs of r(k + 1) model evaluations.Mean value μ, standard deviation σ can now be computed from F i .According to the suggestions of Campolongo et al. 57 , a distribution of absolute values from EE, G i (and whose mean value is called μ*), is used in this study, to avoid nullifying effects of opposite signs when the model is non-monotonic.The coefficient of variation (σ/μ*) can give evidence whether the effect of an input parameter is (i) linear and additive, (ii) non-linear or/and interacting or (iii) negligible.
In this study, parameters were allowed to vary at four levels within given boundaries.For each distribution G i , 1000 bootstrap runs were conducted to calculate a 95% confidence interval (CI) for μ*.Following the approach of Campolongo et al. 57 , 1000 trajectories were created and 50 optimal were chosen for the analysis.'Optimal' in this context means reaching the widest spread over the whole input space (Euclidean distance).As the number of trajectories for original Morris analysis is between 10 and 50 51 , we used 50 optimal trajectories for the enhanced method.To evaluate whether this number of optimal trajectories leads to reliable results we tested different seeds (i.e. starting points for generating pseudo-random numbers) and compared rankings of the respective analysis.Ranking of the parameters was the same for all seeds, considering their confidence intervals, when using 50 optimal trajectories.Final analysis was conducted without setting any seed.Parameter were ranked according to their amount of μ* without considering their confidence intervals.Parameters with equal μ* all got the lowest ranking of this group (e.g.[1, 2, 5, 5, 5, 6, …]).
Overall, 22 parameters were investigated (Table S3).Kattge and Knorr 58 reported a linear relationship between J max and V cmax at 25 °C, but found this ratio to be depended on growth temperature.Thus, we decided to define both, J max and R d as ratios of V cmax and varied V cmax and these ratios within MSA.
We first investigated sensitivity of annual balances while varying parameters within standardized boundaries as well as broader boundaries (i.e.literature observations).Further we analysed the impact of changing site conditions and vegetation composition to flux sensitivities.For all computations regarding Morris sensitivity analysis Pythons SALib package version 1.4.5 59 was used.

Sensitivity of annual balances
To identify the most crucial parameters for modelling of the overall bog CO 2 exchange, annual balances of GPP and R eco were used as output variables in the SA.Parameter boundaries were defined to be ± 30% of the parametrization used to reproduce best the observations of the study site.Using a uniform boundary definition among all parameters minimized impacts on SA results and thus enabled a standardized ranking.In a second step, boundaries were expanded to the broader range of measured parameter values found in a literature review, to investigate whether parameter ranking is likely to change when applying the model to other bog sites.Analogous to model parametrization we derived values from VP if literature did not provide moss parameters.
J max /V cmax ratios are either reported directly or calculated from both values.Same ranges for θ desic and θ symp were used, as most studies report optimal water contents for photosynthesis but not distinguish between desiccation and CO 2 diffusion.

Response of parameter sensitivities to seasonal site condition dynamics
To identify, whether environmental conditions such as WTD or T s affected carbon flux sensitivities, Morris analysis was conducted not only for annual sums, but also for each timestep in half-hourly resolution over the simulation period.Sensitivities of carbon fluxes with respect to concurrent simulated WTD and T s were evaluated to assess how the parameter ranking is affected by seasonal variability in the environmental conditions (research question II).
Thus, average elementary effects could be computed for each timestep.As WTD and T s values were possibly affected by changing parameters, WTD and T s timeseries were sampled from each individual model run within the SA.Following this approach for all trajectories, mean values of WTD and T s distributions for all timesteps were computed.We investigated the five most influential parameters from SA with standardized boundaries.Linear regression analysis was conducted using R software (RStudio version 1.

Response of flux sensitivities to changing composition of plant functional groups
In a case of VP encroachment, increase of VP LAI and decrease of moss biomass is expected due to shading 38 .To investigate the impact of vegetation composition on flux sensitivities and parameter importance, we generated synthetic parameter sets in the range of 0.5-1.5 m 2 m −2 LAI max and 0.15-0.35kg m −2 m dry .To set the upper and lower boundary for LAI max and m dry , respectively, we used measured values from another bog site in Northern Germany with strong encroaching vascular species.The vegetation composition of the Meerkolk experimental site well represents near-natural bog vegetation composition.Observed LAI max and m dry ranges were divided in five even spaced steps and Morris analysis was conducted once for each combination within this 5 × 5 matrix.

Figure 1 .
Figure 1.Average sensitivities of annual balances of (a) GPP and (b) R eco to important model parameters (μ* > 5 mol m-2 yr -1 ) at standardized (± 30% variation; blue dots) and literature boundaries (orange triangles).Inset plots show ranking of variables for both boundary conditions and thus display a change in ranking (ranking 1 = most important).Dashed/dotted grey lines representing isolines of the ratio of standard deviation (σ) and absolute mean (μ*) of the elementary effects.At a high ratio, the respective parameter has likely a non-linear effect and/or is involved in interactions with other parameters.

Figure 2 .
Figure 2. Impact of water table depth (WTD) and soil temperature (T s ) of the top soil layer (1 cm depth) to Morris ranking of the five most important model parameters (according to absolute mean without considering confidence intervals) of (a) GPP and (b) R eco .Negative signs of WTD indicate water table is below surface.

Figure 3 .
Figure 3. Impact of simulated change in vegetation composition (leaf area index of vascular plants, LAI and moss dry mass, m dry ) on annual sums of (a) GPP and (b) R eco while moss ground cover was still assumed to be nearly 100 %.Numbers and colors in panels indicate Morris ranking (1/light = most important to 22/dark = less important). https://doi.org/10.1038/s41598-024-61229-6

Figure 4 .
Figure 4. Climatic conditions at Meerkolk site during modelling period.Figure shows monthly means of air temperature (line) and monthly sums of precipitation (bars).Over the whole 1.5-year modelling period average temperature was 11.1 °C and sum of precipitation was 1079 mm (719 mm yr -1 ).

Figure 6 .
Figure 6.Measured (dots) and simulated (lines) fluxes of GPP (green), R eco (brown) and NEE (beige) for each measurement campaign during the modelling period.
3.959) to find significant relationships.A regression model of the form Y = a + b × WTD + c × T s was used.

Table 2 .
Measured plant properties at Meerkolk site.
Figure 5. Measured (green) and modelled (orange) soil temperatures at 2 cm depth (T s , top, RSME = 2.5 °C) and water table depth (WTD, bottom, RSME = 2.8 cm).A negative sign of WTD indicate water table is below surface.