Blanket bog CO2 flux driven by plant functional type during summer drought

Recent climate predictions for the United Kingdom expect a nationwide shift towards drier and warmer summers, increasing the risk of more frequent and severe drought events. Such shifts in weather patterns impede functioning of global peatlands, especially rare intact blanket bogs abundant in Scotland and representing nearly a quarter of the UK's soil carbon. In this in situ study, carbon dioxide (CO2) fluxes from dominant peatland plant functional types (PFTs) such as Sphagnum spp., graminoids, ericoids and other key cover types (i.e., pools and bare peat) were measured and compared across upland and low‐lying blanket bog margins and centres, immediately before and during a summer drought in 2018, and over the subsequent year. During that period, most sites acted as net sources of CO2 to the atmosphere. Our results showed that net ecosystem exchange (NEE) was limited by water availability during the drought, with ericoid shrubs showing the highest drought resilience, followed by graminoids (which were still limited in GPP in 2019) and Sphagnum mosses. Diverging NEE estimates were observed across centre and margin areas of the blanket bogs, with highest variability across the upland site where signs of active erosion were visible. Overall, our study suggests that estimating growing season carbon fluxes from in situ peatland PFT and cover types can help us better understand global climate change impacts on the dynamics and trajectories of peatland C cycles.

Such changes in precipitation patterns can have direct (i.e., changing water availability and evapotranspiration rates) or indirect (i.e., altered surface water runoff or lateral drainage) effects on ecosystem distribution and functioning (Radu & Duval, 2018).
Quantifying these effects is needed to improve models that include long-term datasets or mean annual variable values to predict future climate, especially since these models often fail to acknowledge the extreme events . This is particularly relevant for ecosystems whose present bioclimatic space is likely to contract under future climate change scenarios, such as northern peat bogs (Chaudhary et al., 2020) and within the Scottish context, blanket bogs Gallego-Sala & Colin Prentice, 2013).
To improve our understanding of the consequences of extreme events, such as droughts, on carbon cycling in blanket bog, one approach is to characterise the individual responses of key plant functional types (PFTs). In blanket bogs, key PFTs consist of ericoid shrubs, graminoids (e.g., sedges), herbs, Sphagnum spp. and other mosses (Berger et al., 2018;Laine et al., 2016). Previous long-term studies using a PFT approach have demonstrated that peatland vascular PFTs (graminoids and shrubs) control gaseous methane emissions (Armstrong et al., 2015;Goud et al., 2021;Robroek et al., 2015), while increased Sphagnum moss coverage reduces fluvial dissolved organic carbon (DOC) effluxes (Dunn et al., 2016). Furthermore, peatland carbon cycling may be constrained by ericoid shrubs presence, as increased gross CO 2 fluxes from graminoids were observed during selective removal experiments of key PFTs, that is, ericoid shrubs (Ward et al., 2009). A climate-driven shift in PFT composition could therefore influence total C sequestration in UK blanket bogs. However, these manipulative studies do not capture the immediate, shortterm effects of actual drought-induced stress on greenhouse gas (GHG) emissions, which combine in situ reduced water availability and higher temperatures.
To date, the short-term effect of drought stress on C sequestration in peatlands has primarily been measured through mesocosm experiments (Dieleman et al., 2015;Fenner & Freeman, 2011;Kuiper et al., 2014). While useful to understand mechanisms and controls on PFT-response to drought, these laboratory-controlled settings may not fully reflect in situ responses, as they cannot account for local variations in microhabitat condition or topographic setting, which may influence GHG fluxes during periods of drought. For example, Sphagnum species located in higher elevation areas (i.e., mountains) suffer a stronger negative effect on CO 2 uptake during periods of increased temperature (Gerdol & Vicentini, 2011). Sphagnum spp. are one of the key peat-forming taxa in blanket bog, owing to their capacity to hold water and intrinsic properties (Bengtsson et al., 2016;Turetsky, 2003). However, Sphagnum spp. are known to be vulnerable to prolonged drought periods, as they are not well adapted to high temperatures and low precipitation rates (Bragazza, 2008;Gerdol & Vicentini, 2011;Van Breemen, 1995). Lowered water table depths in blanket bogs could cause a shift from Sphagnum spp. to other moss species that are better adapted to prolonged periods of droughts such as Polytrichum spp. (Potvin et al., 2014) and Racomitrium lanuginosum (Lindsay, 2010). Such change in vegetation would have potential consequence for long-term C uptake. Identifying the drivers and thresholds in atmospheric and edaphic conditions beyond which different PFTs shift from net CO 2 sinks to sources is therefore essential to improve our understanding of the potential consequences of droughts on blanket bog C dynamics.
In 2018, a persistent high-pressure system resulted in a period of drought across Europe between April and October (Buras et al., 2020;Hari et al., 2020). In Scotland, the months of May, June and July 2018 were nearly 2 C warmer than the long-term monthly mean  and received only between a quarter and a half of the longterm mean monthly (64 ± 15 mm compared to 86 ± 7 mm) rainfall (Supporting information, Figure S2).
This provided an opportunity to measure the immediate and short-term (1 year post drought) effects of the 2018 drought on CO 2 emissions from blanket bog sites in Scotland. For this study, we first aimed to compare in situ small-scale CO 2 fluxes from a range of PFTs (Sphagnum spp., graminoids, ericoids and other mosses) and key blanket bog features (pools and bare peat) in two near-naturel blanket bogs in Scotland. We examine empirical relationships between CO 2 fluxes and edaphic and atmospheric conditions to derive modelled annual CO 2 balance for the dominant PFTs during the 2018 drought and the following growing season of 2019. We report on variation observed between microtopes (margins and centre) and topographical setting (upland and lowland). We hypothesised that CO 2 fluxes would vary amongst the dominant blanket bog PFTs (Sphagnum spp., ericoids and graminoids) because of their functional traits; notably that they would vary in their productivity during drought conditions, and potentially during recovery period in the post-drought year. Differences in fluxes between the topographic settings are expected to relate to water availability variations during the sampling period.

