Peatland-VU-NUCOM (PVN 1.0): using dynamic plant functional types to model peatland vegetation, CH 4 , and CO 2 emissions

. Despite covering only 3 % of the planet’s land surface, peatlands store 30 % of the planet’s terrestrial carbon. The net greenhouse gas (GHG) emissions from peatlands depend on many factors but primarily soil temperature, vegetation composition, water level and drainage


Introduction
Despite covering only 3 % of the planet's land surface, peatlands store 30 % (644 GtC) of the planet's terrestrial carbon (Yu et al., 2010).The present-day global radiative effect of peatlands on the climate is estimated to be between −0.2 and −0.5 W m −2 (i.e. a net cooling) (Frolking and Roulet, 2007) in comparison to a radiative forcing of +2.43 Wm −2 due to anthropogenic greenhouse gas (GHG) emissions since pre-industrial times (WGI, 2021).Future changes to the climate will impact the carbon sequestration capacity of peatlands; however, the net effect of climate change on peatlands is not yet understood (Loisel et al., 2021).Research indicates that some peatlands will form a positive feedback (Dorrepaal et al., 2009), whilst others will form a neutral (Saleska et al., 2002) or negative feedback to warming of the global climate system (Melillo et al., 2002;Lafleur et al., 2003), and the net effect of these complex responses is not yet known.The net warming effect of peatlands on the global climate system, and particularly the potential to both emit and draw-Published by Copernicus Publications on behalf of the European Geosciences Union.
down CO 2 and CH 4 , means that peatlands have a complex and multifaceted relationship with the global climate system.
The net GHG emissions from peatlands depends on many factors but primarily vegetation composition, land management, ground water level and drainage, and soil temperature (Dorrepaal et al., 2009;Tiemeyer et al., 2016).Rewetting drained peatlands is one strategy proposed to combat enhanced CO 2 emissions from peatlands, but this has been documented to both enhance and reduce GHG emissions (e.g.Günther et al., 2020;Boonman et al., 2022) with the majority of studies concluding that rewetting leads to enhanced CH 4 and net GHG emissions (Harpenslager et al., 2015;Knox et al., 2015).Field studies have shown that vegetation restoration in combination with rewetting may reduce GHG emissions (Graf and Rochefort, 2009;Mazzola et al., 2022).Vegetation impacts the net GHG emissions in peatlands by directly influencing the net primary production and organic matter available for decomposition and indirectly and by influencing the substrates available for microbial metabolisation in the soil column (Bansal et al., 2020;Bridgham et al., 2013).
While the effects of the groundwater table on peatland GHG emissions are extensively described (Evans et al., 2021), the impacts of plant type and plant community composition on GHG emissions are less understood (Malmer et al., 2003).Changes in vegetation composition have been observed in long-running water table manipulation experiments (Peltoniemi et al., 2009;Strack et al., 2006).Generally, sedges and mosses establish during wetter conditions, and shrubs and trees develop during drier conditions, with enhanced Sphagnum growth outcompeting shrubs during warming experiments (Dorrepaal et al., 2006).Belowground, changes in vegetation have been accompanied by changes in bacterial and fungal biomass (Jaatinen et al., 2008) and methanogenic and methanotrophic community diversity (Yrjälä et al., 2011;Lippmann et al., 2021).Changes to CO 2 (NPP) have been observed following changes in plant community composition, further impacting root exudation (Ballantyne et al., 2014).Root exudates are a diverse group of organic compounds secreted by plant roots into the nearby soil.The composition and quality of root exudates varies between plant types, influencing microbial community composition and function and CO 2 (Crow and Wieder, 2005) and CH 4 fluxes (Schipper and Reddy, 1996).Plant growth, root exudation, and decomposition of organic matter happen at rates that differ depending on plant type (Dorrepaal et al., 2007).Sphagnum is a primary contributor to the carbon sequestration in many peatlands and decomposes 3 times slower than most vascular plants (Graf and Rochefort, 2009).Spatial variation in the rate of vegetation growth and decomposition, particularly for bryophyte species, leads to the creation of microforms, such as hummocks, hollows, and lawns, which, in turn, impact the water level relative to the surface and spatially variable fluxes (Waddington and Roulet, 2000).Differences in vegetation composition within the same site and with the same water levels have been observed to lead to differences in CH 4 fluxes (Bubier, 2016;Jackowicz-Korczyński et al., 2010).To understand the role of vegetation emissions' feedbacks during peatland restoration efforts, vegetation must thus be treated as a dynamic interactive element of the peatland ecosystem.
Plants with common ecosystem functions or structures can be represented with common model algorithms or parameters in a vegetation model when grouped as plant functional types (PFTs) (Wullschleger et al., 2014).Plant functional types have been found to explain uncertainties in GHG emissions from wetlands in response to warming in a meta-analysis of wetlands exposed to warming (Bao et al., 2023).Dynamic (rather than static) PFTs simulate the inter-seasonal growing and dying of plants that, over a number of years, lead to vegetation succession and are critical to reliably assess the impacts of climate and environmental change on peatland ecosystems (Box et al., 2019).Shifts in community composition lead to feedbacks between species and other environmental parameters such as soil moisture, bulk density, soil organic matter (SOM) content, gas conduit function, rate of growth, rate of decomposition, microbial mineralisation, and aerobic decomposition (De Boeck et al., 2011).Dynamic plant representation is critical to reliably simulate vegetation-environmental feedbacks in models (Toet et al., 2006); therefore, the inclusion of dynamic vegetation classes is critical to reliably estimate C, CO 2 , and CH 4 emissions from peatlands during periods of environmental change (Li et al., 2016;Laine et al., 2022).
Many peatland carbon cycle models have been developed over the preceding decades.The Wetland and Wetland CH 4 Inter-comparison of Models Project (WETCHIMP) evaluated the ability of a variety of models to simulate largescale wetland characteristics and corresponding CH 4 emissions (Melton et al., 2013;Wania et al., 2013).Peatland modelling efforts have made significant advancements to simulate CH 4 fluxes by including CH 4 -specific processes such as CH 4 plant transport and ebullition.However, many models rely on CO 2 fluxes or surface water levels as indicators of CH 4 exchange (Metzger et al., 2015), restricting their capacity to assess feedbacks between environmental change and the peatland CH 4 cycle.There exist two pre-existing models that simulate dynamic vegetation, CO 2 , and CH 4 cycling in peatlands (i.e.PEATBOG (Wu et al., 2016) and LPJ-WHyMe (Wania et al., 2010)), thereby limiting the ability of the modelling community to assess model mechanistic processes.The functionality and scope of current models that simulate peatlands and include either dynamic or static vegetation are compared in Table S1 in the Supplement.The PEATBOG model simulates three PFTs -moss, shrubs, and graminoids -at the Mer Bleue bog site and represents a comprehensive array of peatland processes, including the nitrogen cycle and dissolved gases (carbon, CO 2 , and CH 4 ).LPJ-WHyMe, like its parent model, LPJ-WHy (Sitch et al., 2003;Gerten et al., 2004), includes permafrost Figure 1.Schematic of the movement of carbon in the model.Processes are delineated with rectangles, whereas carbon pools are delineated with curved edges.The pink outline represents non-moss pools and processes, the green outline represents pools and processes applicable only to moss PFTs, and the blue outline refers to pools and processes that are applicable to all plant types.In the background of this figure, the Horstermeer site is shown on the left, and the Ilperveld site is shown on the right.and peatlands, two peatland-specific PFTs (flood-tolerant C 3 graminoids and Sphagnum mosses), a new decomposition scheme when under inundation, and the addition of root exudates.LPJ-WHyMe particularly assesses the impacts of inundation on vegetation composition, net primary production, and the deceleration of decomposition under inundation.
To assess the impact of dynamic vegetation classes on subsequent GHG fluxes in peatlands, we present a new model, Peatland-VU-NUCOM v1.0 (PVN).PVN incorporates features of NUCOM-BOG, a bog ecosystem model (Heijmans and Berendse, 2008), into the Peatland-VU model framework, a process-based peatland model (van Huissteden et al., 2006).The NUCOM-BOG model simulates vegetation competition and C, nutrient, and water cycling in undisturbed bog ecosystems under changing climates using a soil profile divided by an acrotelm-catotelm boundary where plant growth and decomposition are partitioned between plant organs.The Peatland-VU model simulates the CH 4 and CO 2 cycle within a column of peat soil with varying water levels.The Peatland-VU model simulates CH 4 fluxes, gross primary productivity, and CO 2 cycle whilst assuming a constant plant layer, and it does not include a nitrogen cycle.We have developed a model that, with the appropriate site input data, can be used to simulate peatland sites with a wide variety of vegetation types and vegetation management practices.For this reason, there is no limit to the number of PFTs that can be included in a model simulation.The inclusion of dynamic vegetation classes into the PVN model provides a model that is capable of estimating the greenhouse gas balance in response to environmental changes (changes in temperature, radiation, precipitation or evapotranspiration, or water levels) and also different management efforts (changes in harvest regime or vegetation restoration) for peatland sites.The incorporation of features of NUCOM-BOG, a model simulating undisturbed systems, into PVN, a new model simulating disturbed and managed systems, requires that changing environmental conditions and changing management practices both lead to dynamic impacts on vegetation classes.Therefore, this model can serve wetland management by estimating changes in the greenhouse gas balance of peatland sites in response to management decisions whilst considering the effects of environmental change.
2 Materials and methods

