Examining anomalous summer carbonyl sulfide emissions in a boreal forest after thinning

Carbonyl sulfide (COS) is an atmospheric tracer for ecosystem carbon and water fluxes because the diffusive pathway for COS is similar to that for carbon dioxide (CO₂)^1,2,3,4,5. COS, along with CO₂, is taken up through the stomata to the chloroplast surface, where COS is irreversibly hydrolysed by carbonic anhydrase (CA)^6,7,8. While CO₂ is emitted back to the atmosphere during plant respiration, COS emission by leaves is not observed under ambient conditions^7,9,10. The key assumption behind using COS as a proxy for ecosystem-level gross primary production (GPP) is that photosynthesising plants are a major terrestrial sink of COS.

In the boreal forest ecosystem, canopy uptake is typically the largest COS flux, even at night, followed by soil uptake¹¹. The incomplete closure of stomata and the light-independent nature of CA lead to the nighttime uptake of COS^12,13,14, which correlates with stomatal conductance^11,15. During the day, stomatal conductance linearly increases with light intensity in a limited range until it plateaus, and COS uptake could be used to derive canopy and stomatal conductance, given that soil fluxes could be distinguished from separate measurements^13,16,17.

In order to accurately use COS as a GPP proxy, all the non-stomatal sinks and sources have to be quantified. The low seasonal variation of COS soil fluxes in the boreal forest makes it easier to characterise the canopy uptake and account for ecosystem budgets and GPP estimates¹⁸. Soil COS chamber measurements at Hyytiälä forest station have shown an average flux of ~ −3 pmol m⁻² s⁻¹ during the growing season, with weak seasonal changes¹⁸. Soil uptake is linked to the presence of CA in microbial communities^19,20, necessitating optimal soil temperature and moisture conditions for the maximal absorption of COS by soil^{18,21,22,23,24}. The contribution of soil to the global COS sink is estimated to be 26–33%^1,25 but is lower in boreal forests (10–20%)^5,18.

Wetlands and anoxic soils are considered terrestrial sources of COS^20,26,27,28, while oxic soils are generally inferred to be COS sinks from measurements^18,19,29 and modelling studies^24,30,31. However, exposure to high temperatures (above 15 °C) and high insolation in dry conditions can result in COS emissions from oxic soils^{14,21,32,33,34}. Soil COS emissions are primarily linked to the oxic and abiotic degradation of organic matter within the soil³⁴, and in the case of anoxic soils, redox reactions are important^21,26. Furthermore, emission of COS due to exposure to UV radiation was shown to vary with the amount of soil organic matter^34,35. Agriculturally managed temperate mountain grasslands turned into a net source of COS shortly after grass cutting, and the soil contribution could not explain the net ecosystem level emissions³⁶. Commane et al.³⁷ observed unexpected net COS emission over a mid-latitude forest site during the warmest summer weeks.

To better understand and quantify the source and sink processes in various ecosystems, and thereby improve the accuracy of GPP estimates, COS measurements in a variety of ecosystems at longer temporal scales are needed^1,38,39,40. The bottom-up global estimate of COS is not closed, suggesting an unaccounted source of COS in the tropics⁴¹. Inverse modelling studies have also shown a missing source in the tropics and a missing sink in the high latitudes^42,43,44. Parameterising seasonal and interannual variation of the boreal forest COS uptake in the Simple Biosphere Model Version 4 (SiB4) was able to match the missing COS sink to a certain limit in the high latitude, but the regional variation remains unknown⁵. A recent measurement has found boreal wetlands to be a sink for COS⁴⁵, which was previously reported as a source^4,26,27, underscoring the complexity of ecosystem processes in the global COS budget.

COS fluxes at the ecosystem level can be studied using the eddy covariance (EC) method at high temporal resolution. Several EC flux measurements of COS and attempts to estimate GPP from them have been reported in the literature^{5,13,37,46,47}. The most extensive EC measurements of COS flux from a boreal forest were conducted at the Hyytiälä forest station, spanning over 32 months across five years from 2013 to 2017⁵. This multi-year dataset showed that monthly median COS uptake varied by up to 26%. The main controlling bioclimatic factors of COS fluxes at the ecosystem level were photosynthetically active radiation (PAR), vapour pressure deficit (VPD), air temperature (T_a), and leaf area index (LAI).

In this study, we measured and analysed three more years (2020–2022) of COS flux measurements at the Hyytiälä forest station following a thinning conducted in 2019–2020 that removed ~ 40% basal area. Building on the well-studied COS dynamics of Hyytiälä forest station, we leveraged this unique ecosystem manipulation in the boreal forest to better understand the role of different ecosystem components in COS fluxes. Thinning altered the microclimate of the forest floor and created a large amount of cutting residue that was left on the site⁴⁸. The post-thinning COS fluxes exhibited anomalous behaviour, and in the summer of 2021, we observed continuous and long-lasting ecosystem-level COS emissions from the boreal forest. Here, we quantify the emission of COS, identify the various potential sources of COS and understand the drivers that controlled the anomalous COS emission in the boreal forest ecosystem.

