Identifying main uncertainties in estimating past and present radiative forcing of peatlands

Reconstructions of past climate impact, that is, radiative forcing (RF), of peatland carbon (C) dynamics show that immediately after peatland initiation the climate warming effect of CH4 emissions exceeds the cooling effect of CO2 uptake, but thereafter the net effect of most peatlands will move toward cooling, when RF switches from positive to negative. Reconstructing peatland C dynamics necessarily involves uncertainties related to basic assumptions on past CO2 flux, CH4 emission and peatland expansion. We investigated the effect of these uncertainties on the RF of three peatlands, using either apparent C accumulation rates, net C balance (NCB) or NCB plus C loss during fires as basis for CO2 uptake estimate; applying a plausible range for CH4 emission; and assuming linearly interpolated expansion between basal dates or comparatively early or late expansion. When we factored that some C would only be stored temporarily (NCB and NCB+fire), the estimated past cooling effect of CO2 uptake increased, but the present‐day RF was affected little. Altering the assumptions behind the reconstructed CO2 flux or expansion patterns caused the RF to peak earlier and advanced the switch from positive to negative RF by several thousand years. Compared with NCB, including fires had only small additional effect on RF lasting less than 1000 year. The largest uncertainty in reconstructing peatland RF was associated with CH4 emissions. As shown by the consistently positive RF modelled for one site, and in some cases for the other two, peatlands with high CH4 emissions and low C accumulation rates may have remained climate warming agents since their initiation. Although uncertainties in present‐day RF were mainly due to the assumed CH4 emission rates, the uncertainty in lateral expansion still had a significant effect on the present‐day RF, highlighting the importance to consider uncertainties in the past peatland C balance in RF reconstructions.