The PVN model
The new PVN model describes the vegetation, CH 4 , and CO 2 dynamics of a column of an above-and belowground peatland ecosystem (Fig. 1).Carbon dioxide and CH 4 emissions enter the atmosphere by ebullition, transport through plants, diffusion through the soil, and respiration.The aboveground carbon pools are the aboveground living biomass, litter layer (non-moss PFTs only), shoots, and living-moss depth (moss PFTs only).The belowground carbon pools are the peat, lahttps://doi.org/10.5194/gmd-16-6773-2023 Geosci.Model Dev., 16, 6773-6804, 2023 bile organic matter, microbial biomass, litter and dead roots, and root exudates (Table S2).
Plant functional types (PFT) are the key element of NU-COM that is added to the Peatland-VU framework to create the PVN model.Any number of PFTs can be included in a model simulation.PFT attributes (parameters) describe plant physiology and bioclimatic limits.Bioclimatic limits are used by the photosynthesis function (Sect.2.1.1)and the potential growth function (Eq.13).Each PFT is defined as being either a moss or vascular plant type, which impacts the ability of plants to grow vertically or to develop roots.Each PFT is prescribed as having either evergreen or deciduous phenology.For deciduous vegetation, leaf senescence increases when the daily temperature falls below the PFT's minimum tolerated temperature, whereas, for evergreen vegetation, leaf senescence refers to the death of old leaves (Eq.9).Maximum leaf coverage is maintained as long as the daily water level and temperature are within the ideal range.The PFT parameters are defined in Table 1, and the references are listed in Table S3.In this section, the subscript p is used to show that the equation or variable is PFT specific, z is used to indicate that the equation or variable is soil layer specific, t is used to represent time, T represents temperature, and WL represents water level.The convention used in this paper is such that a positive flux represents the movement of gas from the ecosystem to the atmosphere.

Primary production
C3 photosynthesis, leaf respiration (RT), and net primary production (NPP) are calculated using a modified version of the primary production scheme introduced into the Peatland-VU model by Mi et al. (2014), modified from the BIOME3 equilibrium biosphere model Haxeltine and Prentice (1996).The BIOME3 model is based on the premise that GPP and leaf respiration increase with the activity of (Rubisco) photosynthetic enzymes in leaf chloroplasts.Photosynthesis is calculated using stomatal conductance and Rubisco activity of leaves.The net CO 2 fluxes (NEE) for each PFT are the sum of gross primary production (GPP [kg C m −2 d −1 ]) minus plant respiration, CO 2 produced by belowground aerobic SOM decomposition, and CO 2 oxidised from CH 4 (Rox).
In the above equation, JE [kg C m −2 d −1 ] describes the relationship of photosynthesis to photosynthetically active radiation (PAR), and JC [kg C m −2 d −1 ] describes the Rubiscolimited rate of photosynthesis.JE and JC are defined by Eqs.(S2) and (S10) in the Supplement, respectively (Haxeltine and Prentice, 1996;Mi et al., 2014).Interactions among leaf area development, photosynthetic activity, stomatal conductance, temperature, and water availability have been widely recognised (Baldocchi and Harley, 1995;Koebsch et al., 2020).Water stress has a significant impact on plant photosynthetic capacity (Keenan et al., 2010).Studies such as those of Ostle et al. (2009) and Puma et al. (2013) have considered these factors when simulating GPP by introducing water use efficiency terms.Model intercomparison efforts have found improved reproducibility of GPP estimates from models that account for the impacts of water stress on photosynthetic capacity (De La Motte et al., 2020;Churkina et al., 1999).GPP in the PVN model is modified by both a water stress factor (WSF [-], Eq.S1) and a temperature stress factor (φ T , Eq. S4, adapted from Mi et al. (2014)).
In the above equation, NPP represents the net primary productivity In the above equation, Rr [-] is the leaf respiration coefficient (Table 1), and VM [kg C m −2 d −1 ] represents the maximum daily rate of net photosynthesis (Eq.S11, Haxeltine and Prentice (1996); Mi et al. (2014)).The CO 2 flux from each soil layer (B CO 2 ) is calculated before integrating over all layers; it is then summed with CO 2 produced by decomposed litter (LLd [kg C m −2 d −1 ]), and NPP is subtracted.
In the above equation, NEE [kg C m −2 d −1 ] is the net ecosystem exchange, and ] is the CO 2 flux produced by belowground SOM decomposition (Eq.27).

Competition among PFTs
Biomass fraction (BF) is a representation of the ratio of PFT biomass to total biomass (Eq.5).The sum of all PFTs is constrained to a maximum BF of 1.0.All PFTs have a minimum BF of 0.1 and are able to further establish when the conditions become favourable, as adapted from the NUCOM-BOG model (Heijmans and Berendse, 2008).
In the above equation, CB [kg C m −2 ] represents aboveground living biomass (Eq.9).Each plant competes for light  Krinner et al. (2005), which relates vegetation biomass to height.This relationship, initially intended to be used for trees, has since been used to calculate the heights of natural and agricultural grasses in a dynamic global vegetation model (Krinner et al., 2005).Biomass and stem density have been found to respectively explain 98 % and 81 % of the height variance in 65 plots of 29 different species (Gorham, 1979) because most plants are understood to be constrained by self-thinning under crowding in natural stands or by a trade-off between height and foliage growth, reflecting a trade-off between structural and functional physiological development.
In the above equation, H refers to plant height [m]; BD represents biomass density [kg C m −3 ], k 2 [m]; and k 3 [-] are constants with values of 40 and 0.85, taken from Smith et al. (2001).FPAR [-] is the fraction of incoming PAR absorbed by vegetation (Eq.7) and is dependent on LAI and the amount of shading by taller plants.
In the above equation, LEC represents the light extinction coefficient parameter [-].LAI [m 2 m −2 ] is calculated as a function of living biomass and the specific leaf area (SLA In the above equation, CB [kg C m −2 ] represents aboveground living biomass (Eq.9), dependent on shoot growth and biomass senescence lost to the litter layer.
In the above equation, SM represents shoot mass [kg C m −2 d −1 ], calculated using Eq. ( 10), and BS In the above equation, RS [-] represents the ratio of shoot to root growth (Table 1).The allocation of root and shoot growth is a fixed fraction of NPP so that the fraction of shoot and root growth sums to 1.0.
In the above equation, WL refers to water level, min and max refer to the minimum and maximum water levels tolerated for growth, and minopt and maxopt refer to the minimum and maximum optimum water levels for growth.
In the above equation, T refers to daily temperature, min and max refer to the minimum and maximum tolerated temperatures for growth, and minopt and maxopt refer to minimum and maximum optimum temperatures for growth.

Belowground production
The root distribution and the root mass of vascular PFTs are mapped to the layout of the model's soil horizon representation (depth, density, and layer thickness).To account for differences in decomposition rates among roots and exudates, each PFT has designated SOM pools, which are partitioned between the soil layers.Root distribution and root mass decrease exponentially from the surface to the PFT maximum root depth (MRD in Table 1).In general, 30 %, 50 %, and 75 % of roots are observed in the top 10, 20, and 40 cm, respectively (Jackson et al., 1996).Root exudation plays an important role in the rhizosphere by promoting methanogenesis and soil carbon loss through CH 4 production.The production of new roots (Rd) is based on a PFT-prescribed shoot-toroot-growth ratio and NPP.Root exudates (RX, Eq. 19) are a fraction of the calculated belowground root production (Rd).
Exudates develop at a prescribed rate per PFT, dependent on root and shoot growth.Photosynthesis rates are enhanced during spring and summer and are accompanied by the highest levels of root and soil respiration (Högberg et al., 2001).
There is strong evidence to suggest that enhanced photosynthesis fuels exudate production, causing seasonal variation in exudation (Whipps, 1990;Saarnio et al., 2004).The root growth and die-off functions are adapted from van Huissteden et al. (2006).
In the above equation, 1− RS represents the fraction of growth that is root growth, and f (z, p) [m −1 ] represents the exponential root distribution from the surface to the maximum root depth (MRD in Table 1). 0 In the above equation, RM is the root mass In the above equation, DoY represents the day of the year; REX [-] represents the unitless root exudation factor; and f (KSP) [-] is a function depending on the PFT constant, KSP (Table 1), that can be used to determine when stronger exudation occurs during spring.
In the above equation, RSX represents the root senescence rate [d −1 ].

Litter layer production and decomposition
Vegetation composition change directly impacts litter input, which alters the quality and quantity of fresh SOM contributions (Malmer et al., 2005).Senescence of the aboveground living biomass is added to the litter layer for vascular PFTs https://doi.org/10.5194/gmd-16-6773-2023Geosci.Model Dev., 16, 6773-6804, 2023 (Eq.21).Senescence of moss PFTs contributes directly to the belowground SOM pools.Movement of surface litter to SOM pools is an important component of peatlands (Davidson and Janssens, 2006).Carbon dioxide produced from the decomposition of the litter layer and the different SOM pools are summed with NEE (Eq.4).
In the above equation, where  2006)).

Belowground SOM decomposition
Peatlands consists of organic compounds at different stages of decomposition.In the model, these belowground organic components are separated into five SOM pools (peat, humus, microbial biomass, litter and dead roots, root exudates) (Table S2).Each of the SOM pools lose and gain mass, whilst the number and the thicknesses of the soil layers remain constant throughout the model simulation.Biodegradation of SOM leads to the mineralisation of carbon that can be reincorporated into SOM and repeatedly recycled (Basile-Doelsch et al., 2020).This means that some SOM pools are active (microbial biomass, litter and dead roots, root exudates), whilst others are passive (humus, peat).Active carbon pools are available for microbial decomposition and then partitioned between CO 2 and CH 4 , where passive carbon pools decompose very slowly.Vascular plants generally have faster decomposition rates than mosses (Graf and Rochefort, 2009); therefore, vascular plants contribute to only one of the two passive SOM pools (humus, Table S2), whereas moss PFTs contribute to both passive SOM pools (humus and peat).The decomposition of each SOM pool is calculated assuming first-order rate kinetics: where SOM pools are represented by the subscript s, Q [kg C m −3 ] represents the mass of organic carbon in each SOM pool, and ke [d −1 ] represents the decomposition rate for each SOM pool adjusted by an environmental correction factor (Eq. S18, van Huissteden et al. ( 2006)).
In the above equation, SD [kg C m −3 d −1 ] represents the total carbon lost from all SOM pools.A fraction of the decomposed carbon from the SOM pools (litter and dead roots, root exudates, peat) is transferred (mineralised and reincorporated) into microbial biomass and humus, and the remaining fraction of SD is transferred into CO 2 .The CO 2 flux from the decomposition of SOM is calculated per soil layer: where FMI [-] refers to the fraction of SOM transferred to the microbial biomass pool, calculated using the PFT parameter MI [-] (Table 1); HU [-] refers to the fraction of SOM transferred to the resistant humus pool; Rox [µM CH 4 m −3 d −1 ] represents the portion of CH 4 oxidised to CO 2 (Eq.S27); and MC [kg C µM CH −1 4 ] represents the conversion factor (from µM CH 4 to kg C).