Anomalous COS emissions and meteorological conditions

Thinning reduced forest canopy COS uptake, determined from the net ecosystem COS fluxes (F_COS) measured by EC. F_COS was ~19 % lower in 2020 (−19.6 ± 13.4 pmol m⁻² s⁻¹; mean ± standard deviation) and nearly 24 % lower in 2022 (−18.4 ± 13.1 pmol m⁻² s⁻¹) compared to the pre-thinning five-year (2013–17) average of −24.2 ± 13.3 pmol m⁻² s⁻¹. However, F_COS was about 75 % lower in 2021 (−6.0 ± 14.5 pmol m⁻² s⁻¹) due to anomalous COS emissions that were observed from June to mid-August and peaked in June and July (Fig. 1). Consequently, this study focuses on these months to understand and isolate the contributing conditions. The seasonal variation (F_COS) remained consistent in the post-thinning years, with maximum uptake in summer months and minimal uptake in winter months (Supplementary Fig. 1).

Examining the environmental variables, June–July 2021 stands out as warmer and drier compared to other years (Fig. 1 and Supplementary Fig. 2). The monthly diel variation of T_a, VPD and soil temperature (T_s) shows that they stayed well above their interquartile range (IQR) from previous years. June 2020 was also warmer and drier than usual, but July and August returned to previous IQR values (Supplementary Fig. 2). PAR remained within average limits during the study period (Fig. 1 and Supplementary Fig. 2) and it peaked during June and July with values ~ 1700 μmol m⁻²s⁻¹. T_a ranged from −24.7 °C to 30.9 °C during the 3 years of study, and the highest temperatures were measured during June 2021. Similarly, the highest VPD (3.18 kPa) and T_s (16.2 °C) were also measured in the same month. The forest thinning, which happened in January–March of 2020, opened the forest canopy crown, drastically increasing the light availability at the forest floor⁴⁸. The surface layer soil temperature T_s,humus increased on average by 2 °C in the post-thinning summer months. The precipitation in June and July 2021 was lower than in previous years (Supplementary Fig. 2f), which, in combination with the higher T_s, led to a reduced soil water content (SWC) at the site (Supplementary Fig. 2e). In previous years, SWC peaked in April-May and the minimum occurred in August-September. In 2021, SWC was generally low between June and August, with the minimum in July (Supplementary Fig. 2e). The mean summer SWC was 0.20 ± 0.03 m³m⁻³ (mean ± standard deviation) during the study period, but the values dropped below 0.14 m³m⁻³ for a short period in the last 2 weeks of July 2021. However, the water potential at the mineral soil surface dropped to −0.5 MPa on 28 July 2021, which could be considered a mild drought⁴⁹. In Hyytiälä, the forest generally acts as a sink for COS, with interannual variability dependent on environmental conditions. A reduction in uptake was reported in July 2014, when T_a and VPD were higher⁵. Similar COS emissions were observed in summer over Harvard Forest under high VPD and T_a conditions³⁷.

To understand the ecosystem-level responses and isolate the meteorological influence on F_COS during June and July, F_COS is presented as a function of PAR, T_a and VPD (Fig. 2). The relationships between F_COS and environmental variables during 2020 and 2022 were similar to the pre-thinning years⁵, showing a decrease of COS uptake with an increase in T_a and VPD, and remained relatively stable under high PAR (>1200 μmol m⁻²s⁻¹). In 2021, the responses were different, the uptake was lower and net emissions were observed at PAR > 1300 μmol m⁻²s⁻¹, T_a > 24 °C and VPD > 1.8 kPa. A quadratic function was fitted to the F_COS responses (Supplementary Fig. 3), to determine if the meteorological responses differed significantly between years (p < 0.001). The functions were fitted separately for each year and then combined for all post-thinning years. The fitting coefficients for the quadratic functions are different for 2021 when compared to other years (Supplementary Table 1). The regression metrics (Supplementary Table 1) show that the year-specific model fits were able to explain the variance very well (R²: PAR = 0.94, T_a = 0.93 and VPD = 0.91) compared to one fit for all years (R²: PAR = 0.10, T_a = 0.44 and VPD = 0.39). The improved fit highlights that the interannual variation of the flux responses was significant (p < 0.001).

Reduction in canopy uptake during extreme conditions

The changes in the meteorological conditions could alter the canopy conductance (G_C) and therefore the F_COS. The G_C plotted against each season’s low, medium, and high percentiles of PAR from 2020 to 2022 (Fig. 3a) shows that the G_C was slightly lower (~12 − 17%) in the 2021 summer. The COS uptake is proportional to the stomatal conductance¹¹, which was observed for all the years, including 2021. The slope of the line fitted against F_COS and G_C is similar in all the years (Fig. 3b), with only the intercept being positive in 2021. The intercept shows the F_COS when the canopy conductance is zero, i.e., the contribution from forest floor and soil. The positive intercept points towards an emission source, which offset the non-canopy sinks in 2021.