| Site description
The study was conducted at two blanket bogs in northern Scotland, part of the Flow Country peatlands Lindsay et al., 1988) (Hancock et al., 2018) and with peat depths ranging from 0.5 to 4.3 m (Avercamp et al., 2021). Ericoid shrub species include Calluna vulgaris L. and Erica tetralix L., with a graminoid cover dominated by Eriophorum angustifolium L., Trichophorum germanicum L. and Carex panicea L. Common mosses include Sphagnum capillifolium, Sphagnum cuspidatum, Sphagnum fuscum and R. lanuginosum. E. angustifolium colonised peat hags and dry pools in wetter areas.
The lower altitude site Munsary (58 23 0 49.0 00 N, 3 20 0 26.5 00 W, approximately 100 m.a.s.l.) is located 30 km to the east of Knockfin and comprises low-lying blanket bog with peat depths between 2 and 5 m (Marshall et al., 2022). The mire expanse contains a wide range of micro-topographic elements such as hummock and hollows, with an extensive pool system in the centre (Smart, 1982). Sphagnum spp.
include Sphagnum papillosum, Sphagnum medium, Sphagnum austinii and S. capillifolium, while other bryophytes include Polytrichum commune and R. lanuginosum. For vascular species, C. vulgaris, E. angustifolium and Eriophorum vaginatum are abundant, together with T. germanicum and Narthecium ossifragum. Agricultural drains, historically cut in one area close to the southwestern margin of the site, were blocked with peat dams and plastic piling in the early 2000s.
At each site, six replicated sampling plots were set up, with no more than one sampling plot within each 100 Â 100 m subsite to capture variation between different margin and centre microtopes ( Figure 1) as identified by Marshall et al. (2022). Margin plots are associated with generally stiffer peat and are situated near the boundaries of the blanket bog expanse on the site (i.e., streams adjacent to Munsary and erosion gullies at Knockfin Heights). Centre microtope sites are found in the central, wetter and deeper peat areas, often with pool systems. At each plot, after an initial high-level (PFT and dominant cover type) vegetation survey in November 2017, collars were installed around three key peatland PFTs (Sphagnum spp., graminoids and ericoids) and where present, on lichen (Cladonia spp.), other mosses (Pleurozium. schreberi), bare peat (natural and anthropogenic originated erosion [Hancock et al., 2018]) and pools ( Figure 2) (Supporting information, Table S1). In order to allow for efficient CO 2 flux measurements and site-specific environmental condition monitoring, all collars were installed within an area of 5 Â 5 m that was representable for the margin or centre microtope at the plot.

| Flux chamber measurements
Between April-November 2018 and May-October 2019, CO 2 concentration-change was measured on 19 occasions using a mobile infrared gas analyser (IRGA; EGM-4, PP System) attached to a mobile cylindrical non-steady-state chamber (diameter of 19.2 cm, height 20 cm, volume 5.8 L) sealed to the ground-based collars. A domeshaped floating chamber with an air volume of 4 L was used to F I G U R E 1 Top: location of study sites in the upland (>300 m.a.s.l.) and low-lying blanket bogs in the Flow Country, northern Scotland. Map is derived from the Carbon and Peatland map (SNH, 2016), showing all features associated with blanket bog habitat and soil types. Bottom: distribution of sampling plots (100 Â 100 m, after Marshall et al., 2022) of different microtope positions across the upland site Knockfin Heights and low-lying site Munsary, overlain on a modified greyscale ©2022 Bing Maps satellite image for reference to highlight plot distribution measure pools. A fan was installed inside all the chambers to increase homogeneity of the sampled air during deployment.
To estimate ecosystem respiration (R eco ; total CO 2 respiration from plants, roots and soil microbes per unit ground area) and net ecosystem exchange (NEE) fluxes (the CO 2 respiration minus CO 2 fixation of the plant flux per unit ground area), the IRGA recorded CO 2 concentrations at an interval of 1.6 s  from the moment the chamber was put on the collar until the CO 2 concentrations stabilised, which took on average 2.65 min (between on average 2.88 min for Sphagnum spp. and 2.56 min for graminoids). For R eco measurement, the chamber was covered with a reflective shroud (i.e., dark chamber) that was removed for NEE measurement (clear chamber).
Between the two measurements, the chamber was removed from the collar and vented.