Methane processes
The daily CH 4 flux is dependent on the production and oxidation of CH 4 , as well as on the three transport mechanisms (diffusion, ebullition, and plant-transported CH 4 ).The CH 4 flux at the soil surface (Eq.29) is calculated by summing the three transport mechanisms (diffusion, ebullition, and planttransported CH 4 ).The CH 4 concentration of each soil layer (Eq.S19) is calculated before summing all transport mechanisms at the soil surface to obtain the net flux.Methane processes were adapted from the Peatland-VU model (van Huissteden et al., 2006), originally described in Walter and Heimann (2000).] is the ebullitive flux (Eq.S24).Methane production (Eq.S26), oxidation (Eq.S27), ebullition (Eq.S25), and diffusion of CH 4 through the soil (Eq.S20) remain as described in van Huissteden et al. (2006), originally adapted from Walter and Heimann (2000).
Plant-transported CH 4 is calculated for each PFT.There are two mechanisms which determine the amount of CH 4 lost via plant transport.Firstly, the mass and distribution of the root system play a role in determining how much CH 4 is taken up into the plant tissue.Thereby, a dense or large root system enables more CH 4 to enter the plant tissue.When CH 4 passes through the oxic zone around the root tips, a fraction of CH 4 is consumed by rhizospheric oxidation (Schipper and Reddy, 1996).This is represented by the unitless PFT parameter, PlOx (Eq.31).Secondly, the amount of CH 4 transported through the plant tissue and released to the atmosphere is determined by its aerenchyma.Plants with large aerenchyma are efficient transporters of CH 4 .The PFT parameter vP [-] describes the plant's ability to conduct CH 4 through aboveground plant tissue (Table 1).Shrubs and trees generally do not have aerenchyma, whereas grasses and sedges can have large or small aerenchyma (Ström et al., 2005;Walter and Heimann, 2000).The values for these PFT parameters are taken from the literature and are cited in Table S3.
In the above equation, ] is a rate constant with a value of 0.24 (taken from Walter and Heimann (2000)), f (z, p) [m −1 ] represents the exponential root distribution (Eq.17), and C CH 4 [µM m −3 ] represents the CH 4 concentration.The rate of plant-transported CH 4 is integrated over the depth of the root zone to obtain the flux at the surface (Eq.31): where Fpl represents the total plant-transported CH 4 flux [µM m −2 d −1 ].

Harvest scheme
If the harvest scheme is activated in the model input file, PFTs taller than the prescribed harvest height are harvested (mowed) at the prescribed date.This is a relevant feature for agricultural (e.g.Knox et al. (2015)) or other managed peatlands (e.g.Evans et al. (2021)).The harvest height and days are, therefore, optional prescribed model parameters.Living biomass decreases according to the amount of biomass harvested because biomass is assumed to be uniformly distributed with height and is not partitioned into organs.LAI is recalculated (Eq.8), and the PFT height is set to the harvested height.A fixed fraction of the harvested material is assumed to be lost during the harvest process, remains uncollected in the field, and is added to the litter layer.This fraction can also be set to zero.

Two peatland sites
With this study, the PVN model simulates two peatland sites in the Netherlands, the Horstermeer site and the Ilperveld site (Fig. S1).The Ilperveld site (52 ) is currently a nature recreation area that is a former raised bog complex that was drained to be used as agricultural pasture and is frequently exposed to manure fertilisation (van Geel et al., 1983;Harpenslager et al., 2015).Since the early 2000s, the Ilperveld site has undergone restoration efforts which included raising the water level, removal of the fertilised and nutrient-rich top soil, attempts to re-introduce Sphagnum, and water quality management.
The vegetation consists of brown mosses, Sphagnum, and grasses (Poaceae family).Since restoration began, the site has been mown twice a year, in June and September.Vegetation profiles show layers of intact Sphagnum and Carex peat, and, unlike undisturbed peatlands, the top layer has undergone greater decomposition due to land management since drainage (Harpenslager et al., 2015).The Horstermeer site (52 ) lies on the Horstermeer polder and is a former drained agricultural peat meadow that has not been used since the 1990s, when the water level was also raised.It was used for grazing and was exposed to manure fertilisation until the 1990s.The Horstermeer site is now a semi-natural fen containing very heterogeneous vegetation, including reeds, grasses, and small shrubs, and is not subject to mowing or other land management practices (Hendriks et al., 2007).Vegetation consists of different types of grasses and sedges (the dominant species are Holcus lanatus, Phalaris arundinacea, and Glyceria fluitans) and reeds (Phragmites australis and Typha latifolia).The Horstermeer polder is subject to strong seepage of mineral-rich groundwater from surrounding lake areas and Pleistocene ice-pushed ridges (Hendriks et al., 2007).The Horstermeer polder was a freshwater lake that was drained as part of large-scale land reclamation project completed in 1888.

Model calibration
The model was calibrated to reproduce fluxes that fall within the spread of observed in situ chamber measurements measured at the Horstermeer and Ilperveld peatland sites (Sect.2.2).The PVN model simulates processes at a daily time step.We ran the model using 28 years (1990-2017, inclusive) of input data (Sect.2.7) for the Horstermeer and Ilperveld sites.The length of the model spin-up was 5 years, determined by the time taken for the SOM pools, belowground CO 2 , and belowground CH 4 concentrations to stabilise (Fig. S4).Thereby, the first 5 years of model simulations (1990)(1991)(1992)(1993)(1994)(1995) are considered as the spin-up period.
Daily CO 2 and CH 4 fluxes measured at the Horstermeer and Ilperveld sites between 2015 and 2017 were used to calibrate the model.Unfortunately, there were not enough data to split the observational data into separate datasets for calibration and validation.A Monte Carlo analysis was performed separately for each site to calibrate 13 model parameters (Table S4).Parameters without available observational data were included in the model calibration process.The Kling-Gupta efficiency (KGE) metric was used to measure the agreement between simulated and observed CO 2 and CH 4 fluxes (Gupta et al., 2009;Kling et al., 2012).The KGE approach is a threedimensional decomposition of the Nash-Sutcliffe efficiency (NSE) measure and evaluates temporal dynamics, bias, and variability Eq. ( 28).The KGE metric has been used to assess the ability of carbon flux models (Tramontana et al., 2016;Zhao et al., 2020), hydrological models (Dick et al., 2015), and meteorological reanalysis datasets (Chaney et al., 2014;Muñoz-Sabater et al., 2021) to reproduce in situ observations.The calibrated model input values are provided in Tables S6 and S7 for the Horstermeer and Ilperveld site simulations, respectively.
The CO 2 results impact the CH 4 results much more than the CH 4 results impact the CO 2 results; therefore, we first ensured that the parameters impacting the photosynthesis and the above-and belowground growth and respiration schemes reproduced fluxes that fell within the spread of observed CO 2 fluxes (NEE).These were the MolAct, HalfSatPoint, and VegTScalingFactor parameters.Next, the CH 4 scheme was calibrated to reproduce fluxes that fell within the spread of observed CH 4 fluxes.This involved calibrating the remainder of the parameters highlighted in Table S4.Even though the amount of photosynthesis and living biomass does not directly impact the CH 4 production, which primarily occurs in the soil and aboveground litter layers, these processes are precursors to root and shoot growth, respiration, and senescence, which directly impact simulated CH 4 fluxes.After optimisation of the CH 4 fluxes, the PFT parameters (Table S3) were manually adjusted to bring the PFT biomass fractions (PFT biomass as a fraction of total biomass) in line with observed aerial cover fraction ratios.The calibrated model parameters and the necessary input files used to simulate the two peatland sites evaluated in this study are accessible from the Bitbucket repository: https://www.bitbucket.org/tlippmann/pvn_public (last access: 15 November 2023).

Testing the PVN model
To understand the sensitivity of net CO 2 and CH 4 fluxes to PFT-dependent processes, we conducted several model simulations using modified input data.Air temperature, water table, radiation, harvest, and the swapping of PFTs between site simulations were chosen to be used for the sensitivity testing because they are key environmental drivers of CO 2 and CH 4 emissions in peatlands.We tested the sensitivity of PFT processes by varying these inputs one by one (Table 2).
To understand how the new model mechanisms affect emissions, we performed additional simulations with altered model algorithms and compared these to the original model simulations calibrated for the Horstermeer and Ilperveld sites (Table 3).For example, the contribution of competition for shading to the overall simulation result is quantified by comparing an altered simulation where incoming photosynthetically active radiation (PAR) is independent of shading (e.g.fractional par or FPAR = 0.25 for a simulation with four PFTs) to the original model simulations (FPAR_CONST).We calculated the relative difference of the simulation with shading minus the simulation without shading.Similarly, we compared simulations with and without plant-transported CH 4 (CH4_OLD_CF), with and without dynamic BF (CF_CONST), and with and without variable plant height (HEIGHT_CONST).
In order to demonstrate that the PVN model reproduces CH 4 and CO 2 fluxes within the spread of observed fluxes when driven by realistic input data, we compared the calibrated model simulation results and measured CH 4 and CO 2 fluxes for the Horstermeer and the Ilperveld field sites (Sect.2.2).
We compare the CH 4 and CO 2 fluxes of the calibrated model simulation results against the CH 4 and CO 2 fluxes  4) will use the PFTs observed at the Horstermeer site (Typha, tall grass, sedges, brown moss), while the model simulation driven by the Horstermeer input data will use the PFTs observed at the Ilperveld site (short grass, tall grass, Sphagnum, brown moss).4. Attempts to run the Peatland-VU model with new calibrated parameters did not yield results on the same order of magnitude as the observations.Therefore, it was necessary to use different model parameterisations for the PVN and Peatland-VU models.