| INTRODUC TI ON
Peatlands constitute an effective sink of atmospheric carbon dioxide (CO 2 ), but they are also an important source of methane (CH 4 ; Korhola et al., 2010;Petrescu et al., 2015;Yu, 2011) and, under some circumstances, nitrous oxide (N 2 O; Freeman et al., 1993;Langeveld et al., 1997;Martikainen et al., 1993;Repo et al., 2009). The exchange of these greenhouse gases (GHGs) between peatlands and the atmosphere has a dualistic effect on the climate, as one GHG (CO 2 ) is removed from the atmosphere while others are emitted.
The current yearly CO 2 uptake by northern high latitude peatlands amounts to 40-66 Tg (Tg = 10 12 g) of carbon (C; Turunen et al., 2002;Yu, 2011). During the Holocene (the last ca. 11,700 year), northern peatlands have accumulated over 500 Pg C (Pg = 10 15 g), which is equivalent to the pre-industrial atmospheric C reservoir.
Simultaneously, CH 4 losses to the atmosphere account for a significant proportion, up to 25%, of the net ecosystem C balance of peatlands (Limpens et al., 2008), adding up to ca. 15 Tg CH 4 -C year −1 (Mikaloff-Fletcher et al., 2004). The magnitude of the CO 2 sink and CH 4 source has varied throughout the Holocene (Yu, 2011). This is due to the peatland GHG fluxes being regulated by not only autogenic peatland succession (Juottonen et al., 2021;Leppälä, Laine, et al., 2011;Leppälä, Oksanen, et al., 2011) but also to a great extent by climatic factors such as temperature, humidity, insolation, and precipitation (Dorrepaal et al., 2009;Fan et al., 2013;Gorham, 1991). Hence, to understand the role and feedback mechanisms of peatlands in the past global climate system, we need to assess their GHG exchange over the Holocene.
The climatic effect of peatland GHG fluxes can be expressed as radiative forcing (RF; Frolking et al., 2006), which is defined as the change in the global radiative balance (energy flux per area of Earth's surface, W m −2 ). GHG emissions increase atmospheric GHG concentrations reducing the energy radiating to space; hence, this is defined as a positive RF that results in climate warming (Myhre et al., 2013).
Removal of GHGs from the atmosphere has the opposite effect, that is, a negative RF resulting in climate cooling. RF depends not only on the magnitude of the GHG flux in question but also substantially on the radiative efficiency and time scales related to biogeochemical cycling of the GHG. This is especially relevant when comparing CO 2 uptake, which is associated with sustained forcing even when occurring in the short term, and a CH 4 emission pulse resulting in only a transient effect (Frolking et al., 2006).
The majority of northern peatlands initiated during the early Holocene, and their uptake and storage of CO 2 since then have affected atmospheric CO 2 concentrations (Yu, 2011). Several studies have reconstructed northern peatland GHG fluxes since the start of the Holocene (Gorham, 1991;Yu, 2011Yu, , 2012Yu et al., 2010), but few have considered the RF associated with these fluxes (Frolking & Roulet, 2007;Mathijssen et al., 2014Mathijssen et al., , 2017Piilo et al., 2020). Frolking and Roulet (2007) estimated the RF associated with the development of northern peatlands throughout the Holocene, assuming linear peat expansion after initiation dates and various scenarios of CH 4 emissions, resulting in a current RF of −0.22 to −0.56 W m −2 .
In general, peatland RF will turn negative eventually due to sustained CO 2 uptake, overcoming the positive RF of CH 4 emissions.
Uncertainties in reconstructions of peatland C exchange can be divided into uncertainties pertaining to the net CO 2 exchange (as the product of CO 2 uptake and release) and CH 4 emissions. CO 2 uptake reconstructions are commonly deduced from the amount of C stored in peat layers of varying ages, assuming that the organic material that was formed when a certain layer was at the surface stays in place after that layer is buried under younger peat layers. It is common practice not to directly measure the C content throughout a peat profile. Instead, the proportion of organic matter is analysed, which is then multiplied with an assumed C content of organic matter (Loisel et al., 2014). However, the organic matter C content varies from 42% to 57% depending on peat type and decomposition state (Beilman et al., 2009;Loisel et al., 2014).
Estimation of all these C fluxes, that is, CO 2 uptake, CO 2 release and CH 4 emissions, is additionally affected by uncertainties in the areal development of a peatland, which is commonly reconstructed by interpolating areal growth between peat initiation and the present-day peatland size (Korhola, 1994). Variability in expansion rate is then based on the age and distribution of multiple basal peat ages . This approach relies on the accurate localization of the oldest part of the peatland and a sufficiently large number of basal ages evenly distributed across the present-day peat area. Furthermore, any location where the peat layer has disappeared at any point in the past cannot be considered.
In paleoecological peatland studies, the long-term C fluxes are commonly summarized as apparent C accumulation rates (aCAR; g C m −2 year −1 ; Loisel et al., 2014), calculated as the cumulative C stored between two peat layers of known age. However, the actual peat C balance at any given time is the result of both C uptake by vegetation at the peat surface, and C loss from all the peat layers that existed at that time. Even though the actual peat C balance were temporarily negative, that is, C loss from all layers is larger than uptake at the surface, the aCAR would always remain positive and average out any variability of the C balance over time (Frolking et al., 2014), in addition to overestimating recent C uptake rates (Young et al., 2019). However, a paleo-reconstruction of the peatland net C balance (NCB) can be modelled using aCAR as the present-day end result of past C uptake and subsequent decomposition with a constant loss rate (Yu, 2011). Additional uncertainties in C balance reconstructions stem from disturbances, such as peat fires (Loisel et al., 2017;Turetsky & Wieder, 2001), and loss of dissolved organic C through lateral water flow (Evans et al., 2016).
One of the main uncertainties in reconstructions of past peatland C dynamics is the estimation of past CH 4 emissions (Loisel et al., 2017). However, as a consequence of its relationship with vegetation characteristics, it is possible to reconstruct CH 4 emissions by using fossil plant species assemblages as proxy , similar to the use of vegetation as an indicator for current CH 4 fluxes (Bubier et al., 1995;Couwenberg et al., 2011;Gray et al., 2013).
However, the uncertainty in such CH 4 flux predictions often exceeds 50% (Bubier et al., 1995;Mathijssen et al., 2016). Another way to estimate peatlands' past CH 4 emissions is to determine their past type and trophic state and assume a typical flux rate as observed in similar peatlands at present (Mathijssen et al., 2014(Mathijssen et al., , 2017Piilo et al., 2020). The main sources of uncertainty here would be the variability in observed flux rates and their temporal representativeness at the study site.
Several data compilations have been undertaken to synthesize global or continental-scale reconstructions of peatland development and their effect on the C cycle (Kleinen et al., 2012;Korhola et al., 2010;Yu, 2011). However, the Holocene reconstructions of northern peatlands and their climate effects have been performed based on either limited data of peatland C dynamics and lateral expansion patterns (Wang et al., 2009) or simulated data (Frolking & Roulet, 2007). Thus, before upscaling site-scale reconstructions, it would be prudent to investigate the uncertainties involved in them.
This study aims to quantify the relative uncertainties in RF pertaining to the varying assumptions necessary in reconstructions of peatland C dynamics and to assess how these assumptions relate to the interpretation of peatland-climate feedbacks. For that, we have selected three sites for which RF has previously been calculated throughout their entire development based on long-term estimates of C uptake, CH 4 emission and the lateral growth of peat surface.

| Approach
In this study, we focused on three Finnish peatlands for which the C dynamics and RF have recently been reconstructed throughout the Holocene (Mathijssen, 2016;Mathijssen et al., 2014Mathijssen et al., , 2016Mathijssen et al., , 2017. While the sites are characterized by relatively low overall C accumulation rates (Mathijssen et al., 2014(Mathijssen et al., , 2017, they are representative to the variation in northern peatlands in their vegetation and ecohydrology (Rydin & Jeglum, 2013). For our purposes they provide a useful sample because data collection methods and the assumptions inherent in their reconstructions of C dynamics were comparable. Furthermore, these three sites are some of the few for which Holocene-scale C and RF dynamics have been reconstructed (but see Dommain et al., 2018;Piilo et al., 2020).
To address our aim, we limit our study to the empirical data available for our sites; that is, we did not test the effect of having more or fewer measurement data. Consequently, we do not aim at providing 'best practice' guidelines on C dynamics measurements for RF reconstructions that would be better fitted to a modelling study.
The calculation of RF over the development history of each peatland was based on the reconstructed CO 2 uptake and CH 4 emission rates and peat area development, that is, lateral expansion, which were spatially integrated into annual peatland-scale CO 2 and CH 4 fluxes. In this paper, we refer to the CO 2 uptake and CH 4 emission rates expressed per unit area as 'flux densities' (g C m −2 year −1 ) and use the term 'flux' (g C year −1 ) for the rate of peatland-scale C exchange. In the original, site-specific studies, the reconstructed CO 2 uptake was assumed to equal aCAR, and the estimates of reconstructed CH 4 emission and lateral expansion were based on vegetation composition and linear interpolation between basal dates, respectively (Mathijssen et al., 2014(Mathijssen et al., , 2017. However, the uncertainties of these reconstructions were not translated into uncertainties of RF. To understand how the uncertainties involved in the reconstructions of CO 2 and CH 4 fluxes and lateral expansion affect the resulting RF, here the RF from peat initiation to 0 cal. year BP (before present, i.e., 1950 AD) was recalculated for these peatlands using varying approaches. In each approach, the empirical data underlying the original reconstructions were kept constant; that is, C-contents, age-depth models, vegetation composition, number, and age of basal peat samples were the same as in the original studies (Mathijssen et al., 2022). For all three variables (i.e., CO 2 and CH 4 fluxes and lateral expansion), we implemented three different reconstruction methods: the reconstructed CO 2 fluxes were based either on aCAR, net carbon balance (NCB) or NCB plus the C loss due to fires; CH 4 fluxes had 'minimum', 'average,' and 'maximum' emission variants; peat area was assumed to have had an 'early', 'expected' or 'delayed' expansion compared with the best estimate derived from observed basal ages. We followed the micrometeorological sign convention for CO 2 and CH 4 fluxes: a positive sign indicates a flux from the ecosystem to the atmosphere (emission) and a negative sign a flux from the atmosphere to the ecosystem (uptake). 14 ha has a maximum depth of 2.5 m. Peat accumulation started at ca. 10 kyr BP, and the peatland vegetation has resembled a minerotrophic fen throughout its development (Mathijssen et al., 2014).

| Study sites
Additional information on the contemporary C exchange at LOM can be found in Aurela et al. (2009Aurela et al. ( , 2015 and Zhang et al. (2020).

| Carbon dioxide flux
In all our CO 2 flux reconstruction approaches, the CO 2 that was taken up and subsequently lost by the peatland in the form of CH 4 emissions or lateral flow was not considered to be contributing to the atmospheric CO 2 store, that is, we assumed that the C in CH 4 emissions and dissolved in lateral flow originated mainly from the recently assimilated CO 2 (Cooper et al., 2017) and that both rapidly return to the atmosphere as CO 2 (Evans et al., 2016;Frolking et al., 2006). All the approaches used for CO 2 flux reconstruction were constrained in their present-day total accumulated C uptake by the cumulative C observed in peat cores at the time of sampling, and thus the present-day C stocks of the sites did not differ between these approaches.

| Apparent carbon accumulation rate
Peat cores down to the peat base were radiocarbon dated and analysed for C content and bulk density (for details, see Mathijssen et al., 2014Mathijssen et al., , 2016Mathijssen et al., , 2017. Based on these data, the aCAR (g C m −2 year −1 ) was calculated using where PAR (m year −1 ) is the peat accumulation rate resulting from agedepth models of the radiocarbon ages, ρ (g m −3 ) is peat bulk density, and c (g C g −1 ) is bulk peat C content.
The calculation of aCAR was performed for two peat cores from SII, eight cores from KAL, and one core from LOM. In the case of LOM, the aCAR flux densities (g C m −2 year −1 ) were multiplied by the reconstructed peatland area (see below) leading to peatland C fluxes (g C year −1 ) that were converted to CO 2 flux (g CO 2 year −1 ; Mathijssen et al., 2014). In the cases of SII and KAL, with multiple aCAR records from various locations within the site, the aCAR flux densities were averaged and multiplied by peat area, taking into account the different successional stages across the site (Mathijssen et al., 2017); that is, the aCAR values from the bog-vegetation stage at SII were applied only to the bog area present at that time, while the flux densities from the other core from SII containing fen-vegetation were applied to the fen area of SII at that time ).

| Net carbon balance
We applied the NCB model (Yu, 2011) to quantify the C uptake and total C loss through decomposition at any given time during peatland development. This approach was adapted for the use in this study by dividing the aCAR flux densities into 1000-year binned time intervals, multiplying these by the reconstructed peat area during the respective time interval, and using the products as a net peat C pool (NCP, Gg C kyr −1 ). The NCP and a decay coefficient (α) of 0.0001 year −1 (Clymo et al., 1998;Rydin & Jeglum, 2013) were then used to calculate the net C uptake (NCU) and net C release (NCR) following: in which NCU represents the initial C input into the catotelm and NCR represents the summed C release due to decomposition in all peat layers (k) present at time t. The NCB was then calculated as

| NCB combined with fire loss
The third approach to model past CO 2 fluxes assumed that C from organic matter buried in peat would not only be lost due to ongoing decomposition, but also by combustion in peat fires. The presence of charcoaled macrofossil remains was used as an indicator for past fires ( Table 1; Mathijssen et al., 2016Mathijssen et al., , 2017. This approach followed the NCB method with the following modification (NCB-F): time intervals during which a fire took place experienced a C loss of 2.0 kg (3) have been accumulated in the previous time interval. This meant that the burned and lost C was stored in the peatland for 1000 year.

| Methane emission
Our reconstructions of CH 4 emissions were based on contemporary CH 4 fluxes, which were then applied to the past by linking them to the reconstructed vegetation composition. The estimated CH 4 flux densities were also averaged over 1000-year intervals and multiplied by the corresponding reconstructed peatland areas. For SII, a weighted averaging transfer function was used to estimate the past emissions using plant macrofossils as input and chamber flux measurements with corresponding vegetation composition from SII and an additional site as a training set . For SII, the 'average' CH 4 emission variant equalled the predicted emissions, while the 'minimum' and 'maximum' emission variants were defined as the predicted emissions minus and plus the sample-specific errors, respectively. For KAL, the peatland CH 4 flux was calculated from flux densities of subareas for which the development of peat type was described separately (Mathijssen et al., 2017). Each subarea was assumed to have a CH 4 flux density corresponding to peat type, following Minkkinen and Ojanen (2013), who collected contemporary CH 4 flux data of Finnish peatlands. The 'average' CH 4 emission variant used the mean contemporary flux density, while the 'minimum' and 'maximum' emission variants were defined as the mean minus and plus standard deviation, respectively. For LOM, site-specific contemporary CH 4 flux measurements were available (Mathijssen et al., 2014) and adopted for the 'average' CH 4 emission variant. The 'minimum' and 'maximum' emission variants were defined as the mean measured flux density in LOM minus and plus the standard deviation of the fluxes at Finnish-rich fens (Minkkinen & Ojanen, 2013), respectively.

| Peatland area
Reconstructions of peatland lateral expansion were based on the age of basal peat samples spread across each site (18 at SII, 19 at KAL and 9 at LOM). The 'expected' variant of peat area development consisted of the best estimate of expansion based on linear interpolation between basal ages, considering base morphology.
This variant represented the reconstructions published previously (Mathijssen et al., 2014(Mathijssen et al., , 2017. The 'early' expansion variant assumed that lateral growth occurred as early as possible in the lifetime of the peatland while not contradicting the ages of the basal samples. Correspondingly, the 'late' variant assumed lateral growth to have occurred as late as possible.

| RF model
The peatland-scale CO 2 and CH 4 fluxes of the three sites were used to calculate the sites' RF during their development history. We used the sustained pulse-response model described by Lohila et al. (2010), Mathijssen et al. (2017) and Piilo et al. (2020) to calculate the RF resulting from the changes in atmospheric concentrations due to CO 2 and CH 4 exchange at the sites (RF CO2 and RF CH4 ). This model has previously been applied for these sites (Mathijssen, 2016;Mathijssen et al., 2014Mathijssen et al., , 2017, but here RF was recalculated because the model has been updated since these studies. This update included implementation of improved radiation efficiency functions, which increased the RF sensitivity to changes in atmospheric CH 4 concentration (Etminan et al., 2016). In addition, the updated version takes into account the variations in the background concentrations of the GHGs considered, that is, CO 2 , CH 4 and N 2 O (Köhler et al., 2017). For the purpose of this study, the fluxes were calculated in 1000-year intervals, and the yearly RF values modelled were averaged over these intervals. The RF model equations and parameter values are detailed in Supporting Information.

| Peat expansion
The reconstructed peatland expansion, interpolated from the distribution of basal dates, showed different patterns between the three sites. In SII, peat growth started in several small loci between 11 and 10 kyr BP, after which the peat area expanded rapidly until 5 kyr BP and slowly after that (Figure 1a,d). The uncertainty in the expansion pattern in SII was most pronounced during 10-5 kyr BP (Figure 1d).
In KAL, peat formation first initiated across an elongated region in the middle of the current peatland (Figure 1b), after which it expanded steadily over time. The relative uncertainty in the size of the KAL peatland area was highest at 8 kyr BP (Figure 1e). The expansion pattern of LOM contrasted with that of SII in that it contained a rapid peat establishment during 10-9 kyr BP, after which expansion slowed down, TA B L E 1 Fire occurrence based on the presence of macroscopic charcoal in peat samples No charcoal was observed in Lompolojänkkä.
b Large amounts of charcoal were observed indicating multiple fires or a large fire event.
to increase rapidly after 2 kyr BP (Figure 1c,f). The deviation of the 'early' and 'late' expansion patterns from the best estimate ('expected') amounted to a maximum peat area uncertainty at any given time of +20% ('early') to −25% ('late') for SII, +30% to −50% for KAL and +20% to −30% for LOM (Figure 1d-f). The present-day C stock varied according to expansion patterns and reached 692, 615 and 562 Gg C for SII, 50, 44 and 39 Gg C for KAL and 4.0, 3.5 and 2.9 Gg C for LOM after 'early', 'expected' and 'late' expansion scenario, respectively.

| Carbon dioxide uptake
Using the aCAR approach, the CO 2 flux densities of all three peatlands were relatively stable throughout their development  lateral expansion pattern (Figure 3). In SII and KAL, the uncertainty related to the choice of the CO 2 flux reconstruction method increased further after 4 kyr BP, relative to that due to the expansion pattern estimation (Figure 3a,b).

| Radiative forcing
Earlier peatland expansion caused RF CO2 and RF CH4 to decrease and increase faster, respectively ( Figure  during the most recent 1-kyr period ( Figure 6). The NCB-F approach produced RF CO2 trajectories that were very similar to those obtained using NCB (Figure 5a,b), with a maximum difference of 10%-15% during 9-5 kyr BP in KAL (Figure 5b,e).
In SII, the difference between the CH 4 flux reconstruction approaches was larger than that between the methods used to estimate CO 2 fluxes and lateral expansion (Figures 5d and 6). However, during 6-4 kyr BP, the difference between aCAR and NCB ap-

| Differences in original methods between sites
It is important to note, for the interpretation of the results of this study, that the exact methods used to reconstruct C fluxes differed among the three sites considered here. This is due to the fact that the original studies differed in their methods, as they attempted to make the most of the different types of data sets available for these peatlands. The studies of SII and KAL included data from multiple peat cores (Mathijssen et al., , 2017, and thus the reconstructed flux densities were based on multiple locations within the site and took into account the relative distribution of peat types and its changes throughout the Holocene. In contrast, the data from LOM only included one peat core to base flux estimation on (Mathijssen et al., 2014), so we had to assume that this core was representative of the entire peatland.
The largest difference in methods between the sites occurs in the reconstruction of CH 4 emissions. While for KAL and LOM we had to rely on collated measurement data (Minkkinen & Ojanen, 2013) and derive different scenarios from the variation in these data, with a difference between the scenarios ranging up to 3 g CH 4 -C m −2 year −1 (Figure 4b,c), in SII an effort was undertaken to model the past CH 4 emissions using macrofossil composition, which resulted in flux density uncertainties of over 20 g CH 4 -C m −2 year −1 (Figure 4a). The uncertainty in CH 4 flux density in KAL and LOM would be higher if we had taken into account that the mean fluxes from the collated data (Minkkinen & Ojanen, 2013) may not accurately represent these peatlands during their earlier development. As a result of this approach, SII had the largest uncertainties in CH 4 flux densities, even though the most detailed data were available from this site.
We retained these methodological differences between sites as our aim was to investigate how the uncertainties connected to reconstructions of past fluxes affect the reconstructed RF and not to compare the uncertainties among the three sites. However, we argue that each reconstruction approach can be considered a reasonable choice, and thus including various methods provides information on the range of uncertainties that should be considered when interpreting RF reconstructions and on how these uncertainties depend on the underlying site-specific data.

| Relative effect of uncertainties in CO 2 uptake, CH 4 emissions and lateral expansion
The uncertainty in the estimated CO 2 uptake rates had a strong impact on the modelled past RF, as indicated by the fact that during the early Holocene the NCB-based RF CO2 decreased faster than the RF CO2 calculated from aCAR data. Taking into account the temporary removal of CO 2 from the atmosphere increased the CO 2 uptake earlier in peatland development (before the subsequent C loss) and decreased them during the later stages ( Figure 2). However, while the order of CO 2 flux magnitudes in different trajectories changed during the mid-Holocene-with higher early-Holocene fluxes obtained with NCB and NCB-F than aCAR, while the opposite is true for the late Holocene due to ongoing C loss from deep and old peat layers-a similar change did not occur in RF CO2 (Figure 6). This switch and the cumulative nature of the RF due to sustained uptake of atmospheric CO 2 uptake (Myhre et al., 2013) resulted in a convergence of the RF CO2 of the different CO 2 uptake approaches during the late Holocene and in a similar present-day RF CO2 for the aCAR, NCB and NCB-F approaches (Figure 5a-c). This convergence was to be expected because all CO 2 flux approaches were constrained by the present-day C stock.

(b) (c) (a)
switchover occurred in these cases. In KAL and LOM, the CO 2 flux and lateral expansion scenarios had only a minor effect on the switchover, since its timing was predominantly dependent on a rapid decrease in CH 4 fluxes since 2 kyr BP (KAL; Figure 5b,e), or no switchover occurred at all (LOM; Figure 5f).
Although the effect of fires on CO 2 flux was substantial during the millennium in which the fire occurred (Figure 2), this effect was largely absent from the RF CO2 ( Figure 5). The effect of the C lost during a fire on RF was smaller than the effect of the C uptake increase generated into the reconstruction of the previous time interval. For example, the C uptake estimate for SII during 8-9 kyr BP was reduced by 11% when replacing NCB by NCB-F, but the amount of C lost in fires would appear to have been taken up during 9-10 kyr BP, enhancing the earlier RF CO2 and resulting in a RF CO2 at 8-9 kyr BP that is more negative by 3%. In this sense, including C losses by peat fires, but maintaining the same overall cumulative C pool, slightly enhances the estimated cooling effect during a peatland's lifetime because the uptake of CO 2 occurred earlier. In the long term, this effect becomes insignificant after a few thousand years since the last fire ( Figure 5). These findings are in line with those of Dommain et al. (2018), who found that fires would accelerate C losses from tropical peatlands but would not alter the conclusions about the development of RF.
Since a CH 4 emission pulse has a relatively short-term RF effect (Frolking et al., 2006), the instantaneous RF CH4 at any time is close to the equilibrium determined by the atmospheric perturbation time scale of CH 4 and the mean emission during a few decades before that time. This resulted in limited differences in RF CH4 between the different assumptions on peat area expansion ( Figure 5) because all expansion scenarios were constrained to reach the current peatland size at 0 kyr BP (Figure 1). In contrast, RF CO2 was cumulatively affected by the expansion scenarios. The 'early' scenario involved more CO 2 sequestered from the atmosphere over the peatland lifetime, which had the lasting effect of a more negative RF (Figures 5 and 6). In terms of peat area, the 'early' and 'late' expansion scenarios deviated from the 'expected' scenario by 15%-32% at 6-5 kyr BP (ca. midpoint of each peatland lifetime). The resulting RF CO2 and RF CH4 were affected by 21%-44% and 15%-33% at this time, but only by 9%-21% and 1%-13% at 1-0 kyr BP, respectively. The varying assumptions of peat expansion affected the modelled timing of the net RF peak and the switchover time; in SII, for example, the peak shifted by 2000 to 3000 year and the switchover time was varied by 1000 year (Figure 5d).
Comparing the relative importance of different types of uncertainties on net RF, there is a clear distinction to be made between considering the present-day and mid-Holocene net RF. The present-day uncertainties due to CH 4 flux densities had the largest effect, followed by the estimation of lateral expansion and lastly the CO 2 flux approach ( Figure 6). Only in KAL was the present-day uncertainty in RF CH4 reduced due to a recent transition from poor fen to bog (Mathijssen et al., 2017). The lower present-day relative uncertainties of expansion patterns and CO 2 flux can be explained by the fact that these are constrained by the present-day peat area and C stock, respectively, which are easily quantified. During the mid-Holocene, roughly at the midpoint of the study sites' development history, the CO 2 flux approach had the largest effect (NCB vs. aCAR) except in SII ( Figure 6). There the uncertainty due to CH 4 flux scenario was dominant throughout the Holocene (Figure 6a). During the mid-Holocene, the effects of CO 2 flux, CH 4 flux and expansion approaches were of similar magnitude in KAL and LOM (Figure 6b,c), although the largest effect was due to either CO 2 flux (KAL) or expansion (LOM). In SII, the effect of the expansion approach remained smaller than that of CO 2 flux until ca. 2 kyr BP (Figure 6a).
Furthermore, as RF CO2 and RF CH4 have opposite signs and both depend on the past trajectory of peatland area, in some cases the sign of net RF depends on the expansion approach adopted (Figure 5d,e) and thus the related uncertainty becomes insignificant at some point of time (Figure 6a,b).
The uncertainties stemming from the different approaches applied in this study to reconstruct the expansion and CO 2 and CH 4 fluxes of peatlands were larger than the uncertainties in the radiative efficiency parameterizations of our RF model, estimated at ca. 10% for CO 2 and 14% for CH 4 (Etminan et al., 2016), except in the case of the uncertainty in the present-day RF resulting from CO 2 flux approach. A previous RF reconstruction of a northern peatland showed similar results, where uncertain CO 2 fluxes during the early stages of peatland development had a lasting effect on the instantaneous RF induced much later (Piilo et al., 2020). Hence, these uncertainties should be taken into account when long-term RF reconstructions are interpreted.

| Impact on interpretation of peatland role in climate system
Assuming sustained and constant ecosystem-atmosphere exchange of CH 4 and CO 2 , the ratio between CH 4 emission and CO 2 uptake determines how long it takes before the switchover from positive to negative RF takes place (Frolking et al., 2006). When we calculated the molar CH 4 :CO 2 flux ratio (moles of CH 4 emitted per moles of CO 2 uptake) as the mean ratio between CH 4 and CO 2 fluxes from peat initiation until the switchover occurred, we found that with a CH 4 :CO 2 ratio of 0.2 the switchover time was less than 1000 year (SII with aCAR and min. CH 4 ; Figure 5d). This corresponds to the results presented for a hypothetical peatland by Frolking et al. (2006). The CH 4 :CO 2 ratios of 0.5 and 1.1 resulted in switchover times of ca. 6000 and 10,000 year, respectively (SII with NCB and average CH 4 and SII with aCAR and average CH 4 , respectively), while a CH 4 :CO 2 ratio larger than 1.5 would not result in a switchover within the whole peatland history, 11,000 year (SII with max. CH 4 ; LOM all cases). The switchover time of 6000 year is larger than that estimated by Frolking et al. (2006). In KAL, the switch from positive to negative RF was due to the large increase in CO 2 uptake during the last 1000 year (Figure 2b) rather than the cumulative effect of CO 2 uptake over time. Thus, in KAL the switchover time could not be related to the long-term CH 4 :CO 2 ratio.
The uncertainties in RF observed in this study were similar to the uncertainty range presented by Frolking and Roulet (2007). However, when standardized to peat surface area, the individual sites studied here reached a much higher present-day RF of −0.7 to +0.6 nW m −2 per hectare of peatland, than −1.9 to −0.6 nW m −2 ha −1 for northern peatlands collectively, where the range represents various scenarios of constant CO 2 uptake and CH 4 emission over the Holocene (Frolking & Roulet, 2007). Part of this difference can be explained by differences in RF modelling, especially by those affecting RF CH4 (Etminan et al., 2016). This is illustrated by the earlier RF calculations for KAL (Mathijssen et al., 2017) and LOM (Mathijssen, 2016), which resulted in present-day RF values that were ca. 0.4 nW m −2 ha −1 lower than in the current study. The remaining difference can be explained by the fact that SII, KAL and LOM had only slowly accumulated C over the Holocene and thus have lower average aCAR of 12.5, 8.9 and 6.3 g C m −2 year −1 , respectively, than the baseline flux densities of ca.
16 g C m −2 year −1 used by Frolking and Roulet (2007). Additionally, the CH 4 flux densities at our sites were higher than those (ca. 6 g CH 4 -C m −2 year −1 ) of Frolking and Roulet (2007). Our sites also exhibited low C accumulation rates compared with the mean northern peatland rate of 22.9 ± 2.0 g C m −2 year −1 (Loisel et al., 2014) and the mean Finnish subarctic fen rate of 16.8 g C m −2 year −1 (Turunen et al., 2002), and thus may not represent the average northern peatland. However, it should be noted that our RF reconstructions of individual peatlands, with relatively detailed data of C dynamics throughout their development, contain uncertainties of a similar magnitude to those involved in the diverse flux scenarios assumed by Frolking and Roulet (2007). This suggests that the uncertainties in large-scale RF reconstructions derived from peatland data syntheses will be significantly larger than in these simulations.
The persistently positive RF of LOM, resulting from a low C accumulation rate in combination with substantial CH 4 emissions, illustrates that even though a northern peatland has acted as a persistent sink of atmospheric CO 2 , it overall may have a climate warming effect. It must be noted that LOM is a fertile valley fen, in which a constant nutrient flow maintains the CH 4 emissions that are relatively high (Zhang et al., 2020) in comparison with net CO 2 uptake (Aurela et al., 2015). In addition, as shown by our results for SII and KAL, there have been extended warming periods during an early phase of peatland development. Even when taking into account the average fluxes of northern peatlands (Loisel et al., 2014), and their negative RF over the Holocene (Frolking & Roulet, 2007), short-term changes in peatland functioning can have large impacts on RF. For example, seen over a time window of 100 year, permafrost thaw in boreal wetlands is estimated to increase RF by 0.02-0.1 nW m −2 ha −1 (Helbig et al., 2017), and land-use changes may increase RF by as much as 1 nW m −2 ha −1 (Petrescu et al., 2015).
The results of the current study show that such disturbances could potentially counteract the entire negative RF generated by a 10,000 year old peatland.

| Impact on interpretation of peatlandclimate feedback
The present-day RF depends on the present-day fluxes of CH 4 and the cumulative effect of past CO 2 fluxes, but our results show that the method selected to reconstruct CO 2 flux densities (aCAR vs. NCB) has very little effect on the present-day RF, provided we can accurately estimate the present-day total C stock. The different RF trajectories and earlier uptake of CO 2 in the NCB compared with aCAR approach did not lead to diverging present-day RF values. While our present-day RF estimates showed substantial uncertainty also with respect to CH 4 flux estimation, in principle the present-day CH 4 emissions are fairly easy to constrain by flux measurements (Turetsky et al., 2014) and the total C stocks by assessments of peat coverage, depth and C content (Yu, 2012).
Hence, estimations of the peatland-scale present-day RF would mainly depend on an accurate reconstruction of lateral expansion. Contrary to the present-day RF, a reconstruction of peatland RF back in time, and the assessment of its function as peatlandclimate feedback in the past, depends heavily on the approach adopted to reconstruct C uptake rates, and to a lesser extent on the rate of lateral peatland expansion and the assumptions of CH 4 fluxes in the past (Figures 5 and 6).
It is important to note that the three sites studied here contained 9 to 19 basal dates that were used as a basis for reconstructing their lateral expansion. This many basal dates are typically not available for many peatland sites around the world. In a global reconstruction of peat expansion, Korhola et al. (2010) collected basal dates of over 2200 sites, of which 138 had more than three dates, and only 51 sites contained seven or more basal dates. The number of basal dates necessary to accurately reconstruct peatland expansion increases with peat initiation age, peatland size and complexity. This means that for many peatlands with fewer basal dates than were available for the three sites studied here, the uncertainty of lateral expansion rates would play an even larger role in assessing uncertainties in both the present-day and reconstructed RF. From our results, it seems reasonable to suggest that for sites with few basal dates the uncertainty in RF due to uncertain expansion patterns could be as large as 50%, which has not been taken into consideration in previous estimations of northern peatland Holocene RF (Frolking & Roulet, 2007;Wang et al., 2009).

| Unaccounted uncertainties
This study does not offer an exhaustive overview of the uncertainties involved in reconstructing the C-balance and related RF of peatlands but adopts the available measurement data from the study sites. If, instead, one would make the decision to return to the site and take more, or other type of, measurements, more sources of uncertainties could be tackled, which now lie outside the scope of this study. Among these would be the number and distribution of basal peat ages used to reconstruct peatland expansion. As the rate of lateral expansion varies within a site, depending on local factors such as topography (Korhola, 1994), the total number and distribution of basal ages necessary to gain an accurate expansion reconstruction would vary among peatland sites. It would be an interesting modelling exercise to test how the uncertainty in expansion pattern depends on the number and distribution of basal ages.
The number and location of the analysed peat cores will also have an effect on the reconstructed C uptake and CH 4 emissions, since C accumulation rates can vary significantly within a site (Korhola et al., 1995;Piilo et al., 2020), and so can water table depth and the associated vegetation and CH 4 emission rates (Turetsky et al., 2014).
However, peatland studies containing multiple cores that are analysed in detail are rare, since this is a very time-consuming task.
Furthermore, a reconstruction of past CH 4 emissions would ideally also take into account paleotemperatures, given that during warmer temperatures higher CH 4 emissions are expected irrespective of peatland vegetation type (Loisel et al., 2017). Lastly, in all three sites the fate of C which is lost from the peat through lateral flow is not taken into account, assuming that this returns to the atmosphere relatively quickly (Evans et al., 2016). However, if dissolved organic C leached from the peatland would be immobilized subsequently (McKnight et al., 2002), it could be considered to contribute to the long-term C accumulation of the peatland, further decreasing its RF.
The ages of peat layers are commonly analysed using radiocarbon dating and may have uncertainty ranges of up to ±200 year (Nilsson et al., 2001). Additional inaccuracy in peat ages could arise from contamination with 'old' or 'young' C in the analysed material, that is, uptake of CO 2 originating from decomposition of older material or the presence of roots belonging to vegetation growing at the location much later than the deposition of the studied layer, respectively. Varying 14 C fractionation among different species and plant tissue types could further decrease the precision of radiocarbon dates (Nilsson et al., 2001;Väliranta et al., 2014). However, these issues can be avoided by dating similar material across the studied samples, preferably Sphagnum remains (Nilsson et al., 2001). The remaining chronological uncertainty of a few hundred years should not have a major effect on RF reconstructions spanning multiple thousand years.

| CON CLUS ION
In this study we set out to investigate the effects on the modelled RF of differing approaches to reconstruct C uptake, methane emissions and lateral peat expansion throughout a peatland's lifetime. The early-to mid-lifetime RF of peatlands was heavily affected by these choices, except that adding the effect of C lost through peat fires did not result in a major effect. For most of the Holocene a peatland's estimated net RF could be either positive or negative, depending on the approach adopted to reconstruct the C fluxes and peat expansion. Using the NCB model instead of the concept of aCAR, or assuming early versus late peat expansion, could change the estimated timing of the switchover from positive to negative RF by several thousand years. Furthermore, an assumption of high past methane emissions can cause the peatland RF to never turn negative despite sustained CO 2 uptake during the Holocene.
Even for estimating the present-day RF, we need data for the whole development history due to sustained sequestration of atmospheric CO 2 . In this case, however, the uncertainties are mainly limited to the estimation of lateral expansion patterns. This is because cumulative C uptake and present-day methane emissions can be constrained by the present-day C stock and methane flux measurements at the sites, respectively. Hence, if one's aim is to estimate a peatland's present-day RF status only, it seems sufficient to base calculations on the present-day C stock, methane emission rate and peat initiation age, and obtain enough basal age estimates to accurately reconstruct the development of peatland coverage. However, if one is interested in the development of peatlands' RF over time since peatland initiation, more details are necessary concerning the chronology of peat C accumulation and probable past methane emission, and the inherent uncertainty in these factors should be taken into account in reconstructions of RF.

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.6241493. The detailed description of the RF model used in this study can be found in Supporting Information.