The influence of meteorological conditions on the reduced COS uptake can also be explored through a meteorologically-driven ecosystem COS flux model. The parameterisation for the daily scale net COS forest sink presented by Vesala et al.⁵ could be used to quantify the changes in the canopy uptake. The parameterisation considers PAR, T_a, VPD and LAI to estimate COS uptake, reflecting its sensitivity to stomatal conductance (R² = 0.73, p < 0.001). The parameterisation scheme could closely follow the observed fluxes on a daily scale and effectively capture seasonal variations. In this study, the parameterisation was applied to the Hyytiälä forest stand, accounting for the thinning-related drop in LAI (Supplementary Fig. 4). The average parameterised canopy uptake for summer months was estimated to be −12.4 ± 4.1, −11.2 ± 2.8 and −13.7 ± 3.4 pmol m⁻² s⁻¹, for 2020, 2021 and 2022, respectively. The parameterisation generally underestimated the uptake compared to the EC measurements (2020: 34% and 2022: 20%), possibly due to the underrepresentation of forest floor vegetation in the LAI. However, in 2021, the parameterisation overestimated the COS uptake by a factor of two compared to EC measurements. The parameterisation fails to capture the seasonal pattern observed in 2021 (Supplementary Fig. 5). The changes in parameterised fluxes were observed each year, owing to changes in stomatal conductance and, thus, the daily canopy uptake. The average parameterised fluxes for the June-July months demonstrated a small reduction in canopy uptake in 2021 (~1 pmol m⁻² s⁻¹) compared to 2020 and 2022.

The parameterisation and the changes in the G_C both showed that the uptake of COS has undergone slight changes (~5% and ~ 10%, respectively) over the years, which did not explain the observed emissions. When the parameterised COS fluxes are compared with the measured EC daily mean for the 2021 summer (Supplementary Table 2), there was an unaccounted ~ 8 pmol m⁻² s⁻¹ COS flux in the budget of this ecosystem.

Changes in soil fluxes

In the absence of direct soil chamber measurements in the post-thinning years (2020–22), a depth-resolved diffusion-reaction model (COS Soil Model, COSSM) was used to simulate soil COS fluxes from meteorological and soil physical inputs³⁰. The thinning opened the canopy and warmed the forest floor compared to the pre-thinning years (pre-thinning T_s,humus: 12.3 ± 2.6 °C; post-thinning T_s,humus: 13.5 ± 2.2, 15.5 ± 2.7, 13.7 ± 2.6 °C during 2020, 2021, 2022, respectively) which can significantly influence the soil COS dynamics. The average modelled net soil COS fluxes for June-July were −2.3 ± 0.8, −1.7 ± 1.4 and −2.0 ± 1.0 pmol m⁻² s⁻¹, respectively, for 2020, 2021 and 2022. The higher temperature during the 2021 summer led to an increase in the soil COS production rates, thereby reducing the net soil fluxes (Supplementary Fig. 6). The modelled net soil fluxes were still low compared to total ecosystem-level fluxes, which were in the range of −3.7 ± 18.4 pmol m⁻² s⁻¹(F_COS, June-July: 2021). It could not account for the total COS emission observed in 2021.

The changes in the soil COS fluxes in response to the warmer forest floor alone cannot explain the observed emissions. The COSSM used the soil chamber dataset of Hyytiälä from the growing season of 2015¹⁸ for fitting the model. The soil chamber measurements during 2015 were on bare soil without a moss or litter layer, and the temperature throughout the measurement was lower than 18.0 °C. Aslan et al.⁴⁸ showed that in Hyytiälä, light penetration onto the forest floor was much lower in the pre-thinning years (~0.2) compared to post-thinning (~0.4). These factors could contribute to the uncertainty in model-predicted soil COS fluxes. Nevertheless, an increase in soil COS production is most likely in the post-thinning years. At high temperatures and radiation exposure, soils have been found to emit COS^14,32,33. Soil dryness could also enhance abiotic COS production^14,34. Wheat field soil shifted from consumption to production in the dark at high (>25 °C) temperatures, and sterilised soil shifts from uptake to emissions in the dark measured at a temperature of 19 °C³³. Sulfur-abundant agricultural soils emitted COS at ~ 20 °C in dark conditions³⁴. Bryophytes turned from a sink of COS in the dark to a source in light conditions⁵⁰. Even though other ecosystems have demonstrated COS soil production, it is unlikely that the soil alone could explain the ecosystem-level emissions at Hyytiälä.

Spielmann et al.³⁶ have shown that COS emissions from mountain grasslands after harvest hint towards possible COS sources from dead plant matter degradation. Meredith et al.³⁴ found that the microbial cycling of sulfur-containing amino acids, which are COS precursors, was the underlying cause of COS soil fluxes in soil samples from 58 sites. Their study further found that the increase in temperature could explain the apparent light effect, and maximum emissions were observed in boreal peatland soils. The exposure of litter to UV radiation has been found to increase the photodegradation by direct photolysis and also facilitate the microbial degradation by making the organic compounds bioavailable^51,52. Solar radiation can also cause thermal and photodegradation of COS precursors analogous to the processes observed in marine ecosystems^53,54. Differentiating the direct photolysis and microbial activity is difficult, but UV exposure has consistently enhanced COS emissions or decreased uptake by soil and litter samples, particularly those with high organic matter content⁵⁵. COS photoproduction from the soil and litter layer was omitted in the COSSM soil model, which could contribute to the observed anomalous emissions.