Flux measurements
Carbon dioxide and CH 4 fluxes were measured using two to four automated flux chambers (AC) and the Ultra-Portable Los Gatos Gas Analyser Model 915-001.Chambers were cylindrical, 30 cm wide and 40 cm in height, made of transparent acrylate, equipped with a fan, and installed in the field using collars.Where necessary, vegetation was folded gently to fit inside the measurement chambers.Collars were removed from the field between sampling campaigns, which minimises disturbance which can lead to potential biases in the observations.This also potentially introduces uncertainty with regard to the precise measurement location.The CO 2 and CH 4 concentrations were measured for 150 s intervals whilst the chamber was closed.Each chamber was measured on rotation so that a new chamber was measured every 15 min.Measurements were recorded continuously during the day and night for a week at a time, during which time the AC system was moved to another site.We note that, due to the labour-intensive nature of accumulating chamber observations consistently through time, these observational datasets do not offer complete temporal continuity, creating an intermittency bias.From these data, the hourly average CO 2 (net ecosystem exchange) and CH 4 fluxes were calculated for each day.We compared calibrated site simulations against observed daily average CO 2 and CH 4 fluxes.To visualise the daily variability, standard deviations were derived from the hourly fluxes.The values for all GHG emissions are expressed as CO 2 equivalents [kg CO 2eq m −2 yr −1 ] and are calculated as where GWP 20 = 80.8 as 1 kg CH 4 = 80.8 kg CO 2 eq.over a 20-year time horizon, and GWP 100 = 27.2 as 1 kg CH 4 = 27.2 kg CO 2 eq.over a 100-year time horizon (Masson-Delmotte et al., 2021).

Input data preparation
The PVN model is driven by daily air temperature (T ), water level (WL), radiation, a general model parameter input file (Table S4), and a soil parameter input file (Table S5).

