Primary productivity of managed and pristine forests in Sweden

Land use is affecting 70% of global ecosystems and their functioning. Forest management is a regionally dominant land use and affects forest ecosystems by changing both structure and functioning, but its impact on primary productivity is not well known. Here we investigated the effect of forest management on primary productivity by comparing managed secondary forests with relatively pristine unmanaged primary forests in Sweden. As proxy for primary productivity we used the satellite-based vegetation index NIRv which has been shown to be closely and linearly related to primary productivity. We produced a digital map of 390 primary forests across Sweden, and extracted NIRv over these and surrounding secondary forests forming spatially proximate pairs. By comparing the primary and secondary forests NIRv in the pairs we found that secondary forests on average show higher NIRv, but the highest values were found in primary forests. The difference in NIRv between pairs is related to their difference in mean stand age, and at equal stand age the NIRv of primary forests is higher than in their paired secondary forests. Overall, management leads to increased NIRv through regeneration of forests stands that reduce their mean age. However, primary forests show higher NIRv when controlling for age, despite being found on higher altitudes and on steeper slopes with lower soil moisture, which suggests that forest management other than regeneration is not increasing primary productivity of Swedish forests.


Introduction
A majority of the global land surface is managed as croplands, pastures and production forest (Erb et al 2007). Forestry is a regionally dominant land use and is the main land use type in Sweden, covering~60% of the total land area (SLU 2019). The importance of understanding how forestry affects the ecosystem carbon (C) cycle is highlighted by future scenarios designed to limit global warming by large expansion of bioenergy production (Smith et al 2016, Rogelj et al 2018. However, it remains unclear how forest management affects C uptake through gross primary productivity (GPP; the sum of photosynthesis in the ecosystem), which together with carbon use efficiency (CUE) and the turnover time of vegetation determines C storage in vegetation.
In secondary forests, harvesting and re-planting leads to reduced stand age (Schulze et al 2012) and a strong relationship between productivity and forest age has been documented where younger forests are more productive than older forests (Ryan et al 1997, Pregitzer and Euskirchen 2004, Besnard et al 2018. Selection of fast growing genotypes, fertilization, drainage and thinning are common practices in forest management (Kurz et al 1998). While some aspects of management may increase GPP, the management is not designed to increase GPP. It is instead designed to optimize the production of harvestable wood and in some cases various aspects of timber quality, which is reflected in the increased allocation of accrued C to above ground compartments (Noormets et al 2015). The management induced shift in allocation and species composition may also result in higher CUE (Campioli et al 2015), which implies a higher net primary productivity (NPP) to GPP ratio.
Primary forests are an analogue to which managed ecosystems can be compared. They are not managed for wood production and are in a relatively pristine condition, characterized by stands originating from natural succession, with a relatively large share of old trees and dead wood. Primary forests have not been used, often because they are found in rugged, high altitude and inaccessible terrain (Sabatini et al 2018). A recent study analyzing the difference between managed and unmanaged forest globally found little difference in GPP between the management types (Noormets et al 2015), but it is unclear how comparable the primary and secondary forests are in terms of their location and environment. A large-scale assessment of the impact of forest management on GPP is generally lacking, and no such study has been conducted over Sweden to our knowledge.
Officially, Sweden has the largest area of primary forest in Europe (excluding the Russian Federation) making it a suitable study region that may be representative for the much larger global extent of boreal forests. Primary forests are found from the southern temperate zone to the boreal and colder northern regions. However, no official national map or other digital product of the spatial location of primary forests is available to our knowledge.
This study aimed to assess the impact of forest management on GPP. By focusing on GPP we did not analyze productivity of biomass or merchantable wood, but rather the total ecosystem GPP which includes understory and ground vegetation. To study how GPP differ between primary and secondary forests we first created a digital a map of Swedish primary forests and located secondary forests near each primary forest for a spatially paired analysis of primary productivity and forest characteristics. Over primary and secondary forests we analyzed the satellite-based vegetation index Near Infrared Reflectance of terrestrial vegetation (NIRv; Badgley et al 2017). NIRv has been found to be closely related to GPP and show less saturation at higher vegetation cover than other vegetation indices such as NDVI (Badgley et al 2017, Huang et al 2019. In the analyses we included stand-age, altitude, steepness, and soil moisture to isolate the effects of management from differences in environmental conditions.

Study area
Sweden has a total area of productive forest of 235 500 km 2 which covers about 60% of the country (Naturvårdsverket 1982, SLU 2019). The country includes temperate and boreal vegetation zones, but is mainly boreal. The majority of the forested area is planted and managed secondary forest dominated by Scots pine (Pinus sylvestris) and Norway spruce (Picea abies) (Fridman 2000, SLU 2019). More intensive forest management and the utilization of its resource may have started 300 years ago and is today resulting in a landscape dominated by monocultures of even aged stands (Esseen et al 1992). In the last decades the Swedish forest area has remained relatively stable, as a result of conservation policies and little expansion of other land uses. The annual change in forest cover is + 0.01% (1990-2010) (SLU 2019).

Producing a digital primary forest map of Sweden
There are many definitions of primary forest ecosystems. Definitions may be based on stand age, minimum area and known current and historical management (IPCC 2003, Buchwald 2005, FAO 2006, 2010, 2011, Thompson et al 2013. In this study, primary forests are defined as naturally regenerated forests of predominantly native species that are currently (and historically) insignificantly disturbed by human activities such as felling or thinning of trees (Naturvårdsverket 1982). A main indicator of low human disturbance is the existence of larger quantities of dead wood, which implies that the trees have not been harvested and removed but left to die and decompose in the ecosystem. The existence of dead wood is in general a clear and strong indicator of protection and higher stand age of forest in Sweden (Fridman 2000).
Here we digitized an extensive nationwide inventory of primary forest in Sweden consisting of a set of paper reports from the early 1980s and included maps of the forest boundaries (Naturvårdsverket 1982). The reports contain both publicly and privatelyowned forests and were visited by the coordinating working groups. Each forest is described based on known history, continuity, forest type and approximate stand age, as well as the presence of deadwood. From this report series all 423 sites were collected, but 32 sites were excluded in the digitization process due to small areal extent (<0.01 km 2 ). Forests that have experienced known felling or thinning over more than 5% of their total area since the original inventory in 1982 were excluded from the map (supplementary table S1). The database on known felling and thinning used here includes more than 932 000 spatial polygons in vector format of individual thinning and felling events. It is currently compulsory to report the spatial extent of felling and thinning when an area above 0.5 hectares is affected. Although the earliest recorded fellings are from the early 1990s, the database includes few recoded fellings until the late 1990s. This implies that the database on fellings is incomplete between the time of the inventory in the early 1980s until the late 1990s. To account for this uncertainty, we set the relatively low exclusion limit at 5% area felled. 5% felling indicate that the forest has not been preserved, and it may indicate that more felling has occurred before the start of the database. Known felling of less than 5% on the other hand often occur on the edges of forests, which suggests that the boundaries of some of the primary forests are imprecise, and these forests were therefore not excluded. Water and other non-forested areas such as bare rock and wetlands were also excluded from further analysis (supplementary table S1).

Identifying spatially proximate secondary forests
For each primary forest, secondary forests within a 15 km buffer zone from the primary forest edge were identified from the Swedish land cover map (supplementary table S1, see illustration in figure 2(b)). We excluded areas unlikely to be managed forests with several datasets indicating land use and cover, these include areas under various protection statuses and areas that are mapped as having high conservational value. Water bodies were excluded from both the primary and secondary forests, a full list of exclusions can be found in the supplementary table 1. A 15 km buffer was selected as it represents a balance of the amount of secondary forest included and distance from the primary forest. Initial testing of varying buffer sizes around a subset of the primary forests indicated that the analysis was not very sensitive to buffer distance within reasonable distances (5-50 km, not shown here), but that the smallest buffer sizes in some cases included low amounts of secondary forests while the larger buffer sizes may decrease the similarity in environmental conditions between the primary and secondary forests in the pairs.

The satellite-based proxy of primary productivity
We calculated NIRv ((NDVI-0.08) x NIR)/1000 based on surface Reflectance (Tier 1) data recorded by Landsat 8 Operational Land Imager (OLI) (Badgley et al 2017(Badgley et al , 2019, which provides data at a spatial resolution of 30 × 30 m. For this study, we created a 5-year time-series (2014-2018) covering Sweden in Google Earth Engine. All scenes were masked for clouds, cloud shadows, and snow using the included quality flags of the Landsat dataset (LANDSAT/LC08/C01/T1_SR) (Masek et al 2013). For each pixel the per grid cell temporal 95th percentile of the 2014-2018 time series was extracted.

Ancillary data
To characterize landscape and forest characteristics which may influence GPP we include datasets of stand age, soil moisture index (SMI) and topography. Stand age was derived from the SLU forest map at 25 × 25 m spatial resolution, originating from forest inventories and satellite data (SLU 2010). Elevation and slope were based on the Swedish national elevation map at 50 m spatial resolution (Lantmäteriet 2016). The soil moisture index (SMI) is a model product based on depth to groundwater inventories and the topographic wetness index at 10 × 10 m spatial resolution (Naturvårdsverket 2019).

Statistical analysis
Pairwise analysis of NIRv values for each primary forest site and the adjacent secondary forest buffer zone was performed. From the 95th temporal percentile composite, the spatial median was extracted for each primary forest site and each buffer zone of secondary forests. To test whether NIRv differed between primary and secondary forests, we performed a paired sampled t-test. To investigate forest and environmental characteristics affecting differences in NIRv between the primary and the secondary forests, we studied the variability of the ancillary variables using boxplots. A multiple linear regression was also performed between the paired differences in NIRv and stand age, altitude, topographical slope and SMI to test if environmental characteristics could affect differences in NIRv between primary and secondary forests.

Digital Primary forest map of Sweden
The digitized and updated primary forest map of Sweden contains 17 629 km 2 of primary ecosystems located in 391 sites (mean ± one standard deviation, 45 ± 203 km 2 ) (figure 1). Of the 391 sites, 33 were excluded due to later felling or thinning over more than 5% of the total non-water area. The remaining 358 sites have a total area of 16 565 km2 (mean size 46 ± 199 km 2 ) out of which 9 926 km 2 is covered by mature or young forests (after excluding water bodies and non-forested land covers). These remaining forest areas were used in the subsequent analysis. Most of the primary forests are small (mean size 28 ± 107 km 2 ), with the exception of the sites in the northwest mountainous region of Sweden (figure 1).

Forest primary productivity
The pairwise comparison of primary versus secondary forests NIRv shows that secondary forests generally have higher NIRv values (figure 2). On average, NIRv in primary forests is 0.217 ± 0.054, and in secondary forests it is 0.242 ± 0.032, indicating that primary productivity is higher in the secondary forests (t-test, p < 0.0001). However, primary forests show a higher variability in NIRv values, and the highest values are found in primary forests.

NIRv and environmental characteristics
Across all forests, primary forests have significantly higher stand age (p < 0.0001), are found on higher altitudes (p < 0.05), steeper slopes (p < 0.0001), and have lower estimated soil moisture (p < 0.0001), in comparison to their surrounding secondary forests ( figure 3(a)). We also evaluated forest pairs above and below the 1:1 line in figure 2. Forest pairs above the 1:1 line are characterized by higher NIRv in the secondary forest in comparison to the primary forest of the pair, the opposite is true for pairs below the 1:1 line. In general, pairs above the 1:1 line show larger stand age and environmental differences between the secondary and primary forests of the pairs, whereas pairs below the 1:1 line show paired differences closer to zero ( figure 3(b)). This indicates that environmental differences and stand age may explain the observed higher mean NIRv value in secondary forests.
We also analyzed the characteristics of the forest pairs that constitute the 10 largest outliers on each side of the 1:1 line ( figure 3(c)). For outliers above the 1:1 line (secondary forest NIRv ≫ primary forest NIRv) secondary forests have significantly younger stands (p < 0.0001), are found on less steep slopes (p < 0.05), and have higher estimated soil moisture (p < 0.01). Outliers below the 1:1 line (primary forest NIRv ≫ secondary forest NIRv) show a similar pattern, but the differences are smaller and only altitude difference show a significant difference from zero (p < 0.05).
Overall, the analysis of environmental characteristics suggests that the higher mean NIRv of secondary forests may be partly explained by preferential environmental conditions and younger stand ages. To analyze this in more detail we performed a multiple ordinary least-square linear regression with paired NIRv difference as predictand and paired difference in stand age, altitude, slope and SMI as predictors. Analysis of partial relationships to each predictor indicate that stand age difference exert a strong linear relationship on NIRv difference, while the other variables showed no strong relationship in the full sample of pairs (figure S1 (available online at stacks.iop.org/ERL/15/094067/mmedia)). We therefore focused the subsequent analysis on stand age.  The intercept of the linear relationship between stand age difference and NIRv difference is weakly but significantly positive (0.0074, 95% confidence interval 0.0021 to 0.0126) ( figure 4(a)). This indicates that, when using a linear relationship of the full sample, NIRv of primary forests is higher than in the secondary forests (positive NIRv difference) at equal stand age. Given the relatively small positive difference in NIRv, we further investigate the robustness of the positive difference in NIRv at similar stand age by . NIRv difference and stand age differences. (a) NIRv difference and stand age difference for the full sample (n = 346). Some forests are missing due to the incomplete coverage of the forest age map. (b) Mean NIRv difference of stand age bins ranging from ± 2.5 years to ± 15 years. Each boxplot represents the mean NIRv values of 100 randomly drawn samples. The random sampling ensured that an equal number of NIRv values had negative and positive stand age difference.
calculating the mean NIRv difference in samples with increasing positive and negative stand age difference from zero. In other words, we calculate the mean NIRv of forest pairs within bins of increasing width (in years) centered around zero stand age difference. Because there are more pairs where the primary forest is older than the secondary forest, we randomly selected pairs from the positive stand age difference side of the bin (secondary forests stand age > primary forest stand age) so that the same number of pairs were selected with negative and positive stand age difference (all pairs with negative stand age difference within the bin were always used). This procedure was repeated 100 times, and visualized as boxplots ( figure  4(b)). The resulting means are dominantly positive for all bin sizes, and the analysis supports that NIRv of primary forests is higher than the NIRv of secondary forests when controlling for stand age. The mean age of the stands within the ± 15 years bin range from 33 to 130 years, which implies that the positive NIRv difference is true for a large part of the stand age distribution of the sample (mean ± standard deviation: 76 ± 23 years, min: 26 years, max: 148 years).
We also evaluated the potential effect of forest edge effects by comparing primary forest area to the residuals of the age difference to NIRv difference linear relationship (figure 4(a)) by applying a linear regression. A log transformation was applied to the primary forest area because the forest area distribution is skewed with a majority of small forests and a few large forests. The analysis reveals that forest age is weakly related to the residuals (R 2 = 0.002, p < 0.01) and show a positive slope (figure S2). The positive slope implies that the positive difference in NIRv between primary and secondary forests increases with primary forest size. Under the assumption that primary forest area is related to human impacts through edge effects, a large area would indicate lower human impact. Using this assumption, we can extract the linear regression slope value for the largest forest and add that value to the stand age controlled NIRv difference which would increase the positive difference in NIRv between primary and secondary forests by about 0.013 (new zero stand age difference intercept: 0.0206, 95% confidence interval 0.0093 to 0.0319) ( figure  S3 and S4). If instead the mean primary forest area was used, the effect would be smaller but still positive (0.0055, 95% confidence interval 0.0063 to 0.0198). Overall, the analysis suggests that there may be a negative impact of edge effects on the NIRv of primary forests, but also that the overall conclusion of the study would not be affected by the inclusion or exclusion of edge effects in the analysis. However, it is unclear if the observed effect is caused by edge effects alone or if the effects are indicative of human impacts, in part because the largest forests are found in areas of lower population density with lower accessibility.

Discussion and conclusion
In this study we analyzed the impact of forest management on GPP by comparing NIRv between primary and adjacent secondary forests in Sweden. While NIRv is not a direct measure of GPP it is strongly related to GPP (Badgley et al 2017, Huang et al 2019. Remote sensing products have the benefit of high spatial and temporal resolution which makes it possible to perform analysis at an ecosystem scale across Sweden. However, satellite measurements are affected by the atmosphere (atmospheric depth, water vapor, aerosols, haze, and cloud shadows), sun-sensor geometry, slope of the terrain (generating reflections from adjacent pixels and shadowing effects), and illumination variations (Jensen 2007). These uncertainties generally have a negative impact on vegetation indices (Jönsson andEklundh 2004, Tagesson et al 2015). We therefor extracted the 95th temporal percentile of each grid cell because it is likely that the higher percentiles represent cloud free conditions with less atmospheric influence. The pairwise analysis in turn implies that we do not rely on the absolute NIRv values, but rather on local differences between primary and secondary forest where environmental conditions should be more similar than between forests types distributed across Sweden.
A strength of the novel primary forests map constructed in this study is the consistency of the definition of primary forests in the original inventory. Variations in the definition could be assumed to be smaller than if we combined maps from multiple separate inventories. The map is different from existing national and global products describing protection status (UNEP-WCMC, and IUCN 2020), because it is a map of forests in pristine to near pristine condition. These pristine forests are not necessarily protected, and protected forests are not necessarily pristine. They may be protected for reasons other than low historical human impact (protection of cultural landscapes, high biodiversity or other natural and cultural values), and information on historical human land use is not included in maps of protected areas. In contrast to the global Intact Forest Landscapes (IFL) map that is based on global datasets and remote sensing (Potapov et al 2017), our map is based on field inventories conducted by the local authorities. The IFL map also apply a minimum forest area of 500 km2, and therefore only include 13 forests located within or partly within Sweden. As a result of the area requirement, all forests in the IFL are located in the mountainous region of Sweden. The forests in our map are well distributed across Sweden, and they span a large range of climates and potential productivity. While some forests were excluded in the digitalization process due to reported fellings, no newly discovered primary forests were added, and it is unclear how complete the original inventory was. Our mapping was therefore not designed to evaluate the total area of primary forests in Sweden, more research is needed to pursue such efforts. Still, our map is the most complete, consistent and readily available digital map currently known to us.
The reasons why the primary forests have not been affected by forest management may vary; some forests have been preserved for recreation and hunting, others have likely remained unused due to low productivity or low accessibility. The pairwise analysis was designed in part to address this potential bias in environmental conditions, under the assumption that the differences in environmental conditions between a primary and a secondary forests are reduced when comparing spatially proximate forests.
On average, secondary forests have higher NIRv, and likely higher GPP, than primary forests in Sweden. The lower stand age in secondary forests is the main explanation, and at equal stand age our analysis suggests that NIRv in primary forests is marginally but significantly higher than in secondary forests. The analysis of primary forest area strengthens the robustness of this conclusion. Potential reasons for the higher age-controlled NIRv in primary forests include, uneven tree age distributions and height structure, more understory vegetation, and higher species diversity. A study from eastern United States found fourfold higher ground cover in primary forests as compared to secondary forests (D'Amato et al 2009). In Swedish boreal forests understory vegetation increase with forest age (Nilsson and Wardle 2005). Understory vegetation may be responsible for a large share of total ecosystem productivity (Harden et al 1997), accounting for up to 39% of gross primary productivity, with a mean of 14%, across 10 eddy-flux tower sites globally (Misson et al 2007), and mosses alone has been found to account for 50% of ecosystem C uptake in a black spruce forest in Canada (Goulden and Crill 1997). Mosses and lichens have also been shown to fix nitrogen (Crittenden andKershaw 1978, DeLuca et al 2002), and a large share of mosses and lichens could therefore potentially increase primary productivity by their direct contribution, and indirectly by fixing nitrogen that is made available to other species in the ecosystem. Primary forests in Sweden are dominated by the same species that are generally planted in the secondary forests, Scots pine and Norway spruce, but it is also likely that they have a larger share of broadleaved species (Fridman 2000). Together with uneven stand age and tree height, these factors could result in a higher NIRv and potentially a higher GPP. Future research and field inventories could be designed to answer these questions by collecting and analyzing data on potential contributing factors.
The study presented here aimed to compare the influx of C through GPP to primary and secondary forests. Together with downstream ecosystem processes GPP is key to understand how management affects the C cycle and C storage. Other processes of major importance include the efficiency of the conversion of accrued C to biomass (CUE and biomass production efficiency), and the turnover of the C in the system. Studies investigating CUE have found higher CUE efficiency in managed and fertile forests (Vicca et al 2012, Campioli et al 2015, suggesting that a secondary forest with a similar GPP as a primary forest could still have a higher NPP and biomass growth. However, it is not known if CUE differs between the paired and spatially proximate primary and secondary forests analyzed here, and how differences in CUE may affect the storage of C in ecosystem compartments other than live vegetation. Investigation of CUE in these ecosystems would be an interesting future research topic that together with estimates of vegetation turnover time would help in disentangling how forest management affects vegetation C storage. The aim of this study was to assess the impact of forest management on GPP. We produced a novel map of primary forests and analyzed the difference in a vegetation index related to GPP between spatially proximate pairs of primary and secondary forests. Our results suggest that forest management increase GPP by active regeneration of forest stands. This stands in contrast to a global (not spatially paired) assessment indicating little to no increase in GPP due to forest management (Noormets et al 2015). At the same time our analysis indicates that at equal stand age, primary forests have a marginally higher GPP. This suggests that forest management does not increase the potential GPP, and that the combined effects of forest managements other than regeneration rather decrease GPP in Sweden.