| Ancillary environmental measurements
To record environmental data between sampling, at each plot, data were collected at 30-min intervals for volumetric soil water content were recorded manually for all chamber measurements. All recorded soil parameters were measured at 5 cm depth. The EGM-4 also stores relative humidity (RH) and air temperature during deployment. Weather conditions and comments on on-site disturbances (e.g., deer trampling near/on collars) were also noted.

| Flux calculations
2.4.1 | In situ fluxes All CO 2 fluxes (μmol CO 2 m À2 s À1 ) were calculated by applying a linear regression or quadratic regression model in Python (version 3.6.3) using the statsmodels package (Seabold & Perktold, 2010). The convention of a negative rate indicating an uptake of CO 2 (i.e., sink) and positive rates for net emissions (source) has been adopted. The CO 2 concentration (ppm) data were first checked for outliers with an F I G U R E 2 Top: schematic example of collar set-up, indicating plant functional type (PFT) variability within a plot (note: actual PFT distribution and occurrence vary amongst plots). Bottom: margin plot (KH-A) at Knockfin Heights with water level logger (left), soil temperature and water content logger and five collars placed at dominant PFTs and key blanket bog features automated moving-window analysis to subsample the full dataset, retaining as many data points as possible, while omitting concentration outliers. Outliers include anomalously high, low or stable concentrations at the start or end of a measurement, indicating response-lag or equilibrium conditions, respectively (Pirk et al., 2016). Removing these points retained a sufficient period for actual flux calculations where the true rate is identified in the moving-window analysis (Supporting information, Figure S1).
When concentration flux is low, all data points are kept in the subsequent subset analysis, indicative of a net zero rate. A subset of data with the best linear model quality parameters was selected: that is, highest r 2 , lowest normalised root mean squared error (NRMSE) and highest Nash Sutcliffe efficiency coefficient (E). The subset was subsequently used for both linear and quadratic regression analyses, following the flux calculation procedures used by the IRGA (EGM-4, PP Systems). The rates were calculated using a discrete function (Parkinson, 1981;Widén & Lindroth, 2003): where the concentration linear flux rate F l is a function of the change in concentration C at time t i , the volume of the chamber V and the horizontal area of the chamber A. This function contains the linear equation used to determine flux rates, before adjusting for chamber dimensions and water vapour concentration changes: where change in concentration C is a constant linear function over time T.
Non-linear behaviours of both NEE and R eco fluxes have been recognised for chamber measurements (Kutzbach et al., 2007;Murphy et al., 2014;Pirk et al., 2016), resulting in the use of a quadratic function to estimate the CO 2 flux: where the change in concentration C over time T is constant at T = 0, resulting in a linear flux rate value that is not influenced by the derivative of the observer effect (2cT) (Wagner et al., 1997).
Model fit statistics were assessed to determine whether observed flux estimates fit a linear or quadratic function best and hence better represented ambient flux prior to measurement start. For the linear flux estimates, the model's coefficient of determination (r 2 ) generally provides a robust representation of the data when close to 1 (being a perfect linear regression line). Since higher-order polynomial regressions (i.e., the quadratic function) are not properly described by r 2values, they are compared using r 2 adj -values using the statsmodels Ordinary Least Squares function. Both linear and quadratic rates are calculated, and model parameters are stored for the following model-selection steps. Subsequently, rates are corrected for an increase of water vapour concentrations as dilution by water molecules can significantly reduce CO 2 flux estimates (Matsuura et al., 2011).
Confident linear rates are assumed for all NEE and R eco fluxes based on a linear model with r 2 adj -values >0.85 (Huttunen et al., 2002;Silva et al., 2015) (Supporting information, Figure S1). Additionally, linear rates with r 2 adj -values >0.75 are accepted (Strack et al., 2016) unless a concentration change of more than ± 2 ppm of the initial measurement level was observed. The quadratic model rate was used when quadratic r 2 adj -value exceeded 0.75, and the r 2 adj -value of the linear rate is lower than 0.75, together with a r 2 -value ≠ 0 of the moving-window analysis (indicating an absolute rate, greater than zero: jbj ) 0). R 2 adj -values of 0 for linear rates accompanied with a ±2 ppm range (indicative of low fluxes) are also accepted.
In total, 3% of the NEE and R eco fluxes were omitted from the dataset, and all other fluxes are best described by a linear rate. Depending on the final flux type, the value of b from either the linear [2] or quadratic function [3] was subsequently used to define NEE and R eco fluxes at time of the measurement. As the NEE chamber deployment immediately followed the dark R eco flux measurement, gross primary production (GPP; rate of CO 2 uptake through photosynthesis) was calculated as the difference between the two measurements.