Nighttime fluxes

The light independence of plant CA enzyme and incomplete closure of leaf stomata lead to nighttime canopy uptake of COS^11,14. Soil chamber measurements from Hyytiälä have shown that nighttime soil uptake was −2.7 pmol m⁻² s⁻¹ during the measurement in the 2015 growing season¹⁸. The forest was a nighttime COS sink (−12.4 ± 10.6 pmol m⁻² s⁻¹) during the 2013–2017 measurement period. The mean nighttime uptake during 2021 (−3.2 ± 14.4 pmol m⁻² s⁻¹) was less than other years (2020: −10.7 ± 11.7 pmol m⁻² s⁻¹ and 2022: −9.0 ± 14.4 pmol m⁻² s⁻¹), demonstrating reduced ecosystem uptake. Kooijmans et al.¹¹ have shown that the nighttime uptake of COS is controlled by stomatal uptake and reduces with an increase in T_a and VPD. The relation of nighttime F_COS to T_a and VPD (Supplementary Fig. 7) shows this dependence across all years, i.e., reduction of COS uptake as T_a and VPD increase. Even though the median COS nighttime uptake in 2021 was smaller than in other years, uptake still happened at night. The higher average nighttime T_a reduced the stomatal conductance, which would partially explain the reduction in nighttime COS uptake. However, in similar T_a and VPD ranges, the COS uptake was lower in 2021. COS emissions from litter due to the reminiscent heat on the forest floor could not be ruled out, as small emissions were observed at higher temperature conditions at night. Therefore, the emission mechanism could have a thermal component along with a photolytic component that was more active in the daytime.

COS emission from a disturbed forest stand

The thinning in 2020 disturbed the ecosystem, reduced the GPP by 20% compared to the long-term mean and made it a net source of CO₂⁴⁸. The forest stand was able to recover to long-term mean values by 2022⁴⁸. The recovery of the forest was reflected in the response of GPP estimated by partitioning the NEE (GPP_NEE), as shown in Fig. 4a. Further, we compare the response of GPP from COS (GPP_COS) in the post-thinning years with the pre-thinning years (Fig. 4b). The responses were noisier than GPP_NEE, and the recovery of the forest from the thinning was not visible. The GPP_COS during the year 2022 was lower than in 2020, which is clearly not the case in GPP_NEE. Generally, GPP_COS overestimates the GPP compared to GPP_NEE^3,47, but in 2021, GPP_COS deviated from the normal response of the ecosystem and clearly underestimates GPP_NEE.

The disturbance in the boreal forest stand has disrupted the COS budget with an unidentified and potentially significant emission source. In order to estimate the emission in the 2021 summer, the COS uptake by the forest has to be quantified indirectly. The COS uptake by the forest stand can be quantified by inverting the COS-GPP relation (F_COS,GPP)². Assuming that the forest was the major sink, the difference between the F_COS and COS uptake thus gives the COS emissions (equation (11), Supplementary Fig. 8, F_ems:COS,GPP). The mean estimated COS emissions F_ems:COS,GPP thus estimated were −6.0 ± 14.4, 23.3 ± 17.5 and 2.7 ± 14.4 pmol m⁻² s⁻¹, respectively, in 2020, 2021 and 2022 for the June-July months. This partitioning method for COS had its uncertainties associated with using GPP_NEE, LRU and nighttime COS uptake. The emission in 2020 was negative since F_COS,GPP was lower than F_COS, i.e., there was an extra COS sink in 2020 that could not be accounted for by the GPP in that year. Another method to estimate COS uptake is to assume a linear relation between G_C and F_COS and use the slope during 2020 and 2022 to estimate the COS uptake in 2021 F_COS,GC. The difference between F_COS,GC and F_COS (equation (12), Supplementary Fig. 8) would give us an estimate of COS emission (F_ems:COS,GC: −0.8 ± 12.9, 12.2 ± 14.3 and 0.8 ± 12.5 pmol m⁻² s⁻¹ for the years 2020, 2021 and 2022, respectively). F_COS,GC considers canopy-level conductance upscaled from a few tree samples in the forest stand and completely ignores the subcanopy vegetation and soil uptake of COS fluxes (Supplementary Fig. 8). However, F_COS,GC appear closer to F_COS in 2020 and 2022 (Supplementary Table 2). Even though both estimation methods have their own uncertainties, they showed COS emission in 2021 that overwhelmed the average uptake of COS in June–July. The COS atmospheric concentrations (Supplementary Fig. 9) were high during June-July of 2021 (396.6 ± 28.9) compared to the other two years (2020: 281.7 ± 39.4 and 2022: 307.1 ± 37.6). The leaf and soil COS uptake are concentration-dependent³⁹; therefore, the estimated emissions could be underestimated.