Climatological input data
Daily temperature and radiation data, measured at Schiphol, the nearest KNMI weather station, were used as climate input data for both sites (accessed via https: //www.knmi.nl/nederland-nu/klimatologie/daggegevens,last access: 18 May 2022) (Fig. S3).The annual average rainfall at Schiphol was 850 mm yr −1 over the period 1990-2020, with 30 % of the rainfall falling in summer and autumn and 24 % falling in winter, with the remainder falling in the spring.The average daily temperature between 1990 and 2019 was 9.4

Soil profile input data
The model generates a soil horizon representation using soil layers of 10 cm thickness.The generated soil horizon uses properties such as the dry bulk density (DBD), SOM ratio, sand content, and C:N ratio specified in the soil profile input data (Table S5).The number and depths of the site's soil horizons can be adjusted in the soil input file.The PVN model requires input parameters for each PFT, as discussed in Sect.2.3.Soil profile data from the Horstermeer and Ilperveld field sites were collected in 2015 and 2016 and include DBD, C content, SOM content, sand and clay content, and pF curve (Tables S8 and S9 for the Horstermeer and Ilperveld site simulations, respectively).

Water level input data
Water level input data were sourced from the Dutch hydrological model, the Netherlands Hydrological Instrument (NHI) (De Lange et al., 2014), which has a reasonably high spatial resolution (250 m ×250 m).One aim of developing the PVN model is to eventually develop a model of all Dutch peatlands in conjunction with the NHI product.For this reason, the NHI product is used in this application of the model.The NHI water level output was converted to relative surface height using a 5 m ×5 m digital elevation map of the Netherlands, Actueel Hoogtebestand Nederland (Alhoz et al., 2020).It is possible to use in situ water levels as input data for the model, but these data were, unfortunately, unavailable for the duration of the simulation.The input data used for both sites are accessible from the Bitbucket repository: http://www.bitbucket.org/tlippmann/pvn_public (last access: 15 November 2023).

Results
When describing the annual CO 2 , CH 4 , and GHG values, we opt to use the term emissions, e.g. annual GHG emissions, whereas, when describing daily values, we opt to refer to these as fluxes, e.g.daily GHG fluxes.

Model sensitivity to input data
To understand the response of the modelled PFT processes to input data, we ran simulations with modified water levels (Figs. 3 and S6), temperature (Figs. 2 and S5), and radiation (Fig. S7) input and harvest schemes (Fig. 4).The modified input data are summarised in Table 2, and the results of these sensitivity tests are summarised in Table 5.These results are indicative of the model's mechanistic responses rather than being projections of how PFTs might respond under varied environmental conditions.To show how different inputs impact model processes, we present the soil (respiration) CO 2 emissions (Fig. 3), plant-transported CH 4 (Fig. 2), and aboveground biomass (Fig. 4).In the PVN model, the abundance of each PFT varies through time depending on the favourability of growing conditions.Therefore, an increase in CO 2 or CH 4 emissions may be due to increased abundance (i.e.enhanced biomass) or enhanced transport efficiency.To disentangle this difference, the CO 2 and CH 4 emissions for each PFT are plotted as a fraction of litter and root mass.
Increased air temperatures had a positive effect on both plant-transported CH 4 emissions (Fig. 2) and litter and root mass at both sites (Fig. S5).Short and tall grasses showed similar responses to increased air temperatures by producing large CH 4 emissions per kilogram of litter and root mass.Brown mosses showed little variation between the temperature experiments for the Ilperveld site but showed a decrease in emissions with warming temperatures per kilogram of litter and root mass at the Horstermeer site.Sphagnum similarly showed a decrease in CH 4 emissions with warming temperatures per kilogram of litter and root mass at the Ilperveld site.This decrease is because moss PFTs have strict idealtemperature growth limits (T max and T min in Table 1) and were limited by warming temperatures.Whilst belowground, CH 4 concentrations increased with warming temperatures, and the biomass, litter, and root mass of moss PFTs did not increase with warming temperatures.
Belowground CO 2 emissions were impacted by changing water levels (Fig. S6).Previous studies have found that belowground CO 2 production tends to increase with low water levels due to enhanced potential for aerobic CO 2 production (Knox et al., 2015).The results of the Ilperveld site sensitivity simulations showed that belowground CO 2 production increased with low water levels, likely due to enhanced potential for aerobic CO 2 production.However, the results of the Horstermeer site sensitivity simulations showed the converse: that the net CO 2 (Table 5) and belowground CO 2 production increased with high water levels.We simulate that, with high water levels, the reduced aerobic CO 2 production can be exceeded by the enhanced oxidation of CH 4 into CO 2 .The large amounts of CH 4 oxidised into CO 2 in the Horstermeer site simulation are due to the very degraded peat present at the site (represented by low soil OM content in the soil input file) and the strong upwelling of rich groundwater at the Horstermeer site (represented by the calibratable model parameter, MolAct (see Sect. 2.4), which influences the sensitivity of aerobic CO 2 production).The large observed CH 4 emissions at the Horstermeer site are partially due to high CH 4 concentrations in the upwelling water.Furthermore, the large root systems of plants such as Typha, sedges, and tall grasses have greater potential to access and transport stores of belowground gases (represented by the PFT root depth and mass).The conflicting response of the tall grass PFT in the Ilperveld and Horstermeer simulations shows that PFTs may respond differently to changing water levels at different sites.
Increasing the frequency of harvests led to a strong negative effect on vascular plant biomass and a small positive effect on moss plant biomass (Fig. 4).The biomass of nonmoss PFTs is strongly impacted by the occurrence of harvests, as indicated by the pause in biomass accumulation after harvest.However, by reducing tall vegetation, moss species have greater access to sunlight and, therefore, gain an advantage.For this reason, we saw the biomass of moss PFTs increase with more frequent harvests.In the Horstermeer site simulation, the greatest effect on biomass was between no harvests and the once-per-year harvests.In the Ilperveld site simulation, the effects of harvests on biomass increased somewhat linearly according to the frequency of harvest events.We suspect that this is due to the inclusion of different PFTs in the two site simulations.In the Horstermeer site simulation, three PFTs have the capacity to grow above the harvest height (the Typha, tall grass, and sedge PFTs), whereas in the Ilperveld site simulation, only tall and short grasses have the potential to grow beyond the harvest height, thereby limiting the potential effect harvests can have on the PFTs present.Furthermore, the growth of the short grass PFT is height limited to 0.3 m.Overall, total biomass was reduced with more frequent harvest regimes.
It is important to note that, whilst CO 2 emissions were reduced by increasing the frequency of harvests, these emissions do not account for the off-site decomposition of harvested biomass.Methane emissions were slightly enhanced if harvests occurred in comparison to no harvest events for Horstermeer site simulations, whilst the frequency of harvests did not impact emissions (Table 5).Similarly, enhanced https://doi.org/10.5194/gmd-16-6773-2023 Geosci.Model Dev., 16, 6773-6804, 2023  CH 4 emissions occurred with increased harvest frequency for the Ilperveld site simulations.Spikes in CH 4 fluxes transported by vascular PFTs occurred after harvest events in both the Horstermeer and Ilperveld simulation results (not shown), contributing to enhanced CH 4 emissions for both the Ilperveld and Horstermeer site simulations.The impact of fewer or no harvest events led to variable impacts on CH 4 emissions for the Ilperveld site simulations, where a single harvest led to slightly reduced emissions, and no harvests led to slightly enhanced emissions.In the Ilperveld site simulation without harvest events, vegetation became dominated by vascular PFTs that are efficient transporters of CH 4 , leading to enhanced CH 4 emissions.

Assessment of model mechanisms
To understand the role of isolated model mechanisms, we modified the model code to disable the functions responsible for reproducing the vegetation dynamics within the model (Fig. 5).Unlike the other simulations assessed throughout this paper, the simulation results shown in Fig. 5 begin in the year 1990, i.e. without the use of a spin-up period.Removing the spin-up period showed that the modified model simulation results produce similar emissions in the first year of the simulation (1990) and allows an assessment of the trajectory of deviation.Disabling the shading scheme (simulation PVN_HEIGHT_CONST) or biomass fraction scheme (simulation PVN_CF_CONST) led to only slightly enhanced CO 2 emissions, whereas disabling the FPAR scheme (simulation PVN_FPAR_CONST) led to large CO 2 emission differences.Surprisingly, the difference for the PVN_FPAR_CONST simulation is opposite in sign for the two site simulations and is larger for the Ilperveld simulation.This means that maintaining constant FPAR led to a small enhancement in CO 2 emissions in the Horstermeer simulation but a large reduction in CO 2 emissions for the Ilperveld simulation.These results show that FPAR plays a large role in simulated CO 2 emissions.The results of the Ilperveld PVN_FPAR_CONST simulation also showed that the FPAR function has the potential to introduce large variability into the emission results.This is interesting to note because the PVN model showed limited skill in reproducing the CO 2 emissions at the Ilperveld site.These results indicate that the function calculating FPAR plays a driving role in CO 2 emissions but particularly at the Ilperveld site.Further model developments may investigate ways to improve the representation of FPAR in the model.The PVN_FPAR_CONST simulations also led to enhanced CH 4 emissions for the Ilperveld simulation.It is likely that CH 4 production was enhanced due to increased stores of CO 2 .
The use of the Peatland-VU CH 4 scheme (PVN_CH4_OLD_CF) led to large differences in CH 4 emissions for both the Horstermeer and Ilperveld simula-tions in comparison to the PVN model results.The CH 4 emissions of the model simulations that use the Peatland-VU CH 4 scheme (simulation PVN_CH4_OLD_CF) were small when compared to the CH 4 emissions of the PVN model for both model simulations.This indicates that the PFT modifications to the CH 4 scheme have led to substantial impacts on modelled CH 4 emissions.

Assessment of calibrated model simulations
Here, we describe the simulation results of the model calibrated at two field sites, Horstermeer and Ilperveld.We describe the net annual CH 4 and CO 2 emissions and GHG budgets (Fig. 6), as well as the simulated PFT dynamics as indicated by changes to LAI, aboveground biomass, litter mass, and PFT height and/or depth (Figs. 7 and S8).All net GHG values are expressed as CO 2 equivalents (CO 2 eq.).
The model simulation results indicate that the simulated annual mean net GHG emissions from the Ilperveld simulation were approximately half the emissions of the Horstermeer simulation.However, these model emission estimates do not consider off-site decomposition of harvested biomass.The model estimated that the 2015-2017 annual average net GHG emissions were 2.5 and 8.9 kg CO 2 eq.m −2 yr −1 for the Ilperveld and Horstermeer simulations, respectively (Table 6), calculated using the 20-year GWP.Using the 100year GWP, the 2015-2017 annual average net GHG emissions were 2.3 and 5.6 kg CO 2 eq.m −2 yr −1 for the Ilperveld and Horstermeer simulations, respectively.The model estimated that the 1995-2017 annual average net GHG emissions were 2.4 and 8.0 kg CO 2 eq.m −2 yr −1 for the Ilperveld and Horstermeer model simulation results, respectively (Fig. 6), calculated using the 20-year GWP.Using the 100year GWP, the 1995-2017 annual average net GHG emissions were estimated to be 2.3 and 5.2 kg CO 2 eq.m −2 yr −1 for the Ilperveld and Horstermeer model simulation results, respectively.
Assessment of the Horstermeer simulation showed that, on average, CH 4 contributed approximately half (52 %) of the net annual GHG emissions of the Horstermeer simulation, where CH 4 contributed 4.2 kg CO 2 eq.m −2 yr −1 and CO 2 emissions contributed 3.8 kg CO 2 eq.m −2 yr −1 , on average.Assessment of the Ilperveld simulation showed that CO 2 was the primary contributor to net GHG emissions, where CO 2 contributed the majority (92 %) of the annual GHG emissions (2.2 kg CO 2 eq.m −2 yr −1 of the total 2.4 kg CO 2 eq.m −2 yr −1 net GHG emissions).These model emission estimates neglect the off-site decomposition of harvested biomass.Therefore, CO 2 and CH 4 emissions contribute equally to the net GHG emissions in the Horstermeer simulation, whereas CO 2 emissions dominate the GHG emissions in the Ilperveld simulation results.

PFT dynamics
Here, we describe the living biomass, LAI, litter layer, biomass fraction, and height changes of the PFTs of the calibrated Horstermeer and Ilperveld model simulations (Figs. 7  and S8).Assessment of the aboveground biomass (top row of Fig. 7) shows that the tall grass (blue line), Typha, and sedge PFTs (red line) were abundant during the Horstermeer simulation.whereas the Ilperveld simulation was dominated by the short grass (green line), Sphagnum (pink line), and tall grass PFT (blue line).All plants showed seasonal variability.
The ratio of the litter layer to biomass is between approximately 1 : 4 and 1 : 3 for most PFTs [kg C].The Typha PFT is an exception, and the ratio is approximately 1 : 1.Overall, the sedge PFT showed comparable seasonal variability to the tall grass PFT whilst maintaining less biomass, smaller LAI, and shorter height throughout the Horstermeer simulation.The similar behaviour of the Typha, sedge, and tall grass PFTs was expected because the PFT input parameters repre-sent similar plant phenologies.Assessment of the size of the litter layer (first row of Fig. S8) showed that, in the Ilperveld simulation, the PFTs reached peak litter during autumn (September), whilst in the Horstermeer simulation, which is not mown, the litter continued to accumulate until January, where rates of decomposition exceeded accumulation.The LAI (second row of Fig. S8) displayed strong seasonal variability.Each year, the LAI of the short grasses reaches its maximum LAI value of 1.2.The tall grass PFT, whilst very competitive in the Horstermeer simulation, is less competitive in the Ilperveld simulation, partially due to the occurrence of harvests and partially because it is outcompeted by the fast-growing short grass PFT.Assessment of the Ilperveld simulation reveals that the short grass PFTs were constrained by the maximum-height parameter, MaxCanopyHeight.The tall grass PFT was not limited by MaxCanopyHeight in the Ilperveld simulation but was instead limited by the biannual mowing regime.PFT height showed strong seasonal variability for both simulations (third row of Fig. S8).The tall grass PFT constituted the tallest plants in the Horstermeer simulation until 2009, and its height was frequently limited by the PFT MaxCanopyHeight parameter.However, as the Typha PFT grew in biomass, the tall grass PFT appeared to have less access to sunlight as height and biomass values were reduced.The Typha and sedge PFTs were not limited by their maximum-height parameters.These changes in biomass fraction are also evident in the emissions.
The relative contributions of each PFT to the net annual CH 4 , CO 2 , and GHG emissions are shown in Fig. 6, where the CH 4 emissions refer to only the plant-transported CH 4 .The net CO 2 emissions for each PFT are the sum of the https://doi.org/10.5194/gmd-16-6773-2023 Geosci.Model Dev., 16, 6773-6804, 2023 photosynthesis minus respiration, the CO 2 produced by belowground aerobic decomposition of SOM, and a portion of CH 4 oxidised to CO 2 .The tall grass (red boxes), sedge (orange boxes), and Typha (purple boxes) PFTs are large transporters of CH 4 emissions in the Horstermeer simulation results.However, only the tall grasses and Typha compose the net CO 2 emissions in the Horstermeer simulation.Thereby, the tall grass PFT was the largest contributor to the net annual GHG emissions, followed by the Typha and sedge PFTs.The Ilperveld simulation results showed that the short grass PFT was the largest contributor to the net annual CH 4 , CO 2 , and GHG emissions.

Comparison of modelled and observed plant dynamics
We compare simulated PFT biomass fractions against observed aerial plant cover fractions (Fig. 8).For assessment against observational data, we compare model simulation results against observed fluxes by comparing time series, boxplots, and 1 : 1 scatterplots for CH 4 (Fig. 9) and CO 2 (Fig. 10).Gaps in observational data exist due to measurement collection limitations; therefore, the model comparison against observational data can only be shown for the days where observational data exist.Unfortunately, this means that the model was not assessed equally across all seasons or on the same days of the year at the two sites.A simple linear regression is used to compare the model simulation results and observational data using all days with available measurements.For these reasons, the 1 : 1 plots and R 2 linear regression results may only give a flavour of model performance.To understand the degree of uncertainty in the observational measurements, daily standard deviations were derived using the hourly fluxes (plotted as black error bars in Figs. 9 and 10).In each case, the model simulation results generally lay within the spread of observational uncertainty.
The observations indicated that both sites are annual sources of CH 4 and CO 2 and, therefore, net annual sources of carbon to the atmosphere.The Horstermeer site produced large annual mean CH 4 and CO 2 emissions in comparison to the Ilperveld site (Fig. 9 and Fig. 10).

Evaluation of plant composition dynamics
Plant cover fraction observations were made at the location of the chamber measurements and were not representative of the site's complete plant community composition.Although aerial cover fraction and biomass fraction (the ratio of PFT biomass to total biomass) are not the same, changes in plant composition are depicted in both representations.
In 2006, the chamber measurement location at the Horstermeer field site was composed of tall grasses (50 %), sedges (40 %), Typha (5 %), and brown mosses (5 %) (left panel in Fig. 8).The Horstermeer simulation results show good agreement with the observations but overestimated the amount of tall grasses (66 %) and underestimated the amount of sedges (40 %).In 2016, a decade later, the amount of tall grasses remained consistent, whilst the amount of Typha had increased by 10 %; 1 year later, in 2017, the vegetation had not undergone changes proportionally.Parallel to the observations, the Horstermeer simulation results estimated that the tall grass PFTs decreased to 60 % from 2005 onwards, whilst the biomass fractions of the Typha and sedge PFTs increased.Overall, the Horstermeer simulation overestimated the biomass fraction of the tall grass PFT and underestimated the proportion of the sedge and Typha PFTs.Model estimates of year-to-year PFT biomass changes were of the same sign and similar magnitude as in situ observations.
In March 2016, the chamber measurement location at the Ilperveld field site hosted short grasses (50 %) and tall grasses (50 %).The model overestimated the amount of short grasses (80 %), underestimated the amount of tall grasses (5 %), and overestimated the amount of Sphagnum (10 %).The Omhoog met het Veen (Raising the Peat) project delivered on-site management attempts to initiate Sphagnum growth by hand dispersing living fragments of Sphagnum spp.from a nearby donor site between 2013 and 2015 (Geurts and Fritz, 2018).For this reason, we expected that the model may not match the development of Sphagnum at the Ilperveld site.In October 2017, the vegetation shifted to be composed of short grasses (50 %) and tall grasses (25 %), Sphagnum (15 %), and brown mosses (10 %); 1 month later, in November 2017, the Sphagnum was no longer visible (0 %), brown mosses remained (10 %), and the site was dominated by short grasses (80 %).The model estimated that the short grass and Sphagnum PFTs remained consistent into 2016 and 2017, whilst the tall grass PFT was reduced, and brown mosses increased slightly.Whilst the model simulations ended in 2017, we saw that, in October 2018, the vegetation remained constant at both sites.

Evaluation of simulated CH 4 fluxes
The time series presented in Fig. 9 shows the behaviour of the Horstermeer simulation CH 4 flux results (purple line), the observed mean daily fluxes (black dots), and the spread of the hourly observed fluxes (black error bars).Whilst the Horstermeer simulation reproduced the seasonal variability of the observed CH 4 fluxes, the boxplots showed that the simulation results (purple box) tended to overestimate the CH 4 fluxes.Overall, the Horstermeer simulation showed a robust pattern of variability when compared with the observations (R 2 = 0.7) whilst overestimating the magnitude of observed fluxes.Assessment of the Ilperveld model simulation showed that the model was able to reproduce the observed CH 4 fluxes and followed the pattern of variability when compared with the observations (R 2 = 0.8).The summer of 2015 is an exception where the simulated results showed an increase in CH 4 fluxes larger than the observed CH 4 fluxes.Assessment of the boxplots showed that the simulated CH 4 fluxes (green box) are of a similar mean and spread to the observed fluxes (purple box).

Evaluation of simulated CO 2 fluxes
The boxplots showed that the PVN Horstermeer simulation reproduced the median and range of observed daily CO

Comparison to the PEATLAND-VU model
To understand the impact of including vegetation dynamics, we compare the results of the new PVN model against the results of the pre-existing Peatland-VU model (Fig. 9 and CO 2 in Fig. 10).The simulation results are summarised in Table 6.Overall, the PVN model estimated the net annual CH 4 , CO 2 , and GHG emissions to be larger than the emissions estimates made by the Peatland-VU model.The Peatland-VU model estimated the annual mean 2015-2017 GHG emissions to be 1.3 and 5.9 kg CO 2 eq.m −2 yr −1 for the Ilperveld and Horstermeer simulations, respectively, calculated using a 20-year GWP.When calculated using a 100-year GWP, the Peatland-VU model GHG emission estimates for the Horstermeer simulation were 3.8 kg CO 2 eq.m −2 yr −1 (for both periods of 2015-2017 and 1995-2017).The Peatland-VU GHG emission estimates for the Ilperveld simulation were 1.

Discussion
We have developed the PVN model, a new dynamic vegetation-peatland-emissions model capable of understanding the role dynamic PFTs play in CO 2 and CH 4 emissions in peatlands.We tested the sensitivity of simulated PFT processes to changing environmental parameters and investigated the impacts of the new schemes introduced into the model that attempt to replicate competition between vegetation types.Here, we discuss potential sources of uncertainty, https://doi.org/10.5194/gmd-16-6773-2023 Geosci.Model Dev., 16, 6773-6804, 2023  Model 2015-20171995-20172015-20171995-20172015-20171995-2017 Horstermeer PVN 8.88 (5.56)  both in the observational data used to evaluate the model results and in the chosen model input parameters.Secondly, we discuss the processes in the model that allow the representation of dynamic vegetation and the ability of these processes to respond to changing environments.Lastly, we discuss how the new PVN model compares to its two parent models, the NUCOM-BOG model and the Peatland-VU model, as well as the one other site-specific GHG emissions peatland model that uses dynamic PFTs.

Site
4.1 Sources of uncertainty

Input parameters
It is important to note that the Peatland-VU, NUCOM-BOG, and PVN are heavily parameter-dependent models.The Peatland-VU model has been shown to reproduce observed fluxes using widely different parameter sets, which means that the Peatland-VU model has a strong equifinality of parameterisations (van Huissteden et al., 2009) because there are simply not enough data available to constrain all model dynamics.One aim of introducing PFTs into the Peatland-VU model was to develop a model with a greater dependence on observational data (measurable PFT traits) and less dependence on optimised parameters, reducing the equifinality of the model.It is important that improvements to model processes capture the critical processes but as simply as possible to minimise problems that arise due to the equifinality of parameterisations (Beven and Freer, 2001).The introduction of PFTs allowed several Peatland-VU parameters that were previously calibratable to become observation-informed parameters whilst introducing few new parameters; thereby, the net result is a reduction in the breadth of the parameter space.

Site heterogeneity and chamber measurements
We compare the findings of this study against other studies that have assessed observed CH 4 fluxes at the Horstermeer site and discuss uncertainties accompanying the chamber measurement technique.The sites simulated in this study pose challenges because they are degraded peatlands where easily decomposable carbon is likely to have been mineralised (Dorrepaal et al., 2007;Järveoja et al., 2013), peat has been artificially removed for centuries (Erkens et al., 2016) (Hendriks et al., 2007).
The different chamber measurement locations used by the two studies may account for some of the observed differences.Heterogeneous vegetation and heterogeneous water levels relative to the surface are known to impact both automated and manual flux chamber measurements.For this reason, observational measurements are impacted by the physical placement of flux chambers in the field, leading to potential measurement bias (Speckman et al., 2015;Baldocchi, 2003).At very heterogeneous sites, such as the Horstermeer site, flux strengths vary due to micro-topography (Wania et al., 2010), and chamber measurements have been reported to vary significantly within one site, which may explain differences between studies.
The Horstermeer site has vegetation standing taller than 1 m.At times, it was necessary to consider the vegetation height when selecting the chamber location to ensure vegetation (even when folded) could fit within measurement chambers.Field measurements that exclude areas covered by tall vegetation may result in a significant underestimation of CO 2 or, particularly, CH 4 fluxes.The absence of tall-vegetation measurements limits the capacity to test model representations of tall-vegetation processes, restricting the ability to predict changes in CO 2 and CH 4 fluxes in the presence of tall vegetation (Pangala et al., 2013).Due to the labour-intensive nature of accumulating chamber observations consistently through time, these observational datasets do not offer complete temporal continuity, creating an intermittency bias.The high cost of AC meant that sites could not be measured simultaneously, leading to an interrupted sampling regime that may bias CO 2 and CH 4 flux estimates (Morin et al., 2014a(Morin et al., , 2017)).Most chamber measurements are taken during the plant-growing season, assuming that the winter fluxes are negligible, which has been shown to not always be the case (Morin et al., 2014b).Future studies can benefit from continuous AC measurements.

On the efficacy of simulating dynamic vegetation
The PVN model was developed by building upon the functionality and structure of the Peatland-VU model whilst incorporating vegetation dynamics from the NUCOM model.The model has incorporated vegetation dynamics and enhanced the Peatland-VU model's existing carbon-cycling processes.Competition is based on water table depth, temperature, vegetation height, and shading.To verify that the model dynamics are robust and to understand the sensitivity of the PFTs, we performed model sensitivity simulations.
Considering that the short grass, Sphagnum, and brown moss PFTs share similar PFT parameters, these three PFTs can respond somewhat similarly.Whilst the short grass PFT is a non-moss PFT, its parameters are not dissimilar to those of moss PFTs.However, the short grass PFT quickly increases in biomass due to its broad range of temperatures and water levels for growth.This means that the short grass PFT provides strong competition against other PFTs.Even though the short grass PFT is height limited, its quickly increasing biomass allows increasing access to PAR, which leads to large amounts of plant respiration, root growth, and net CO 2 fluxes when compared to the Sphagnum and brown moss PFTs.With only a shallow root system (maximum 0.1 m), moss PFTs have limited abilities to transport belowground CO 2 , and, as expected, the total belowground CO 2 flux is small for mosses.Whilst mosses do not have root structures in reality, we allocated moss PFTs to have a presence in the top 10 cm of the soil layer because, in the presence of bryophytes, there is often no clear separation between the living-moss layer and the soil surface.In this way, we intended to replicate a transition zone.Key differences in the parameters between short grasses and brown mosses are that short grasses are not considered to be a moss PFT (relevant for height growth and light interception).In the model, moss PFTs have large CH 4 vP values and low leaf maintenance respiration coefficients and biomass senescence values (Table 1).Whilst these differentiations have been somewhat effective, future model versions might consider further ways of distinguishing moss PFTs (especially Sphagnum).The presence of Sphagnum in SOM increases the acidity of the soil.By influencing the acidity of the soil and limiting the nutrient availability, Sphagnum gains an advantage over other plant types because Sphagnum flourishes in nutrient-poor conditions (Moore et al., 2007).A useful addition to future model versions may be to adapt the living-moss layer to be incorporated into the soil layer, altering the height of the land sur-face and simulating hummock and hollow microtopography, as well as impacting corresponding soil properties (e.g.pH, DBD).
Largely, decomposition of the peat reservoir led to enhanced CO 2 fluxes due to a thick aerobic layer with low water levels.Modelled photosynthesis and plant respiration are dependent on both temperature and water levels.This enables an assessment of the impacts of water availability and extreme temperatures on plant type.Future model applications may consider the relationship between water availability and plant dynamics and, particularly, the impacts of drought on plant photosynthetic capacity, respiration, soil respiration, and CH 4 production and oxidation.

Impacts of changing temperature input
Studies show that, whilst both CH 4 production and oxidation rates are enhanced by warming, the net CH 4 flux increases with warming because CH 4 production increases at a rate faster than oxidation (Granberg et al., 1999).As expected, the PVN model simulated enhanced CH 4 emissions under simulations driven by warmer temperatures and simulated reduced CH 4 emissions under simulations driven by cooler temperatures.Sphagnum, tall grasses, and brown mosses showed unexpected results because they released fewer CH 4 emissions under warmer simulations.This may be indicative of the narrow temperature limits of Sphagnum moss (Gunnarsson et al., 2004).The impacts of temperature on model processes are 3-fold.Firstly, the amounts of photosynthesis and plant respiration performed are dependent on the ideal and tolerated PFT growth temperatures.Secondly, the amount of litter converted to belowground SOM reservoirs is dependent on soil temperatures, where warmer soil temperatures lead to larger amounts of litter being converted to belowground reservoirs.Thirdly, decomposition of belowground SOM is dependent on soil layer temperature (as well as pH, saturation, etc.), where soil layers closer to the surface are warmer.Thereby, temperature influences the PFT abundance, the size of litter, and the belowground SOM reservoirs available for decomposition, along with the efficiency of belowground SOM decomposition in the model.The results of our sensitivity analyses are in agreement with field studies which have found that CH 4 emissions are typically higher when dominated by Carex compared to Eriophorum or Juncus (Ström et al., 2005;Jackowicz-Korczyński et al., 2010).This is likely partly due to the presence of aerenchyma and partly due to differing litter quality and rates of carbon turnover (Christensen et al., 2003;Ström et al., 2015).

Belowground decomposition
Enabling different PFTs to contribute to, oxidise, and decompose different belowground SOM pools impacted simulated CO 2 and CH 4 fluxes.Decomposition in the PVN model is dependent on the decomposition rates of different https://doi.org/10.5194/gmd-16-6773-2023 Geosci.Model Dev., 16, 6773-6804, 2023 PFTs.Decomposition rates have been found to differ between forbs, graminoids, deciduous shrubs, and evergreen shrubs (Dorrepaal et al., 2006(Dorrepaal et al., , 2007(Dorrepaal et al., , 2009)).The peat SOM pool of moss PFTs contributes to CO 2 and CH 4 fluxes because (Sphagnum) mosses are the primary peat-contributing plant, and mosses (especially Sphagnum) have slow decomposition rates (Hobbie et al., 2000).Moss PFTs are the only PFTs able to contribute to the peat SOM pool, which means that the CH 4 fluxes arising from decomposition of the peat SOM pool are only transferred to the surface by moss PFTs.Future modelling efforts could work to improve the representation of peat decomposition, whereby CO 2 fluxes resulting from the decomposition of peat can be transferred to the surface by both moss and non-moss PFTs.Mosses are prescribed to have roots of 0.1 m at maximum the model is initialised, and this remains constant throughout the model simulation.Mosses do not have an aboveground litter layer, and instead, their living biomass after senescence is added directly to the belowground SOM.

Root distribution representation
Plant-transported CH 4 and aerobic CO 2 production processes are dependent on root mass and are independent of aboveground biomass.In the model, the belowground CO flux is comprised of CO 2 produced by peat, root exudates, litter, roots, microbial biomass, humic matter, and CH 4 oxidation.Root traits play an important role in species competition and processes such as leaf : root allocation, turnover, root stocks, and root distributions have been shown to be dependent on climate, species, and land cover type (Smithwick et al., 2014), particularly in Arctic and boreal systems (Iversen et al., 2015).Root exudation plays an important role in the supply of substrates that can later be metabolised into CH 4 (Aulakh et al., 2001;Waldo et al., 2019), where the fraction of belowground production that consists of exudates (REX) was an important parameter impacting CH 4 production in the model.Next to this, the parameter representing root aerenchyma (PlOx) played a role in the oxidation of CH 4 (Walter and Heimann, 2000).These processes, as well as plant-transported CH 4 , are only possible from soil layers with roots present (Bansal et al., 2020;Walter and Heimann, 2000).For this reason, the parameter representing maximum root depth (MRD) played a role in the production, oxidation, and transport of CH 4 , where the relative impact of each of these processes on surface CH 4 fluxes is dependent on PFT properties.
Root distribution structural representation is important to reliably simulate CO 2 and CH 4 fluxes in the model.Land surface models have, for the most part, used exponential relationships to describe root distribution (Smithwick et al., 2014;Zeng, 2001).Advances have been made in developing the knowledge and observational data of root distributions in boreal peatland systems.Whilst the exponential relationship is representative for several peatland plant types, an alterna-tive root representation for the exponential relationship may be relevant for certain peatland plant types (Clemmensen et al., 2013;Iversen et al., 2015).Future model versions may consider introducing alternative root representations.

The impact of harvests on plant competition
The inclusion of harvest has proven to be necessary to reproduce the seasonal variability of emissions in grasslands and croplands, where crop harvests occur (Van den Hoof et al., 2011).Whilst CO 2 emissions were reduced with increased harvest frequency, these emissions do not consider off-site decomposition of harvested biomass.The harvest method implemented in the PVN model was similar to the instantaneous harvest method featured in other dynamic vegetation models (such as JULES, Littleton et al. (2020)), where the plant is reduced to a certain set height, and living biomass and LAI are subsequently adjusted accordingly.JULES assumes that 100 % of lost biomass is harvested whilst killing off a proportion of belowground biomass that is converted to litter.The PVN model assumes that 20 % of harvested biomass is lost to litter and does not account for root death.The increased litter layer leads to enhanced emissions resulting from the decomposition of the litter layer.The PFT living biomass is reduced by the proportional biomass lost, assuming the plant's biomass is uniformly distributed with height, and LAI is recalculated.Root mass observational measurements over time, as well as observational data on the impact of harvests on plant productivity, would further improve model representations of harvests.A further assessment may investigate in what ways the photosynthesising and gas conduit capacities of plants are further reduced in the days after harvest and how this can be better captured by the model.

Comparison to other site-specific peatland GHG emission models
Here, we compare the functionality of the new PVN model against its parent models, the Peatland-VU and NUCOM-BOG models.We then also compare the functionality of the PVN model against the functionality of PEATBOG, the one other site-specific peatland GHG emissions model that includes dynamic vegetation (Table S1).We have developed a new model capable of understanding the role dynamic PFTs play in CO 2 and CH 4 emissions in peatlands.The PVN model simulation results estimated the 1995-2017 annually averaged net GHG budget to be larger than that of the Peatland-VU model at both sites.We suspect that there are two reasons for this, the first being a trade-off between enhanced CO 2 emissions or enhanced CH 4 emissions.In both the Peatland-VU and PVN models, the CO 2 processes are calculated first.Calibration of the photosynthesis-and plantrespiration-related parameters impacts the amount of CO 2 available for CH 4 production.Photosynthesis and leaf respiration mechanisms were the greatest cause of uncertainty in the model's ability to reproduce the net GHG budget.Future model versions may consider ways to constrain the net CO 2 flux by improving the response of photosynthesis to environmental variables.To improve upon this in future model versions, it may be useful to consider the representation of belowground carbon decomposition.The belowground CH 4 pool in the Peatland-VU model increased consistently during the model simulation; therefore, an increasing quantity of CH 4 was released from the soil profile throughout the simulation, indicating that these fluxes were likely underestimated early in the simulation.The PVN model prescribes each PFT to have root and shoot mass and root depths.This enables each PFT to access different soil layers and belowground CH 4 and carbon pools, potentially impacting the long-term variability of CH 4 emissions.When compared to observed fluxes, the results indicated that the CO 2 scheme in the PVN model may have limited skill when applied to peatland sites of certain physical properties.These results cannot be compared with previous modelling studies because the Peatland-VU CO 2 production scheme results have not been published since the CO 2 production scheme was introduced by Mi et al. (2014) for assessment of the impact on simulated CH 4 fluxes.
The NUCOM model was developed to assess the impact of climate change on bog ecosystems by analysing simulations lasting 200-500 years.Running the model over time periods similar to the NUCOM's 1760-2000 simulation period can be used to assess the model's ability to reproduce shifts in vegetation in response to climate variability.This would require model evaluation using multi-centennial observational data, such as macrofossil evidence.To further investigate the impact of climate change on peatland ecosystems, future studies may consider using macrofossil data in combination with forward or backward multi-decadal or multi-centennial climate projections.
The PEATBOG model (Wu and Blodau, 2013) is the one other site-specific peatland model that simulates CO 2 and CH 4 fluxes and includes competition between moss, shrub, and graminoid PFTs.The PEATBOG model has simulated the Mer Bleue bog in Canada, a pristine (untouched) raised acidic ombrotrophic bog, over a 6-year period.The Mer Bleue bog is a nutrient-poor bog unlike the two sites assessed in this study.The net annual GHG emissions for the Mer Bleue bog site were small, approximately 0.02 % of the GHG emissions observed at the Ilperveld field site in the Netherlands.Peat has been accumulating at this site since 8400 cal yr BP and has developed a peat depth of 6 m in the centre.The PEATBOG model is a complex model that simulates many of the same processes as the PVN model but beyond this also includes representation of the nitrogen cycling, electron-accepting processes, dissolved inorganic and organic carbon, and subsequent CO 2 and CH 4 run-off.The PEATBOG model underestimated the annual net GHG emissions (net ecosystem carbon balance) by approximately half of observed field observations.Wu and Blodau (2013) noted the sensitivity of the PEATBOG model to temperature, re-porting that 1 • C of temperature change was enough to initiate a model bias, swaying the model from a source to a sink.This is concurrent with the results of the sensitivity testing performed in this study, which showed that changes in air temperature had large impacts on both CO 2 and CH 4 emissions.Plot-scale model inter-comparison efforts could help improve the representation of small-scale processes in peatland models.However, the breadth of observational data required to run and test site-specific models makes site-specific model inter-comparison efforts cumbersome and difficult.

Conclusions
Here, we present Peatland-VU-NUCOM v1.0 (PVN), a new site-specific peatland dynamic vegetation emissions model.By including plant-environmental feedbacks, the model can serve wetland management by estimating changes in the GHG balance of peatland sites in response to environmental changes -such as changing air temperatures, water levels, or precipitation and/or evapotranspiration -or new management decisions -such as raising the water table, vegetation restoration, or modifying mowing regimes.The PVN model was designed to simulate plant competition above-and belowground whilst developing carbon pools for the production and oxidation of CH 4 and CO 2 .PFTs compete for light where production and respiration are dependent on ideal temperature and water levels.Structural differences in vegetation root, exudation, and stem representation impact CH 4 production, oxidation, and transport.Peatlands are one of the most important carbon-storing ecosystems.The challenges facing our understanding of the carbon balance and CH 4 dynamics subsequent to the rewetting of previously managed peatlands are numerous.One challenge is the ability of sitespecific peatland models to reproduce CH 4 fluxes, particularly in relation to plant functioning.This question is particularly timely because there exists an urgent need to restore drained peatlands to reduce land subsidence whilst limiting GHG emissions.We show that the PVN model is able to reproduce plant biomass fractions and CH 4 and CO 2 fluxes under changing environmental conditions.This confirms that the model provides the capability to understand the relationship between peatland plant dynamics and CH 4 and CO 2 emissions.The PVN model is a relevant tool that can be used to optimise vegetation management with the goal of reducing GHG emissions.

Figure 2 .
Figure 2. The results of the sensitivity tests show the relationship between different temperature inputs and the mean daily plant-transported CH 4 for each year for each of the PFTs at the Horstermeer site (top row) and Ilperveld site (bottom row).Temperature input was increased and decreased by 1 and 3 • C, respectively.The legend shows the input change [ • C] where ± signs in front of the legend labels show the direction of change.Note the different y axes between the top and bottom panels.

Figure 3 .
Figure 3.The results of the sensitivity tests show the relationship between different water level inputs and the mean daily soil CO 2 flux for each year for each of the PFTs at the Horstermeer site (top row) and the Ilperveld site (bottom row).Water level input was decreased by 0.1 and 0.2 m and was increased by 0.1 and 0.2 m, respectively.The legend shows the input change, where ± signs in front of the legend labels indicate the direction of the change.Note the different y axes between the top and bottom panels.

Figure 4 .
Figure 4.The results of the sensitivity tests show the relationship between different harvest schemes and biomass for each day of year (shown as a fraction of litter and root mass) at the Horstermeer site (top row) and Ilperveld site (bottom row).Vegetation is cut to 0.15 m (×0.15) at the moment of harvest.The legend shows the harvest input scheme, and the vertical dotted lines indicate the four possible harvest days (days 120, 186, 220, and 268).Harvest was set to either not occur (H0.0) or to occur once per year (H1x0.15)on day 268, twice per year (H2x0.15)on days 186 and 268, three times per year (H3x0.15)on days 120, 220 and 268, or four times per year (H4x0.15)on all harvest days.

Figure 5 .
Figure 5.The CH 4 and CO 2 emissions for various isolated model mechanisms compared against the PVN model result.We investigated maintaining constant fractional PAR (PVN_FPAR_CONST), maintaining constant plant height (PVN_HEIGHT_CONST), maintaining constant cover fraction (PVN_CF_CONST), and including the original Peatland-VU CH 4 module multiplied by the PFT cover fraction (PVN_CH4_OLD_CF) at each time step.

Figure 6 .
Figure 6.Relative contributions of each PFT to simulated annual average net GHG (a, d), CH 4 (b, e), and CO 2 (c, f) emissions.The results of the Horstermeer site simulation are represented in (a)-(c), and the results of the Ilperveld site simulation are represented in (d)-(f).

Figure 7 .
Figure 7. Vegetation dynamics.The results of the Horstermeer site simulation are represented in (a), (c), and (e), and the results of the Ilperveld site simulation are represented in (b), (d), and (f).Note the differing y axes.

Figure 9 .
Figure 9. Simulated and observed methane fluxes at the Horstermeer (a, c, e) and Ilperveld (b, d, f) sites.The R 2 values are provided for comparison between the new PVN, the Peatland-VU model, and the observations.In (a) and (b), the 1 : 1 line is plotted in grey.The black dots are in situ flux chamber observational measurements in (c), (d), (e), and (f).Note the differing x and y axes.

Figure 10 .
Figure 10.Simulated and observed carbon dioxide fluxes (NEE) at the Horstermeer (a, c, e) and Ilperveld (b, d, f) sites.The R 2 values are provided for comparison between the new PVN, the Peatland-VU model, and the observations.In (a) and (b), the 1 : 1 line is plotted in grey.The black dots are in situ flux chamber observational measurements in (c), (d), (e), and (f).Note the differing x and y axes.
3 and 1.2 kg CO 2 eq.m −2 yr −1 for the 2015-2017 and 1995-2017 periods, respectively.The comparison of modelled and measured CH 4 emissions showed that the PVN model performed well, reproducing CH 4 emissions within the spread of observations, in comparison to the Peatland-VU model.The PVN Horstermeer simulation results estimated large mean annual CH 4 emissions (5.1 kg CO 2 eq.m −2 yr −1 ) in comparison to the Peatland-VU model (3.2 kg CO 2 eq.m −2 yr −1 ) for the period 2015-2017.The R 2 value of the PVN model results in comparison to the observations was 0.7 for the Horstermeer simulation and 0.8 for the Ilperveld simulation.In comparison, the Peatland-VU model results produced R 2 values of 0.3 and 0.6 for the Horstermeer and Ilperveld simulations, respectively.The Peatland-VU model showed good skill in reproducing the CO 2 fluxes at the Horstermeer site (R 2 = 0.7) and less skill at the Ilperveld site (R 2 = 0.6).Similarly, the PVN model showed good skill in reproducing daily CO 2 fluxes at the Horstermeer site (R 2 = 0.8) but less skill at the Ilperveld site (R 2 = 0.6), as indicated by the linear regression results.Overall, assessment of the linear regression results showed that the behaviour of the PVN model performed well against the observations when compared to the Peatland-VU model.
Code and data availability.All model code has been written in C++.The model code is publicly available from the Bitbucket repository (https://www.bitbucket.org/tlippmann/pvn_public,last access: 15 November 2023) under the GNU General Public License version 3 or any later version.Users are welcome to contact the authors for technical support.All input data used to generate the model simulations presented in this study can be accessed through this Bitbucket repository.This includes site model paramhttps://doi.org/10.5194/gmd-16-6773-2023Geosci.Model Dev., 16, 6773-6804, 2023

Table 1 .
Name, units, description, and values of PFT input parameters.Associated references are listed in TableS3.In the left column, each PFT parameter is tied to its relevant model mechanism.Note that some PFT parameters are, at times, used by multiple model processes.
Poaceae and are grass-like plants with elongated, long, blade-like leaves without aerenchyma.Typha PFTs represent a genus of about 30 species of monocotyledonous flowering plants in the family Typhacea with large aerenchyma.The short grasses PFT is representative of forbs and agriculturallike grasses with shallow root systems.The Sphagnum PFT is representative of hummock Sphagnum species which are generally more drought tolerant.Brown mosses represent all non-Sphagnum mosses but have similar but slightly broader temperature growth ranges.The SOM evolved from short grasses decomposes more easily than the SOM evolved from brown mosses, which decomposes more easily than the SOM evolved from Sphagnum.The six PFT input parameter sets used in this study are accessible from the Bitbucket repository: https://www.bitbucket.org/tlippmann/pvn_public(last access: 15 November 2023).
. Sedges, tall grasses, and Typha all represent graminoids with deep root systems that can grow at a range of water levels but have different aerenchyma and growing ranges.Sedges are from the families Cyperaceae and Juncaceae and are grass-like, monocotyledonous flowering plants with aerenchymae.Tall grasses are from the famhttps://doi.org/10.5194/gmd-16-6773-2023Geosci.Model Dev., 16, 6773-6804, 2023 ily

Table 2 .
A summary of the varied input data used to understand the sensitivity of the model.To compare the PFT dynamics, both simulations use the no-harvest regime.The exchange of PFTs means that the model simulation driven by the Ilperveld input data (Table

Table 3 .
A summary of the simulations with altered model algorithms.VU CH 4 module multiplied by the PFT cover fraction CF_CONST Biomass fraction is constant for all PFTs, i.e.BF = 0.25 FPAR_CONST FPAR is constant for all PFTs, i.e.FPAR = 0.25 HEIGHT_CONST Constant plant height for each PFT simulated by the Peatland-VU model to understand the impact of introducing PFTs on the simulation of CH 4 and CO 2 fluxes.These model simulations are summarised in Table over the same period.The average daily temperature for the warmest month, August, was 22.1 • C, and the lowest daily monthly temperature for the coldest month, January, was 0.8 • C.

Table 4 .
A summary of the model simulations using both the new PVN model and the pre-existing Peatland-VU (PV) model.Model input parameters for the Horstermeer and Ilperveld site simulations are provided in TablesS6 and S7, respectively.

Table 5 .
The results of the sensitivity testing.The CH 4 and CO 2 columns indicate how much the respective emissions changed when the input changed relative to the results of the respective default Horstermeer and Ilperveld PVN simulations described in Table4.A dash[-]indicates that the simulation is the default site simulation.An overview of the sensitivity tests can be found in Table2.
2003rier-Uijl e(van Huissteden et al., 2009))ock grazing(Schrier-Uijl et al., 2014).It remains unclear what impacts these events continue to have on present-day CO 2 and CH 4 fluxes.Unfortunately, at the time of publication, there were no published studies investigating the CO 2 or CH 4 fluxes measured at the Ilperveld site.The CH 4 fluxes observations (0-17 mg CH 4 m −2 h −1 ) presented in this study compared well to reported chamber CH 4 fluxes measured at the Horstermeer site from2003  till 2008(van Huissteden et al., 2009)), in the range of 2-15 mg CH 4 m −2 h −1 at an area of the site with a varying water table.Interestingly, the CH 4 observations presented measured in a wet area of the Horstermeer site were more than double the measurements measured in dry areas of the Horstermeer between 2004 and 2006 using the manual chamber technique eterisations, site soil profiles, climate data, water level data, and PFTs.The exact version of the model source code used to produce the results presented in this paper is archived on Zenodo (https://doi.org/10.5281/zenodo.7701698,LippmannandvanHuissteden(2023)).Supplement.The supplement related to this article is available online at: https://doi.org/10.5194/gmd-16-6773-2023-supplement.Author contributions.TJRL, KvH, and MMPDH developed the theoretical framework of the model.TJRL performed the model developments, composed the PFTs, made the figures, conducted all analyses, and wrote the paper.TJRL and KvH collected the observational data and developed the model parameterisation scheme and soil profiles.KvH processed the observational data and offered valuable suggestions on the development and calibration of the model.YvdV offered valuable suggestions on the testing of the model, the presentation of analyses, and the writing of the paper.HD participated in the writing of this paper.HD, MMPDH, KvH, and YvdV participated in the revision of this paper.DMDH, HD, MMPDH, and KvH acquired the funding and administered this project.Competing interests.The contact author has declared that none of the authors has any competing interests.Disclaimer.Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper.While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.