| Potential maximum photosynthesis (P max )
Since all GPP fluxes (derived from NEE-R eco ) were measured under varying light conditions, P max was calculated for the dominant PFTs across both growing seasons using the Michaelis-Menten function (Laine et al., 2016): with P max as the potential maximum rate of photosynthesis (μmol CO 2 m À2 s À1 ), the half saturation constant k and the PAR measurement associated with the GPP flux during year i, on collar h at site g.
A non-linear least-squares (nls) model was used to estimate P max and k coefficients across the data, grouped by collar and year (function nlsList package nlme [Pinheiro et al., 2020] in R version 3.6 [R Core Team, 2014]). Subsequently, negative k coefficient estimates are omitted, along with the paired P max values; negative light levels at which half of P max is reached are not used in further analyses. Only two extreme P max estimations (À81 and À119 μmol CO 2 m À2 s À1 ) were left out due to their very high k estimates ()30,000 μmol m À2 s À1 ). A total of 34 P max and k coefficients of the collars are then used to describe functioning of the different PFTs between years, sites and microtopography.

| Estimation of growing season fluxes
The number of flux estimates allows for modelling of hourly GPP and R eco rates across the sampling periods in 2018 and 2019. Soil temperature, water table and PAR can be used to predict diurnal variability of fluxes  and offer an effective way to model annual NEE based on flux chamber measurements (Huth et al., 2017;Swenson et al., 2019). To be able to compare the carbon sequestering functioning to similar blanket bogs and peatlands, carbon fluxes are converted to g C-CO 2 m À2 h À1 . The focus was on Sphagnum spp., graminoids and ericoids, which had the highest number of measurements available for this modelling. Furthermore, to assess both temporal and spatial PFT flux behaviour, field measurements are divided based on site and microtope. To allow for appropriate siteand microtope-specific GPP and R eco estimations, each subset was accompanied by their respective hourly average soil temperature ( C), water level (cm below surface) and average PAR (μmol m À2 s À1 ).
where PAR is the photosynthetic active radiation (μmol m À2 s À1 ), T soil is the soil temperature ( C), W level is the water level (cm), which are all measured empirically. In the empirical model, Julian day of the year (J) was used to account for any seasonal variability of plant phenology, such as changes in green leaf area (Wilson et al., 2007) and model parameters a, b, c, d and e are specific to each PFT-site-microtope combinations (3 PFTs Â 2 sites Â 2 microtopes = 12 individual models) (Supporting information, Table S4).
A similar non-linear regression empirical approach was used to fit measured R eco fluxes: where model parameters a, b and c are unique to a PFT-sitemicrotope combination and the T soil and W level are the same as Equation (5) (Supporting information, Table S5).
Using both functions, subsets of the environmental data are used to predict hourly GPP and R eco for the months covering the sampling period (April 1st-November 31st) for both years to capture the majority of the growing season. As with the in situ flux processing, all positive GPP (3.9%) and negative R eco (1.2%) estimates were omitted from the dataset. Complete diurnal measurements were not obtained, and therefore, only daytime fluxes are modelled; the daytime (PAR > 25 μmol m À2 s À1 ) flux estimates, accounting for 53% of the modelled GPP and R eco , are used for growing season NEE estimates (NEE = GPP + R eco ). In order to obtain a confidence range of predicted fluxes, the residual standard error (RSE) of each model (predicted values are converted back to μmol CO 2 m À2 s À1 to allow for direct comparison with the in situ measurements) was used to calculate lower and upper estimates for the fluxes and subsequently the cumulative GPP and R eco , as well as NEE.
2.6 | Statistical analysis 2.6.1 | Comparison of in situ NEE, R eco , GPP and P max fluxes All statistical analyses were conducted using R (version 3.6) and the most recent versions of the packages at time of writing. Linear mixedeffects models (LMEMs) (function lmer from package lme4 [Bates et al., 2015]) with plot ID as a random effect were used to test whether NEE, R eco and GPP are different between PFTs and ground cover types, both within and between both years and between sites (upland vs. lowlying). For each LMEM, an ANOVA was applied to check the contribution to the models' variances for each fixed effect (ANOVA type III, function anova from stats package [R Core Team, 2014]). PFTs and ground cover types significantly contributed to the variance (p < 0.001) of the three LMEMs (NEE, R eco and GPP). Even though the interaction between PFTs and years is not significantly contributing to the variance in the R eco LMEM, the interaction factor was kept in all subsequent post-hoc tests as it is a significant contributor in the NEE and GPP models (p < 0.05). A pairwise comparison test (function glht from multcomp package [Hothorn et al., 2008]) was applied to each LMEM using a post hoc (Tukey) analysis to investigate differences in marginal means for all 'PFTs and ground cover -Year' pairs and was accompanied by letter-based representations of significantly different marginal means (p < 0.05) (function cld from multcomp package).
Only P max estimates of Sphagnum spp., graminoids and ericoids were used for testing differences between years, sites and microtopography. Small sample sizes (n < 30) of P. schreberi and Cladonia uncialis/impexa-mix impeded reliable coefficient estimates and were excluded together with the other ground cover types (i.e., pools and inundated/dry bare peat). Means of yearly P max estimates per collar were compared between Sphagnum spp., graminoids and ericoids using an LMEM with collar ID as a random effect (function lmer package lme4). A complex model (with the fixed-effects PFT, year, sites and microtopes) was fitted and compared to models with univariate models and interactions (i.e., PFTs across years). For all model configurations-including and excluding some predictors-the Akaike information criterion (AIC) was calculated (Burnham & Anderson, 2004). The most complex model has the best fit   [Lüdecke et al., 2020]). Furthermore, p-values for fixed effects are used as a guide to model selection: Simple models (<3 predictors) often include a higher number of fixed effects that significantly contribute to model variance compared to more complex models.