Despite the aforementioned uncertainties, the main drivers for F_ems:COS,GPP and F_ems:COS,GC could be clearly identified (Supplementary Fig. 10). The correlation was significant only in the 2021 summer, as the F_ems:COS,GPP and F_ems:COS,GC are close to zero in 2020 and 2022. The univariate analysis against PAR, T_a, VPD, T_s,humus and surface layer soil water content (SWC_humus), shows that T_a and T_s,humus and PAR primarily drive the emissions in both daily and weekly scales. Various combinations of the drivers, avoiding intercorrelation, were used as drivers in the multivariate analysis. A combination of PAR + T_a and PAR + T_s,humus was able to best explain the emissions. It is likely that the emissions were driven by a photothermal component, and we are unsure if it is the temperature that drives the effect of radiation or vice versa.

The nature of the emission sources, whether they are occurring in hotspots or spread evenly through the footprint of the EC measurements, could be inferred from a flux-weighted map of positive F_COS (Supplementary Fig. 11) and weighted mean polar plots (Supplementary Fig. 12). Both analyses point towards the unknown COS source distributed heterogeneously within the footprint area. Thinning in the Hyytiälä forest stand has opened up the canopy and doubled the amount of radiation reaching the forest floor (PAR_ff in pre-thinning years: 100.2 ± 112.9 μmol m⁻²s⁻¹ to 242.8 ± 245.8 μmol m⁻²s⁻¹ in 2020). The forest floor vegetation recovered steadily in 2021 and 2022, as demonstrated by the forest floor GPP⁴⁸. The forest floor (R_ff) and ecosystem (R_eco) respiration fluxes during June were higher in 2021 compared to 2020 and 2022, which was evident from the nighttime NEE fluxes (Supplementary Fig. 13) and modelled respiration (Supplementary Fig. 14). The total respiration for the entire growing season in 2021 was lower than in 2020 or 2022⁴⁸, but the warmer June-July months have seemingly spiked the ecosystem respiration in 2021 (Supplementary Fig. 14). A similar pattern can be found in the forest floor NEE (Supplementary Fig. 14).

The presence of S-containing amino acids (SAAs) in soil was found to be a potential precursor of COS^34,56. These SAAs could be introduced into the soil through the excess litter left after the thinning. The amount of leaf litter of various species (Supplementary Table 4) over the years and the estimate of S content^57,58,59,60 in them can give a rough estimate of available S in the soil. The yearly S-content from litter was 0.41 ± 0.1, 0.46 ± 0.1 and 0.44 ± 0.1 g S m⁻² in 2020, 2021 and 2022, respectively, which may drive the observed COS emissions.

The COS emission during the senescing period of wheat soil was associated with the mobilisation of SAA and its metabolism¹⁴. It should be noted that the surface moss layer in Hyytiälä was considered a COS sink¹⁸, due to the consumption of COS by CA enzymes in mosses. The thinning resulted in ~120 g C m⁻² decomposable fine-cutting residue, including leaves and needles. Aslan et al.⁴⁸ estimated that ~50% of the residue already decomposed in 2020–2022 post-thinning years, assuming an exponential decomposition rate⁶¹. An initial delay in the decomposition is possible, as the cutting residue was relatively intact and was not incorporated into the litter layer, as it was not compacted by snow right after the thinning in winter 2020. The cutting residue was partially degraded over the following winter and was more susceptible to photothermal degradation over the summer of 2021. The decomposition rate depends on environmental variables such as soil temperature, which was above the seasonal averages, especially in the summer of 2021. The combined effect of the delay in the decomposition of pine needles and the higher temperature possibly resulted in elevated decomposition rates in the 2021 summer. Consequently, the elevated decomposition likely contributed to COS emission from the forest floor.

The source of the COS emission observed in summer 2021 can only be speculated at this point. The global COS budget remains unbalanced, and we have demonstrated that the boreal forest could be a source when disturbed. A complete understanding of the COS sources and sinks is necessary for its effective use as a GPP tracer. The continuous observations from the forest floor in a disturbed forest stand could advance our understanding of COS fluxes. Identifying and quantifying this newly observed yet unknown source is required to close the COS budget, reducing a key uncertainty in using COS as a tracer for global photosynthesis. This improvement will help refine carbon cycle models and enhance our ability to monitor and predict ecosystem responses to climate change, with implications for science and policy.

Site description and thinning

The flux measurements from the SMEAR II station in Hyytiälä, Finland (61°51’N, 24°17’E; 181 m a.s.l.) were used in this study. The forest stand is dominated by 50–60-year-old Scots pine (Pinus sylvestris L.) and consists of other tree species, including Norway spruce (Picea abies (L.) Karst.) and deciduous trees (e.g., Betula sp., Populus tremula, and Sorbus aucuparia). Understorey vegetation is characterised by dominant shrub species, notably lingonberry (Vaccinium vitis-idaea L.) and bilberry (Vaccinium myrtillus L.), and prominent mosses such as Schreber’s big red stem moss (Pleurozium schreberi (Brid.) Mitt.)⁶². In 2018, the basal area was 30.7 m² ha⁻¹, basal area weighted mean canopy height was 19.9 m, one-sided leaf area index (LAI) 3.9 m² m⁻², and stem density 1141 stems ha⁻¹. The stand was thinned in two phases, i.e., understorey removal and main thinning in March–April 2019 and January–March 2020 (Supplementary Fig. 15), respectively, reducing the basal area by ~ 40%, and foliage by ~ 45%, stem density by ~ 60%^48,63. Post-thinning measurements indicated a basal area of 18.0 m² ha⁻¹, LAI of 2.1 m² m⁻², and stem density of 465 stems ha⁻¹. The species composition remained unchanged after the thinning, as the thinning ratio was almost uniform. Due to canopy opening, light availability almost doubled at the forest floor, and wind speed and daytime air temperature increased at the trunk space in 2020 compared to pre-thinning years⁴⁸. The daytime duration varied greatly over the months, from ~19.5 h in summer to ~8.3 h in winter. All results were presented in Finnish winter time (UTC + 2). Daytime is defined as the solar elevation angle greater than 0 degrees, and nighttime is less than 0 degrees.

EC and other measurements

F_COS were measured during 2020–2022 at the height of 23 m, which consists of an ultrasonic anemometer (uSonic-3 Class A, METEK Meteorologische Messtechnik GmbH, Elmshorn, Germany) for measuring wind speed in three dimensions and sonic temperature; an Aerodyne quantum cascade laser spectrometer (QCLS; Aerodyne Research Inc., Billerica, MA, USA) for measuring COS, CO₂, carbon monoxide (CO) and water vapour (H₂O) mole fractions. The processing and gap-filling of the COS fluxes were described by Kohonen et al.⁶⁴. The setup was similar to that reported by Vesala et al.⁵, which used the data from 2013–17 from Kohonen, K.-M. & Kooijmans, L.⁶⁵.

The EC data were processed using the EddyUH software⁶⁶ following the procedure recommended by Kohonen et al.⁶⁴ for COS flux processing. Raw data were despiked so that the difference between subsequent 10 Hz data points was a maximum of 200 ppt for the COS mixing ratio and 5 m s⁻¹ for the vertical wind velocity component (w). The coordinate system was set using a 2D-coordinate rotation, and linear detrending was used to separate the time series into mean and fluctuating components. The time lag of COS was determined by maximising the cross-covariance between CO₂ mixing ratio and w. High-frequency spectral losses were corrected using the experimental approach⁶⁷ according to Mammarella et al.⁶⁸ so that the response time of CO₂ was used for COS spectral correction. Low-frequency correction was done according to Rannik and Vesala⁶⁹. Finally, only 60% of all the half-hourly fluxes that met the following criteria were accepted: the number of spikes in the 30-minute COS mixing ratio and w less than 100, the second wind rotation angle less than 10 degrees in absolute value, flux stationarity less than 0.3, kurtosis between 1 and 8, and skewness between −2 and 2 (Flag 0 and 1 in Supplementary Fig. 16).

The net ecosystem exchange of CO₂ was separately measured from above the forest canopy (NEE_eco) and at the forest floor (NEE_ff) using EC. The NEE_eco was measured at a height of 27 m above the ground using a closed-path nondispersive infrared sensor (LI-7200, LI-COR Biosciences, USA) paired with an ultrasonic three-dimensional anemometer (HS-50, Gill Ltd, UK) measuring at 10 Hz. This setup is part of the Integrated Carbon Observation System (ICOS) Ecosystem Station network^70,71 and the measurements are according to the ICOS recommendations⁷². (NEE_ff) was measured at 2.4 m height with an EC system that comprised of a closed-path nondispersive infrared sensor (LI-7200, LI-COR Biosciences, USA) along with a three-dimensional ultrasonic anemometer (model USA-1; METEK GmbH, Elmshorn, Germany). Gas sampling was conducted using an Eaton Synflex tube, with diameters of 6 mm (0.35 m length) and 12.5 mm (1.67 m length), at a flow rate of 12 LPM.

The temperature was measured at two heights, 16.8 m and 33.6 m, using Pt100 sensors inside ventilated radiation shields. Similarly, RH was measured at 16.8 and 35 m heights using Rotronic MP102H RH sensors. The average of the measurements at these two heights was used to represent the conditions at the flux measurement height of 23 m. Photosynthetically active radiation (PAR) was measured at 35 m height using a LI-COR LI-190SZ quantum sensor. T_s was measured at five locations in the following depths: organic layer(soil surface), A1 (2–5 cm), B1 horizon (9–14 cm depth), B2 horizon (22–29 cm) and C horizon (42–58 cm) of the forest floor using Philips KTY81-110 temperature sensors. Soil water content (SWC) was also measured using a Delta-T ML3 sensor at 4 locations at above mentioned depths. The all-sided leaf area index (LAI) was measured for a three-day interval using an array of four under-canopy PAR measurement sensors in four different plots by inverting the Beer-Lambert equation (Supplementary Fig. 4). The average of all the sensors was fitted to a bell curve using the scipy.optimize.curve_fit⁷³ package, and the wintertime LAI was interpolated (October-March).