| Generalised additive mixed-effect models
To overcome any potential nonlinearity of environmental variation throughout the year that would limit the usability of LMEMs, generalised additive mixed-effects models (GAMMs) were used to test interactive effects of changing edaphic and atmospheric conditions. In contrast to the LMEM approach, variables are not scaled to allow for direct interpretation of the model output and seasonality of environmental conditions (i.e., for temperature and water content) and applied smooth terms to the variables to run the GAMMs using the mgcv package (Wood, 2017)  was used in all models, as univariate models performed poorer when using smooth constructions in the GAMMs. In each model, site location, plot and year were used as random effects to allow for individual slope estimates across both sites; microtope was left out as both fixed and random effects, as it did not improve the model fit criteria significantly either way. All models were fitted using the gamm4V function from the mgcViz package (Fasiolo et al., 2020), as it can handle both nested and crossed random effects (1jsite/plot) + (1jyear) (Wood, 2019) and allows for greater graphical freedom when plotting the output. Visual interpretations of the model output were produced with ggplot2 (Wickham et al., 2019), and covariate interactions were presented as fitted effects of NEE, R eco and GPP fluxes (Supporting information, Table S6).

| Comparing NEE
The growing season NEE estimates (g C-CO 2 m À2 y À1 ) were contrasted with results from previous work by comparing sites, across a range of peatland types, summarised by Swenson et al. (2019). Both flux chamber and eddy covariance (EC) measurements were included to provide a comprehensive overview of CO 2 balance estimates across boreal and temperate peatlands. These peatland types include (near-)natural blanket bogs; other bogs (such as oligotropic and ombrotropic peatland sites that also have not been impacted by any disturbances); a range of poor fens; sites with bare peat (resulting from peat harvesting or erosion); and peatlands which have undergone some kind of restoration management (e.g., drain-blocking or rewetting). As both Knockfin Heights and Munsary are identified as 'near-natural' blanket bogs, comparisons with afforested sites (including restoring sites after tree-felling), rich-fens and tropical peatlands have been excluded. Mean annual water table was used to compare between sites and years as this aspect reflects annual hydrological conditions on a site and microtope scale and can be linked to annual NEE estimates (Gažovič et al., 2013;Helfter et al., 2015;Kritzler et al., 2016).  Figure S2). The precipitation deficit was not limited to the drought period as lesser precipitation (38%-79%, 35%-83% and 29%-76%) was recorded across the spring, summer and autumn months, respectively. Mean summer (June-August) air temperatures recorded for Northern Scotland in 2018 (13.1 C) and 2019 (12.8 C) were comparable, although mean daily maximum temperatures recorded over these months were almost 1 C higher in 2018 (17.0 vs. 16.3 C) (Supporting information, at both sites (January-October) was lower in 2018 than in 2019 (7.9 vs. 8.9 C for Munsary; 6.6 vs. 6.9 C for Knockfin Heights).
The precipitation deficit observed in 2018/2019 impacted soil water content and water level distributions across both years. At Munsary, the decline in water content was greater in margin sites compared to centre plots during the drought period (Figure 3e). At Knockfin Heights, higher water content was recorded in the margin sites compared to centre plots after the drought period. Munsary F I G U R E 3 Environmental data overview for Munsary and Knockfin Heights: daily average (±SE) of soil and water temperature (a,b); monthly cumulative photosynthetic active radiation (PAR) (c,d), daily average (±SE) soil water content for centre and margin plots (e,f), monthly cumulative precipitation (g,h) and daily average water level (i,j). Drought period (24 May Table S2). At both sites, water levels showed less variability during the drought period than during the spring, autumn and winter months (Figure 3i,j).

| In situ flux dynamics
A total of 1114 unique CO 2 flux measurements obtained from the blanket bog study sites and GPP flux rates were derived from the paired NEE and R eco flux estimates across the three dominant blanket bog PFTs ( Figure 5) and other mosses, lichen, pools and bare peat (Supporting information, Figure S3).
There was no significant effect of study site and microtope on the variance across the fluxes for all PFTs and ground cover types (F F I G U R E 6 Inter-annual mean (±SE) for all CO 2 fluxes from both sites, as: net ecosystem exchange (NEE); ecosystem respiration (R eco ); and gross primary production (GPP) for both years across the plant functional types (PFTs) and ground cover types. Amongst flux types, shared letters are used for groups not significantly different (p > 0.05), number represents sample size.
F I G U R E 7 Mean estimated potential maximum photosynthesis (P max ) ± SE for the three dominant plant functional types (PFTs) (Sphagnum spp., graminoids and ericoids) for both years. Shared letters indicate no significant differences (p > 0.05) between and across PFTs and years (Tukey post-hoc test). and combined effect (F(2, 9.015) = 9.548, p = 0.006)).

| Growing season GPP and R eco
The predicted GPP and R eco values showed strong linear correlations with the measured CO 2 fluxes (Supporting information, Figures (Table 1).

The increased range in NEE estimates on Knockfin Heights for
graminoids is the result of higher predicted maximum net uptake of CO 2 in 2019 compared to the year before on the margin and centre plots (Figure 9). Similarly, at Munsary centre plots, predicted annual net CO 2 flux for graminoids was lower in 2019 but the margin estimates show an opposite trend with an increase in net CO 2 emissions in 2019. This decrease in CO 2 sequestration is also observed in ericoid NEE estimates, but only on Knockfin Heights. However, both margin and centre plots are still likely to be net C-CO 2 sinks.
Switches between overall source and sink behaviour are observed on Knockfin Heights, where Sphagnum spp. on the margin plots moved towards a stronger C-CO 2 sink, but at the centre plots, the opposite occurred, with higher emission in 2019. Consistent CO 2 sources for both years are found in the Sphagnum spp. NEE estimates for Munsary, T A B L E 1 Ranges in cumulative predicted net ecosystem exchange (NEE) across each site for both the drought (2018) and post-drought (2019) years using the maximum and minimum NEE (Figure 9). Change in ranges quantified as percentage increase or decrease in 2019 compared to the NEE ranges in 2018 with almost no changes between both years across margin and centre plots. Conversely, in the ericoids at Munsary, centre plots acted as sources and margin plots net sinks for both years (Figure 9).

| Environmental controls on in situ fluxes
GPP, R eco and NEE fluxes in this study were driven by environmental conditions, and this relationship differed between the PFTs. The best performing models with Year included as a fixed effect had a soil temperature component to explain CO 2 exchange variability (Table 2).
Fluxes in the models are generally related to increases in soil temperature, or daily soil temperature range in the case for ericoid NEE (pseudo r 2 = 0.40) and GPP (pseudo r 2 = 0.52). Only the graminoid NEE model did not include a temperature component, but rather a single PAR covariate to account for variability in net CO 2 exchange (pseudo r 2 = 0.32). The daily mean and range in peat water content also accounted for variance amongst flux estimates for Sphagnum spp.
(R eco pseudo r 2 = 0.59; GPP pseudo r 2 = 0.65), with similar, but stronger effects on R eco fluxes for graminoids (pseudo r 2 = 0.65) and ericoids (pseudo r 2 = 0.59). Changes in daily mean water level accounted for graminoid GPP variance, with increasing uptake of CO 2 associated with higher water levels. However, variation was still predominantly controlled by soil temperature (pseudo r 2 = 0.96). As a F I G U R E 1 0 Interactive fitted effects of environmental covariates on net ecosystem exchange (NEE), ecosystem respiration (R eco ) and gross primary production (GPP) for Sphagnum spp. (a, b, c), graminoids (d, e, f) and ericoids (g, h, i). Only fitted flux values within the range of in situ environmental measurements are shown. Covariates include photosynthetic active radiation (PAR) (in situ or loggers), soil temperature (in situ or daily range from loggers), water content (daily mean or range from loggers) and water level (daily average).
fixed effect, Year was only a significant predictor of graminoid NEE and ericoid R eco and GPP fluxes but was kept in all other models to account for variability across both sampling campaigns. Interactions between fixed effects did not improve model predictions with microtope type only applied to the Sphagnum spp. NEE model, where centre plots appear to be a more important predictor than Year for the period of sampling (pseudo r 2 = 0.31).
An optimal moisture regime is not clear in ericoids R eco but trends towards higher fluxes with drier soil conditions (Figure 10h). These relationships between temperature and water availability are replicated by GPP; however, the effect of drier conditions is more pronounced. GPP rates for Sphagnum spp. generally increase with soil temperature (above 12 C) when daily fluctuations in water content are above 0.04 m 3 /m 3 d À1 (limited to days of high soil temperature and low daily precipitation) (Figure 10c). Graminoid GPP was highest during days with shallow water levels (<À0.15 cm) and high soil temperatures (>12 C) (Figure 10f). However, these rates are restricted to the warmer seasons, during which almost no inundation occurred (daily soil temperature >11 C). Ericoid GPP was less well defined by the best model's covariates, soil temperature and average water con-  (Saunders et al., 2021) and Sweden on an ombrotrophic blanket bog (Keane et al., 2020), although emissions recovered in 2019 to pre-drought rates. Though not measured in this study, higher temperatures likely also contributed to increased microbial respiration (Järveoja et al., 2020;Keiser et al., 2019).
For non-vegetated sites (bare peat collars), the steady overall net CO 2 emissions measured were in line with other sites in the UK (Clay et al., 2012;Dixon et al., 2014;Gatis et al., 2019). The slight increase in CO 2 emissions observed in inundated bare peat collars and floating chamber measurements in 2019 could reflect increased mineralisation associated with biogeochemical and microbial cascades associated with the 'enzymic latch' observed in peatlands following droughtrewetting cycles (Fenner & Freeman, 2011). In pools at Munsary, increased mineralisation rates have been associated with high DOC concentrations and correlated with higher dissolved CO 2 concentrations (Turner et al., 2016) and net CO 2 emissions. For vegetated sites, earlier onset of growing season (associated with higher spring temperatures) is not expected to influence cumulative GPP for northern peatlands (Kross et al., 2014). However, it can negatively affect annual NEE, when short, severe drought periods in peak growing season increase R eco , but not GPP (Lund et al., 2012). This is observed for  Figure S2), Knockfin Heights experienced a smaller anomaly towards the end of the sampling period, mitigating some of the early drought impact from the dry spring months. This difference is in line with the east-west gradient of increasing rainfall across the Flow Country region (Lindsay et al., 1988). Importantly, blanket bogs rely on occult precipitation, which can contribute up to 20% of their water input (Lapen et al., 2000). It is therefore possible that at 350 m above sea level, Knockfin Heights experienced more frequent fog and increased dewfall, effectively reducing the site water deficit and avoiding Sphagnum spp. 'shut down'.
The drought of 2018 led to surface subsidence at both Munsary and Knockfin Heights of between 5-10 and 1-2 cm, respectively Marshall et al., 2022). Lowering of the surface following water level drops in the peat is predicted to mitigate against drought conditions by allowing higher water content to be sustained at the peat surface (Price, 2003) and increasing water availability for the blanket bog vegetation (Lapen et al., 2000). Peat surface subsidence and its potential to buffer against water level drawdown could explain why some of the PFTs were still overall CO 2 sinks during the 2018 period, as water deficits may not have reached critical thresholds, (e.g., graminoids at margin plots on Knockfin Heights [ Figure 9]).
Since water level is used for the modelled cumulative GPP and R eco estimates, high subsidence rates at Munsary resulting in the peatland surface effectively tracking the water table could explain similar annual estimates for both drought and post-drought years. For Munsary, this explains NEE estimates for all but the graminoids in the margin plots, where annual estimates switch to net CO 2 sources.
Although Knockfin Heights experienced more variable water levels and marginally drier conditions, Munsary margin plots experienced more frequent and prolonged lower soil water contents during 2018 (Figures 3 and 4). This supports the observation that blanket bog margins show higher CO 2 flux due to higher environmental variability (i.e., water level and subsidence rates) in response to drought. It is possible that graminoids in margin areas at Munsary, identified as areas of high subsidence rates (Alshammari et al., 2018;Marshall et al., 2019;Marshall et al., 2022), could be responding (lower GPP, higher R eco ) to structural alterations of the peat surface (i.e., cracks and bare peat exposure) in the post-drought year.
4.2 | Relating drought effects on CO 2 fluxes to environmental conditions Drought effects on fluxes across the PFTs can be summarised by the effects of increased soil temperatures and changes in water availability.
Ericoids showed a stronger increase in both CO 2 sequestration and respiration rates (GPP and R eco ) compared to Sphagnum spp. and graminoids in the post-drought year. Acknowledged is that the increased R eco fluxes could also be explained by increased methane oxidation during drier upper peat condition (Freeman et al., 2002), though this was outside the scope of this study. The resilience pattern of ericoids > graminoids > Sphagnum spp. compares well with results from other ombrotrophic bog studies (Bubier et al., 2003;Goud et al., 2017), due to the physiological advantage of vascular plants to better stabilise their CO 2 fluxes during drier conditions by actively taking up water and limiting transpiration. Conversely, increased water availability promoted Sphagnum spp. GPP and graminoid GPP (dominated by Eriophorum spp.) in line with their preference for overall wet conditions in blanket bogs  and peatlands in general (Cooper et al., 2014;Van Breemen, 1995). Increased water availability and higher water levels also increase R eco , something also observed in shallow peat cores from a site 10 km west of Knockfin Heights (Hermans et al., 2019) and a rewetting experiment on similar PFTs (Kuiper et al., 2014).
Accounting for different light intensities, drought conditions did not have a significant effect on the productivity estimates (P max ) for Sphagnum spp. and ericoids. Uptake estimates in 2018 were lower for graminoids ( Figure 7) with modelled maximum GAMM GPP compared to P max estimates (À6.7 against À12.8 μmol CO 2 m À2 s À1 ), indicating that the graminoids had not returned to optimal condition in 2019 and remained limited in their productivity (Strachan et al., 2016). Similarly, both measured and modelled ranges for Sphagnum spp. and ericoid GPP in 2019 remained below P max estimates. The general lower light saturation point of bryophytes (Marschall & Proctor, 2004) and comparable PAR regimes across 2018 and 2019 (Figure 3c,d) could explain the consistency of maximum productivity estimates. The uniformity of in situ fluxes and P max estimates suggests that conditions were not (yet) limiting Sphagnum spp. productivity, evidenced by the limited Sphagnum bleaching observed in the field. Overall, the response of Sphagnum spp. to reduce water levels was likely less severe than hypothesised, for example, a substantial reduction in productivity (Jassey & Signarbieux, 2019;Strack et al., 2009).
PAR was only a dominant driver for ericoid GPP, suggesting that ericoid photosynthetic activity was not limited by water availability.
NEE estimates for a range of C. vulgaris age/height ranges show that annual sinks are not common   ± 48.0 g C-CO 2 m À2 yr À1 ), other bogs (À21.0 ± 91.9 C-CO 2 m À2 yr À1 ) and fens (À70.3 ± 39.0 g C-CO 2 m À2 yr À1 ) ( Figure 11). Most (blanket) bog NEE estimates are from EC measurements, with large footprints, whereas CO 2 estimates from chamber-based studies include microscale targets. Nonetheless, the multi-annual estimates from the studies involved allow for comparisons with MAWL variability in northern peatlands (Swenson et al., 2019).
MAWL for the upland and low-lying site (Supporting information, Table S2) is also comparable to that of intact/near-natural sites (À9.6 ± 8.3 cm), mostly staying above the MAWL of bare peat sites (À24.5 ± 17.5 cm) and drained/restored sites (À26.0 ± 20.7 cm). Net C-CO 2 sink and source behaviour related to MAWL is visible in the other peatland types where annual and seasonal NEE (g C-CO 2 m À2 yr À1 ) estimates show a general increase in net emission with lower MAWL ( Figure 9). Generally, lower water levels-in these cases linked to peatland type and restoration or management history-correspond with higher overall C-CO 2 emissions. This relationship was estimated to be strong for all vegetated sites and less so for the bare peat sites, being overall CO 2 sources (Swenson et al., 2019). However, MAWL between the drought and post-drought year only explains a minor part of the water deficit differences across the plots in this study, and water-stress estimates for PFTs are better reflected in observed soil water content relationships (GAMM results).