Sap flow density was measured with constant heat-dissipation probes⁷⁴ in eight Scots pine and four silver birch trees. Sensors were installed at breast height on the northern side of the stem, with the reference and heated needles 10 cm apart. The sensors were covered with a reflective aluminium shelter to prevent temperature artefacts. Data were recorded every minute for trees on an RMD680 channel transmitter (Nokeval, Nokeval Oy, Nokia, Finland). Data were averaged over 15-minute periods, and sap flow density was then calculated according to the original calibration⁷⁴. Maximum temperature difference under zero flow conditions was estimated according to Oishi et al.⁷⁵ with the following thresholds: nighttime conditions (PAR < 50 μmol m⁻² s⁻¹), low vapour pressure deficit (<0.2 kPa) and a stable temperature difference (coefficient of variation < 0.5% over two hours period). When those conditions were not met, the maximum temperature difference under zero-flow conditions was linearly interpolated from neighbouring days.

Sapflow upscaled stomatal conductance

The sap flow measurements from twelve trees (eight Scots pine and four birches) distributed around the forest stand are used to estimate the transpiration and upscale to the ecosystem level. The diameter class technique for upscaling the sap flow observations is utilised, which is based on the assumption that sap flow rates depend only on the diameter at the breast height (DBH) of the tree^76,77,78. The DBH distribution for different species of trees per hectare in the sampled area is given in Supplementary Fig. 17. An exponential function is fitted to the relationship between sap flow rates (Q_tree) of the sampled trees and their DBH^78,79 for Scots pine and birches separately. The coefficients from the exponential function fit are used to obtain the average sap flow for all given tree DBH classes. Further, sap flow for all the DBH classes was aggregated accordingly for sample areas (on a scale of 1 ha) and divided by the average sap flow for sample trees measured within the sample areas. This allowed us to obtain a specific scaling factor (F_s) for sap flow measurements at thinned and control areas, transforming the measurements performed on sample trees to an ecosystem scale. The stand transpiration (TSF) for sample areas was obtained equation (1).

The stomatal conductance upscaled from sapflow was calculated by dividing TSF by VPD (equation (2)).

Upscaling sap flow measurements from tree level to ecosystem level is challenging and brings uncertainty to achieved results^78,80. The comparison of upscaled sap flow measurements, which on a daily or longer time scale represent transpiration with fluxes measured by the EC system, is affected by the footprint dimension and tree distribution within the fetch area. It has also been reported that the uncertainty of upscaled sap flow measurements could also result from the systematic underestimation of sap flux by the sensors (especially thermal dissipation probes, operated according to the constant power principle^81,82, but also for the Tissue Heat Balance method⁷⁸). Upscaling the sap flow measurements is more challenging for mixed forests than monocultures, which also causes uncertainties⁸⁰. Nevertheless, it is a known and widely used method for assessing transpiration at the ecosystem level based on point sap flow measurement in the tree scale.

COS uptake parameterisation

The COS fluxes from the Hyytiälä forest stand were parameterised using the PAR, T_a, VPD, and all-sided LAI from the long time series data shown by Vesala et al.⁵. The LAI was reduced by almost 40% after the thinning in 2020. The parameterisation of COS flux uptake is given as equation (3).

The four functions for the parameterisation and the parameter values were taken from⁵, where FPAR is the function of stomatal response to PAR; FS is the phenology of biochemical reactions; FVPD gives the stomatal regulation, and FLAI is the function of foliage and canopy light penetration (Supplementary Fig. 18). The functions and parameterisation description can be found in Vesala et al.⁵. The parameterised COS fluxes are plotted against the measured EC COS fluxes in Supplementary Fig. 5.

Mechanistic soil COS flux model

We used a depth-resolved diffusion-reaction model (COS Soil Model, COSSM) to simulate soil COS fluxes from meteorological and soil physical inputs³⁰. The model represents both abiotic COS source and enzymatic COS uptake activities and evolves the COS concentration profile prognostically to solve for soil-atmosphere COS fluxes. Here, key improvements have been made for simulating soil COS fluxes at Hyytiälä. Firstly, the model framework has been reimplemented in a differentiable manner using the JAX library in Python⁸³, which facilitates parameter optimisation against soil chamber observations and improves computational efficiency through just-in-time compilation. Secondly, the vertical discretisation of soil layers now follows the Community Land Model version 5 (CLM5)⁸⁴, with 10 layers of increasing thickness down to 1.5 m depth, fewer than the 25 layers in the original model³⁰. The reduced total number of soil layers and tightly packed soil layers near the surface make simulations faster without aliasing the COS concentration profile.

We used chamber measurements of soil COS fluxes at Hyytiälä in 2015¹⁸ to parameterise the model. Between the two soil chambers in 2015, we used data from soil chamber #1 (SC1) as the training set to optimise COS metabolic parameters and data from soil chamber #2 (SC2) as the test set to assess model generalisability. The parameters that were optimised include COS production capacity (V_max,P), temperature sensitivity of COS production as indicated by Q₁₀, carbonic anhydrase (CA) COS uptake capacity (V_max,CA), Gibbs free energy, enthalphy change, and entropy change associated with enzyme activation/deactivation (ΔG_a, ΔH_d, and ΔS_d), and the optimal soil moisture for COS uptake (θ_opt). The Michaelis constant for CA-mediated COS uptake (K_m) is fixed at 0.039 mol m⁻³ following the experimental value for β-CA^7,24. See Supplementary Table 5 for prior and optimised parameter values. The optimised model achieved a root-mean-square error (RMSE) of 1.12 pmol m⁻² s⁻¹ on the training set and an RMSE of 1.53 pmol m⁻² s⁻¹ on the test set for hourly flux measurements and better performance on daily fluxes (Supplementary Fig. 19). The reduced performance on SC2 data is not surprising given that the root-mean-square deviation between quasi-concurrent SC1 and SC2 flux measurements is 1.50 pmol m⁻²s⁻¹, reflecting small-scale heterogeneity.

We then used meteorological and soil physical observations at Hyytiälä during 2013–2022 to drive the optimised model for half-hourly soil COS flux simulations. June–July soil COS fluxes for 2013–2017 and 2020–2022 are summarised in Supplementary Fig. 6.

GPP estimate from COS

The measured EC COS fluxes comprise both vegetation and soil flux components. There were no COS soil chamber measurements during the 2020–2022 period. Hence, average COS soil flux (−2.7 pmol m⁻²s⁻¹) measured by Sun et al.¹⁸ from 2015 in Hyytiälä was used for this analysis. The soil contribution was subtracted from the total EC COS fluxes to estimate the vegetation contribution (F_COS,canopy) and the GPP from COS fluxes (GPP_COS) was estimated using the method described by Sandoval-Soto et al.² (equation (4)).

where ({[C{O}_{2}]}_{a}) and [COS]_a are mixing ratio of CO₂ and COS measured by the QCL spectrometer. The leaf relative uptake (LRU) was estimated as a function of PAR (equation (5))³.

GPP from partitioning of NEE

The GPP_NEE was estimated by traditional NEE flux partitioning (equation (6)) :

where R was estimated from the nighttime NEE fluxes (equation (7)).

where R_c is the respiration at a reference temperature (T = 10 °C) and Q₁₀ is the temperature sensitivity at which the R increases for every 10 °C and T is the driving temperature for the ecosystem. The mean of the soil temperature at 5 cm depth and air temperature at 16.8 m height was found to be a better driving temperature for the entire ecosystem respiration^85,86, in the case of the forest floor respiration, only the soil temperature at 5 cm depth is used. The Q₁₀ was 2.0 for the above-canopy fluxes and 2.5 for the forest floor, while R_c was estimated in a moving time window of 11 days (above-canopy fluxes) or 15 days (forest floor fluxes) and seasonal variability ranged from ~ 0.5 to 2.5 μmol m⁻² s⁻¹ for above canopy and ~ 0.2 to 1.2 μmol m⁻² s⁻¹ for the forest floor. A non-linear regression model for GPP (GPP_NLR) is used when the NEE measurements are not available (equation (8)).

where α and P_max are fitting parameters that were estimated in a moving time window, f(T_a) is an instantaneous temperature response that reduces the GPP to zero near freezing temperatures (equation (9)). The seasonal variability of P_max was 0–25 μmol m⁻² s⁻¹ and 0–5 μmol m⁻² s⁻¹ at the canopy and the forest floor, respectively. α ranged from 0 to 0.04 and from 0 to 0.015 (dimensionless) at the canopy and the forest floor, respectively.

where the inflection point is at T₀ = − 2 °C⁸⁷.

Estimation of apparent COS emissions

The GPP-COS relation described by Sandoval-Soto et al.² can be inverted to estimate the COS uptake by the ecosystem (equation (10)). The GPP_NEE was found to follow the natural restoration of the forest stand post-thinning⁴⁸ (Fig. 4). If there were no emissions in the canopy, the F_COS would have also shown increased uptake as the forest recovered. Therefore, GPP_NEE is input in the GPP-COS equation to estimate the COS uptake of the ecosystem (F_COS,GPP) as if the emissions were absent. This method assumes that COS uptake is happening only during photosynthesis and does not consider the COS uptake during the nighttime.

The COS uptake thus estimated was used to find the emission as given in equation (11).

The uptake of COS during the day and night is a direct function of the ecosystem stomatal conductance (G_C)¹¹, which can also be used to estimate the COS uptake in the forest stand. Assuming the COS emissions were negligible in the summer of 2020 and 2022, the linear relation between G_C and F_COS was estimated (F_COS = −0.21× G_C – 12.4) and was used to estimate the canopy uptake of COS (F_COS,GC) in the 2021 summer (Supplementary Table 2). The COS emission was estimated as given in equation (12).