| Future blanket bog drought resilience
The CO 2 balance of the PFTs and key blanket bog features in this study shows that the Flow Country peatlands are sensitive to drought conditions (Fenner & Freeman, 2011;Lund et al., 2012). Changes in the long-term functional state (net sink vs. net source) of blanket bogs are expected to depend on the frequency and duration of predicted droughts in Scotland (Grillakis, 2019). Repeated droughts are likely to trigger a shift to more drought resilient vascular plant species and Sphagnum spp. replacing other species occupying similar ecological F I G U R E 1 1 Left: comparison of annual/seasonal net ecosystem exchange (NEE) estimates plotted against mean annual water level for global studies on boreal and temperature peatlands, including bare peat sites, (near-)natural blanket bogs, fens and other bogs, drained and restored sites (adapted from Swenson et al., 2019) and the modelled estimates from this study. Number in brackets in legend signifies data point count for each category in the plot. Right: seasonal NEE estimates for Sphagnum spp., graminoids and ericoids for both sites, across margin and centre plots for 2018 and 2019: C-CO 2 ranges are left out for simplicity, please refer to Figure 9.
niches, which are unable to adapt (Robroek et al., 2017). For example, shifts from Sphagnum spp. to more resilient bryophyte species during prolonged periods of droughts such as Polytrichum (Potvin et al., 2014) or R. lanuginosum (Lindsay, 2010) could also be expected.
Changing blanket bog vegetation composition during and after drought conditions (i.e., upland graminoids and ericoids) also results in higher biomass and litter load (Grau-Andrés et al., 2018), which alongside changing climate will make peatland more vulnerable to wildfires  and fire-related emissions. This coupled with drought-induced structural peat surface alterations-that is, cracks and drier conditions (Li et al., 2018)-also helps facilitate prolonged smouldering, increasing peat fire hazard (Davies et al., 2013). This risks creating positive feedback loops where drought and fire frequency increases at the cost of peatland resilience (Field et al., 2007), offsetting any gains in productivity.
Together with novel remote sensing techniques (Alshammari et al., 2018(Alshammari et al., , 2020Fiaschi et al., 2019;Marshall et al., 2022;Tampuu et al., 2020) and accurate land-cover estimates, the flux estimates from the margins and centre plots at the upland and low-lying site provide insight into the blanket bog response characteristics across a range of dominant vegetation types during and after an extreme drought. The ability to differentiate between PFT and key blanket bog features highlights the unique response each type is likely to have in response to a predicted shift towards drier and warmer summers and the impact on carbon fluxes on blanket bogs across the United Kingdom.

| CONCLUSION
Overall, our study shows that estimating CO 2 fluxes of PFTs and key blanket bog features is valuable to understand different responses during and after a drought event. With these growing season flux estimates, we were able to show a drought resilience pattern of ericoids > graminoids > Sphagnum spp. at both upland and low-lying sites. The relation to water availability and their topographic setting suggests that if drought frequency and severity increase, blanket bog vegetation is expected to shift towards assemblages that can cope with these conditions (i.e., ericoid shrubs). Since these hydrological conditions are expected to become more frequent, it is important to establish baseline studies to understand the long-term impact of these extreme events (e.g., drought and wildfires) on carbon fluxes. This will require upscaling using vegetation as a proxy for GHG emissions to assess landscape scale drought-responses, post-fire recovery studies or even national and global climate impact predictions and requires integration with remote sensing measurements, including the deployment of unoccupied aerial vehicles (UAVs) for increased detail and more long-term monitoring on peatlands to assess the effectiveness of conservation management.

ACKNOWLEDGEMENTS
HPS was funded through a PhD project supported by the European Social Fund and Scottish Funding Council as part of Developing Scotland for supporting us with field preparations. We also thank Jasmijn Sybenga, Natalie Isaksson and Julia Avercamp for helping in the field.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are openly available in NERC Environmental Information Data Centre at https://doi.org/10. 5285/1be4eef2-0591-4073-bae2-e00b6ff4462f.