Developing a spatially explicit modelling and evaluation framework for integrated carbon sequestration and biodiversity conservation: Application in southern Finland

• A modelling framework allowing spatial anddynamic analysis of carbonprocesses and biodiversity aspects was developed. • Modellingof forestrymeasures andemission reductions showed that carbon neutrality could be achieved in the area


H I G H L I G H T S
• A modelling framework allowing spatial and dynamic analysis of carbon processes and biodiversity aspects was developed. • Modelling of forestry measures and emission reductions showed that carbon neutrality could be achieved in the area by 2030. • Application of space and airborne measurements for mapping and monitoring biodiversity and ecosystem processes is described. • Optimal allocation of set-aside areas for conservation would contribute to preserving both biodiversity and carbon values. • Biodiversity gain in the area can be increased without a loss of carbon-related benefits.

G R A P H I C A L A B S T R A C T
a b s t r a c t a r t i c l e i n f o

Introduction
The recent IPBES (Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services) assessment shows that humans are increasing pressure on biodiversity (BD) and ecosystem services (ES) at a truly planetary scale (Diaz et al., 2019). The incipient sixth mass extinction has erased over 300 mammal species and, with them, more than 2.5 billion years of unique evolutionary history . Accelerating anthropogenic land use and greenhouse gas (GHG) emissions -induced climate change is projected to cause further significant deterioration and loss of natural ecosystems and BD, as well as negatively impact regional carbon (C) budgets by depleting the C sequestration and storage capacity of ecosystems (Le Quéré et al., 2018;IPCC, 2019;Overland et al., 2019). Climate change will also severely impact global climatic conditions needed for human thriving (Xu et al., 2020). The challenges posed by climate change, BD loss, and harmful land use are thus deeply interconnected problems (Diaz et al., 2019).
It is increasingly understood that effective combatting of, and adaptation to, these detrimental trends requires development of regional-scale multidisciplinary approaches. In such approaches, spatially explicit data on valuable aspects of BD, and factors determining the regional C budgets and the spatial variation in C processes, can be jointly assessed with state-of-the-art methods. Due to the linkages between BD and C processes, and the limited resources for conservation, cost-effective land use planning also requires methods that can optimally prioritize sites for the co-management of BD and C values. An integrative approach based on multiple data sources, including robust modelling tools to fill data gaps and rapidly updated and ecologically meaningful Earth Observation (EO) applications, is thus needed (Forsius et al., 2013;Kujala et al., 2018a;Holmberg et al., 2019;Buotte et al., 2020;Heikkinen et al., 2020). Such multi-source information enables detection of sites with maximal potential for current and future BD conservation and C sequestration/storage.
In this context, a good understanding of the regional C budget, the magnitude of C storage, and its vulnerability to change, is vitally important (Nabuurs et al., 2015;Bustamante et al., 2016;Akujärvi et al., 2019;Hutchins et al., 2020). Estimation of complete regional C budgets, consistent with national accounting rules, for different land use classes is an integral step to predict responses and potential feedbacks to a changing climate regime, and for the design of different management and mitigation actions/measures (Buffam et al., 2011;Morecroft et al., 2019). Globally, the land sector (agriculture, forestry, and other land uses) is responsible for 10-12 GtCO 2 eq a −1 (carbon dioxide equivalent), about 25% of net anthropogenic greenhouse gas (GHG) emissions, with approximately half from agriculture and half from land use, land-use change and forestry (LULUCF) (Le Quéré et al., 2018;Roe et al., 2019). LULUCF emissions represent the net balance between emissions from soils, land-use change and C sequestration from the regeneration of vegetation and soils. Transforming the land sector and deploying mitigation measures in agriculture, forestry, wetlands and bioenergy could feasibly and sustainably contribute 15 GtCO 2 eq per year (about 30%) of the global mitigation needed in 2050 to deliver on the 1.5°C target (Roe et al., 2019). Recently, governments have also introduced concepts such as 'C neutrality', and 'climate neutrality' where the aim is to achieve a balance between sources and sinks of CO 2 or GHGs by a target year (e.g. EU Commission, 2018). Also, the term 'net zero emission' is increasingly used to describe a broader commitment to climate action, meaning that any GHG emissions are balanced by absorbing an equivalent amount from the atmosphere by natural or technological options (e.g. Rogelj et al., 2019).
There are great expectations on EO systems providing standardized, spatially complete, and cost-effective information for detection, quantification, and forecasting of BD and C processes at large scale (Lausch et al., 2016;Vihervaara et al., 2017). New concepts such as Essential Biodiversity Variables (EBVs) have stimulated progress to unify BD and ES monitoring globally (Pereira et al., 2013;Pettorelli et al., 2016). EO techniques are also in the core in EU policy processes for mapping the landuse changes. Yet, the potential of EO data in BD conservation and other ecological research remains underutilized (Pettorelli et al., 2016). Current EO-based systems provide accurate information on, e.g., forest resources in terms of growing stock, annual growth and distribution of main tree species (Tomppo et al., 2008). However, this information is designed mainly for the needs of forest sector, which greatly affects its usability in other fields. The estimates rely on common features of forest environment that are rather evenly dispersed within the landscape. Focusing on the common forest metrics, the current methods provide reliable data on abundant species, but are unable to provide accurate information on rare species which is essential for BD assessment.
As intensive species-level mapping of BD is infeasible at large scale, various surrogates are used when assessing the state of BD. For example, in boreal forests the amount of dead wood and presence of keystone species have shown to indicate the overall level of BD (Siitonen, 2001;Kouki et al., 2004). The surrogates are typically rare and unevenly dispersed in forest landscape, and mapping their distribution introduces special challenges for EO-based systems. Spatial resolution of a sensor must be detailed enough to capture small details like individual trees. On the other hand, spectral resolution must be sufficient for classifying individual species by spectral reflectance. Finally, temporal resolution needs to be suitable for monitoring the surrogates over time. Thus, new thinking and technologies are needed to capture and detect different features and monitor their changes to support land use and management planning at landscape-level. For example, new EO methods can provide detailed regional information on boreal keystone species such as European aspen (Populus tremula L.) trees Viinikka et al., 2020). Furthermore, biophysical and biochemical properties, and genetic diversity of trees can be studied with recent EO methods (Madritch et al., 2014).
A key challenge for spatial planning is to find a balance between the use of land resources, BD conservation, actual and potential future C sequestration/storage in the ecosystems, and other environmental impacts. A solution to this task is to apply integrated adaptive land use and conservation planning methods which can support sustainable resource management, and lead to win-win situations regarding both conservation and mitigation. Such methods use mathematical formulations and optimization to analyse large amounts of spatial data, and to identify sets of locations that jointly meet the different objectives (Thomas et al., 2013;Moilanen et al., 2014;Kujala et al., 2018a;Buotte et al., 2020;Reside et al., 2020). By identifying areas with high present and future BD values, C storage and C sequestration rates, forestry and other management actions can be planned around these areas to minimize negative impacts on the regional BD and to maximize climate change mitigation. The BD and C storage and sequestration policies can be implemented by a combination of anthropogenic emission mitigation actions, municipal and regional land-use planning, and voluntary actions by private land-owners. The latter can be incentivised by developing economic instruments for joint BD protection and C-storage/ sequestration (Kangas and Ollikainen, 2021;Kosenius, 2021).
The aim of this study is to (i) describe the development of a quantitative, spatially explicit modelling and evaluation framework that integrates results on C storage/sequestration and C neutrality in the forested landscape, mapping of BD values, and joint spatial planning and prioritization of BD and C benefits; (ii) advance the use of diverse EO techniques in producing improved spatial datasets for assessing and monitoring the state of BD; (iii) give an overview of the background assumptions, methodology, data derivation and calculation and mapping of the selected ecosystem variables and BD indicators; (iv) demonstrate the approach and methods using data from an intensively studied large region in southern Finland (Kokemäenjoki river basin), and (v) synthesize key results and discuss uncertainties as well as needs for further developments.
We focus on the forested ecosystems because they cover the largest part of our study area and have the highest potential for C sequestration and management. We use spatio-temporal C-and GHG-modelling and future scenarios based on formal national forestry policies and climate change mitigation strategies for anthropogenic emissions. Highresolution EO data was collected with a particular aim to improve detection of the keystone species European aspen and explore methods for upscaling this information for large-scale spatial prioritization and monitoring.

Study area
The Kokemäenjoki river basin (27,024 km 2 ), located in southwestern Finland, has mixed land-use with forest and seminatural areas covering about 68% of the area (Fig. 1, Table 1). The mean and maximum elevation for the study area is 118 m and 254 m, respectively. The geomorphology of the region has mainly been formed by the latest Ice Age processes. Superficial deposits are mainly composed of gravelly and sandy till (56%), clay and silt deposits (17%) together with pre-Quaternary bedrock exposures, peat and glaciofluvial deposits (GTK, 2014).
The dominant tree species are Norway spruce (Picea abies), Scots pine (Pinus sylvestris), silver birch (Betula pendula) and downy birch (Betula pubescens). Approximately 83% of forest land is on mineral soil, 14% on peatland, and 3% on rocky soil (SYKE, 2019a). Protected areas (PAs) cover 1.8% of the region. Lakes and other water bodies form an important part of the landscape (ca. 11% of total area). The second largest city region of Finland, Tampere, with ca. 330,000 inhabitants is located within the area.
There are two LTER (Long-Term Ecosystem Research) heavily instrumented research stations located within the area, and data from these stations (e.g. C flux data from eddy covariance tower of Hyytiälä SMEAR II station) has been used for model calibration/validation (Fig. 1). The detailed EO-based studies were conducted at Evo test site, which belongs to Lammi LTER area.

Conceptual framework for modelling and evaluation
Our integrated conceptual framework connects different databases, GIS, EO data sets, models and tools, scenarios and socio-economic evaluations, and facilitates the flow of data between different sections (Fig. 2). The methodology, model systems and databases for each section are summarised below (Sections 2.3.1-2.3.4). A large effort was devoted to collect the required input and background information by intensive field studies and from the research sites in the study region ( Fig. 1), national and regional databases, and literature sources (e.g. for emissions coefficients and model parameters). The dynamic model systems were calibrated using site-specific data and other data sources, and the impacts of the different scenarios on the C budget/processes assessed.
We focused here on the methodological development and integrated application on four main sectors: i) scenario-based modelling of C processes/neutrality, ii) derivation and mapping of BD data, iii) integrated land-use and conservation planning based on the derived spatial BD and C data, and iv) the use of EO-based technologies for enhanced data derivation, assessment and monitoring of changes (Fig. 2). The fifth sector, development of economic incentives for land owners is described in detail elsewhere (Kangas and Ollikainen, 2021;Kosenius, 2021). Utilizing our integrated conceptual framework, we show how a systematic approach linked with spatial and quantitative multi-source information can assist in finding optimal land management solutions for climate change mitigation and conservation planning. Calculation of present-day C pools and fluxes for the different landuse classes of the study area is described in detail in Holmberg et al. (2021). Here, we focus on describing the scenario-based modelling for assessing C and GHG fluxes and the concept of C neutrality. The net emissions of C from the forested areas in Kokemäenjoki basin were based on the results of the forest growth and gas exchange model PREBAS, which combines a process-based forest growth model (Valentine and Mäkelä, 2005), and a daily canopy gas exchange model (Peltoniemi et al., 2015). The model has been calibrated using information from Nordic eddy covariance sites and Finnish growth experiments (Minunno et al., 2016. PREBAS was used to simulate forest C balance, growth and management in the study area. The initial state of forests in forest C balance simulations was determined based on MS-NFI (Multi Source -National Forest Inventory) maps, which describe the forest parameters in the form of thematic maps across Finland at 16 × 16 m resolution. MS-NFI maps are developed by combining the information from NFI field measurements, satellite imagery and digital map data (Tomppo et al., 2008). To decrease computational load, similar forest areas were segmented, with the size of segmented area varying between 256 and 505,008 m 2 and a mean size of 3926 m 2 . PREBAS was used to calculate the net ecosystem exchange (NEE, gC m −2 a −1 ) for current climate conditions, and the amount of harvested biomass (gC m −2 a −1 ). The harvested amount for each simulation year was specified by realized removals statistics. Removals data defined the total amount of harvests in the study area. Regarding forests growing on mineral soils, the PREBAS model is linked to the soil C model YASSO07 (Tuomi et al., 2009), which estimated soil respiration. For forests on drained peatland soils, soil respiration estimates were based on measured soil respiration, which includes both peat decomposition and litter decomposition (or accumulation) (Minkkinen et al., 2018). Total NEE was calculated as a sum of stand and soil fluxes, and it acquired positive values when the C flux from the decomposition of soil organic matter was larger than the assimilation of C into growing vegetation, and negative when the assimilation of C into growing vegetation was larger than the C flux from decomposing soil organic matter. The net emissions of forests were calculated by adding the NEE to the amount of harvested biomass. Net emissions are positive when the forest acts as a source of CO 2 to the atmosphere, and negative when the forest is a sink.
The forest harvesting scenarios and their assumptions were based on the national climate and energy strategy (Koljonen et al., 2017) and forestry planning, where total annual removal of biomass in the three main scenarios at the country-scale amount to: Low ca. 40 Mm 3 a −1 , Policy ca. 80 Mm 3 a −1 and Maximum Sustainable Cuttings (MaxSust) ca. 85 Mm 3 a −1 . The scenarios and related forest modelling have been described in Kalliokoski et al. (2019). The PREBAS scenario results for the Kokemäenjoki region were extracted from these national-scale simulations and aggregated using GIS-techniques.
Emissions from other ecosystems (arable land, surface waters and undrained mires) were also estimated although we here focus on forests which are by far the most important ecosystem for C pools and fluxes. A detailed description of methodology regarding the processes and net emissions from other land-use classes is given in Holmberg et al. (2021; see also Supplementary information S.1.1).
The anthropogenic CO 2 , CH 4 and N 2 O emissions of the study area under current and future conditions were calculated with the Finnish Regional Emission Scenario (FRES) model (Karvosenoja, 2008). The emissions are given separately for point sources on municipal level and area sources on 250 m × 250 m resolution. The point source emissions indicate the emissions from major energy production and industrial plants, calculated based on several year average emissions. The point sources were aggregated to municipal level, including only sources within the Kokemäenjoki basin for each municipality. The FRES model uses several proxies to estimate the spatial distribution of the area source emissions (Paunu et al., 2013;Karvosenoja et al., 2018). The area sources are aggregated to 5 sectors: traffic exhaust (CO 2 ), machinery and off-road (CO 2 ), small scale wood combustion (CH 4 , N 2 O), other small-scale combustion (CO 2 ) and agriculture (CH 4 , N 2 O). The emissions from agriculture originate from fertilization (N 2 O) and cattle (CH 4 as enteric emissions). The emissions of methane (CH 4 ) and nitrous oxide (N 2 O) were converted to CO 2 eq to account for the climate warming impact of each substance. Results of two scenarios for future conditions, based on assumptions on energy use and mitigation measures of the national climate and energy strategy were used and derived for the Kokemäenjoki area: WEM (with existing measures) and WAM (with additional measures). Assumptions of these scenarios are described in detail in Koljonen et al. (2017).

Assessing biodiversity values
We developed three types of spatially explicit data on forest BD values in the study area: (i) distribution of forest stand types likely to include valuable BD elements, (ii) observation frequency of National Red List forest species, and (iii) modelled nesting suitability for four forestdwelling bird species indicative of forest stands with conservation value. In addition, we included (iv) spatial data on the dead wood potential of forest stands, described in Mikkonen et al. (2018) and Mikkonen et al. (2020). All the source data for these different BD data layers were converted or resampled using ArcGIS 10.5.1 tools (Esri) to a uniform 96 × 96 m grid covering the study area.
The distribution of four types of valuable forest stands was determined based on the spatial distribution of ecologically important forest features. The data for this mapping were extracted from three national sources of forest survey data, The Finnish Forest Centre, Metsähallitus Parks & Wildlife and the multi-source forest inventory conducted by Natural Resources Institute Finland. These three data were combined to cover the whole study area at 16 × 16 m resolution, providing detailed data on the mean forest stand age and key forest structure elements such as stand volume and dominant tree species (for details see Supplementary information S.1.2.1.). The first three forest types included mature and old-growth forests, delimited using the following age and main tree species criteria: 1) pine-dominated forests with stand age > 120 years, 2) spruce-dominated forests with stand age > 100 years and 3) deciduous-tree-dominated forests with stand age > 60 years. Due to the slow growth rate of boreal forests and intensive forest management in Finland, where stands are clear-cut at the age of 60-100 years, old-growth forests and their species have become nationally rare and/or endangered. The three mature and old-growth forest layers were derived in two steps. First, a 96-m grid cell was assigned into the respective dominant tree category if at least half (18/36) of the 16-m pixels within it met one of the above criteria (e.g., if ≥18 pixels had spruce as dominant tree species and average age > 100 years, it was assigned as spruce-dominated mature and old-growth forest cell). Second, the final value of the 96-m cell was calculated as the summed age of those 16-m pixels meeting the age criteria (e.g., for a sprucedominated cell, we summed the age of all pixels >100 years). This approach allowed us to identify those 96-m cells that had large proportion of old forest stands, while also differentiating between the ecological value of mature and the even rarer old-growth forests.
As the fourth type of forest stands, we delineated the locations of the 96-m grid cells which were spruce-dominated with a mean stand age > 60 years, but which also included deciduous trees higher than 20 m. Such forest stands represent sites where large deciduous trees of conservation importance, such as European aspen , occur as mixed stands within coniferous forest. Such forests are potentially valuable to many red-listed species, such as the Siberian flying squirrel (Pteromys volans) (Santangeli et al., 2013).
To assess the observation frequency of forest-dwelling threatened species in the study area, we extracted records of Critically endangered (CR), Endangered (EN), Vulnerable (VU), Near threatened (NT) and Data deficient (DD) species from the national HERTTA database for endangered species maintained by Finnish Environment Institute (see S.1.2.2.). We excluded records of extinct species and those made before the year 1990. In addition, we only included records with spatial accuracy of 100 m or finer. Each occurrence was assigned to a 96-m grid cell, and then a weighted sum of species observations, using the following weights: DD 1, NT 2, VU 5, EN 10, CR 20, was calculated for each cell.
For each species, we selected a set of explanatory variables that described (i) forest stand structure and quality, and (ii) type of land use at the focal cell, together with (iii) landscape level forest structurequality and land use around the nest (500-m and 1-km buffer for the Eurasian three-toed woodpecker and the three hawk species, respectively) (S.1.2.3.). We accounted for observation bias by building a density layer of all bird ringing points and used this to weight sampling of background points. We ran Maxent with auto-features and 10,000 background points, using a one-time data-splitting model validation and calculated the area under the ROC curve (AUC) in test data to measure model performance. We used the logistic output of final (full) model as predictions of nesting suitability for the four species.
There are well known dependencies between threatened forest BD and the amount, type and composition of dead wood parcels in boreal forests (Siitonen, 2001;Junninen and Komonen, 2011;Lassauce et al., 2011). Thus, dead wood volume provides a useful indicator of both forest naturalness and the presence of rare and threatened species. Dead wood potential (DWP) was estimated for each forest stand using a mechanistic forest growth model (MOTTI v 3.3, Salminen et al., 2005;Hynynen et al., 2014;Hynynen et al., 2015) and spatial data on forest site type and tree species (diameter at breast height (DBH) and volume), resulting in 20 combinatory data layers across 5 forest site types (herb-rich, herb-rich like, mesic, semi-xeric, xeric) and 4 tree species (Scots pine, Norway spruce, birches, other broad-leaved trees). These describe the potential DBH corrected volume of dead wood on site, summed across tree species and aggregated (summed) to each 96-m cell (Fig. 4C). The data have been used together with information on forest condition to identify important forests for BD conservation and have been shown to correspond well with ecological forest indicators (Lehtomäki et al., 2015; more details see Mikkonen et al., 2018 andMikkonen et al., 2020).

Landscape prioritization for integrated biodiversity and C evaluation
We used the spatial prioritization tool Zonation (v.4.0) (Moilanen et al., 2005;Moilanen et al., 2014) to identify jointly important areas for BD protection and C processes. The premise of such analysis is that tree harvesting in a currently forested location has negative impacts for both BD and C processes through the destruction of habitat and resources for forest-dwelling species, reduction of local C storage, and by shifting C net ecosystem exchange (NEE) from a sink to a source.
As input data, we used maps of C pools and fluxes (Section 2.3.1) and the different BD values (Section 2.3.2, Fig. 4). In the prioritization, we focused only on forests on mineral grounds (Supplementary S.1.4.1), as impacts of forestry actions on the net C fluxes of peatland forests are more complex and currently less well understood (despite these uncertainties, peatland forests were considered in the PREBAS modelling of the Cprocesses for the whole study area, Section 2.3.1). Here, to introduce the prioritization approach, we focussed in our analyses on data on C stored in forest trees (including roots) and the net ecosystem production (NEP), which includes soil respiration (role of stored soil C will be examined in future work). NEP describes the total amount of organic C available for storage (gCm −2 a −1 ) after accounting for gross primary production and ecosystem respiration in the PREBAS model. We only included positive values for NEP (setting negative values to zero) as we are only interested in including C sinks, not sources, to the priority areas. In total, we included 31 data layers: two for carbon (storage and NEP), and 29 for biodiversity values. As our aim is to demonstrate a joint prioritization for C and BD values, we weighted all input data equally.
Zonation produces a hierarchal priority ranking of each spatial unit (here 96 × 96 m grid cells) through an iterative process in which, at each iteration, the value of retaining (not harvesting) a unit is assessed based on input data, and the least important remaining spatial unit is removed from the analysis. The removal order defines the priority rank of each unit, the most important units being removed last. The amount of values remaining for each input feature is tracked throughout the prioritization, so that all features are captured in the solution in a balanced manner (Moilanen et al., 2011;Moilanen et al., 2014). Consequently, a set of top ranked priority areas together typically capture high value areas for all input features (Supplementary S.1.4.1).
In the prioritization, we accounted for past forestry actions (logging and management, Supplementary S.1.4.2) at each pixel by penalizing locations where forest management and drainage of wet forests have changed tree volume and the naturalness of forest. The penalty values (a multiplier with values 0-1) were constructed separately for C storage, NEP and two groups of BD values, i) dead wood potential on herb-rich soils; ii) all other BD values and were based on the intensity and frequency of management and time since management actions (Supplementary S.1.4.2). The ecological connectivity of sites was accounted for in two ways: For all features, we used kernel connectivity transformation at a 400 m mean decay distance to highlight areas where many high value sites occur near each other. We also accounted for connectivity between similar forest habitats, based on the spatial patterns of dead wood potential in different forest types (Section 2.3.2) (Supplementary S.1.4.1). The prioritization focused on unprotected areas, but we also included information from protected areas, so to identify those priority forests on mineral grounds that best complement the current reserve network.

Derivation of EO-based indicators for spatial assessment and monitoring
A large variety of EO data were acquired with different sensors and platforms and in different scales to produce species level information of forest BD (Table 2). Comprehensive field data were collected from the Evo test site (83 km 2 ), located in the south-eastern part of the Kokemäenjoki basin (Fig. 1), to support the EO approaches. These data consist of 400 circular sample plots, 2256 RTK-GNSS positioned trees, and leaf and soil samples, used as ground truth when developing methods for tree species detection and recognition of their biochemical characteristics.
Particular attention was given to detecting the keystone species European aspen, together with three main tree species (Norway spruce, Scots pine, birch). We used three different data sources with different spatial scales: Unmanned aerial Vehicles (UAV), aeroplane-operated sensors, and satellite data. First, the tree-level detection with highest spatial resolution was conducted using high-resolution photogrammetric point clouds and combination of spectral features derived from two different multispectral UAV cameras using linear discriminant analysis (LDA) and support vector machine. The results were validated with the fieldmeasured reference data. Second, to extend the geographical area for aspen detection, we used high-resolution airborne hyperspectral imaging and airborne laser scanning (ALS) data that have earlier been successfully utilized for individual tree detection and species classification in different forest environments (Dalponte et al., 2014;Fassnacht et al., 2016;Tuominen et al., 2018). Here, we assessed the role of different wavebands (455-2500 nm), principal component analysis (PCA) and vegetation indices (VI) in tree species classification based on canopy spectral reflectance using two machine learning classifiers: support vector machine (SVM) and random forest (RF) . In addition, the usability of convolutional neural networks (3D-CNNs) in tree species classification was assessed (Mäyrä et al., 2021). Third, we tested area-based estimation of species occurrence using Sentinel-2 satellite imagery that further extends the spatial scale of the study.

Results
3.1. Regional C budget and estimation of C neutrality of the study area Regional C budget patterns for forested areas were simulated with PREBAS, under both present and future conditions (Fig. 3, Tables 3 and 4). These results were used also for the landscape prioritization for integrated BD and C evaluation (see Sections 2.3.3 and 3.3). FRES model results for anthropogenic emissions for both present and future (WEM and WAM mitigation scenarios) were also derived for the area (Tables 3 and 4).
Forested areas formed the largest sinks in the present-day GHGbudget of the study area (Table 3). Forests were also the largest source together with anthropogenic emissions from artificial surfaces, while other ecosystems (arable land, surface waters and undrained mires) were lesser but quantitatively important sources. Lakes and rivers form a source of about 1.1 TgCO 2 eq a −1 (Holmberg et al., 2021; Tg = 10 12 g), due to decomposition of terrestrially derived C in the waters and are thus affecting the regional C budget (Kortelainen et al., 2006). Emissions from agricultural areas were estimated at ca. 0.6 TgCO 2 eq a −1 (Holmberg et al., 2021). Forests dominate estimated present-day sinks of the area and form a sink even after accounting for harvesting removals and decomposition/respiration processes (net emissions −3.7 TgCO 2 eq a −1 ). However, the anthropogenic emissions from artificial surfaces, emissions from forests (harvesting + decomposition/ respiration) and other ecosystems greatly exceed the total sinks in the area at present, with estimated total net emissions of ca. 4.2 TgCO 2 eq a −1 (Table 3). Thus, the Kokemäenjoki area is far from C neutrality, a general policy target set to be reached by 2035. Note that we do not account for C stored in wood-based products because our analysis does not include life-cycle calculations.
Mitigating anthropogenic emission sources indicated a limited potential for significant reductions by year 2030, with net emission of 5.1 and 4.8 TgCO 2 eq a −1 for the WEM (with existing measures) and WAM (with additional measures) scenarios, respectively. Longer-term policy scenarios for the anthropogenic emissions were currently not available. However, the study areas' forests presently provide a large theoretical sink of −11.5 TgCO 2 eq a −1 (Table 3). Under the Low harvesting intensity scenario, net emissions for forests amount to −5.8 TgCO 2 eq a −1 by the mid-century (Table 4), compared to the Policy (base) scenario outcome of −1.4 TgCO 2 eq a −1 . The net sink of all scenarios is affected by changes in the forest age structure caused by already implemented cuttings in the area. Harvesting according to the MaxSust scenario would convert forests into an emission source already by 2026-2033 (net emissions 1.6 TgCO 2 eq a −1 , Table 4).

Assessing biodiversity values
The three types of old forest stands (spruce, pine and deciduous forests) occur sporadically in the study area and their cover is very low (see Supplementary information S.2.2.4). The most abundant type of the three, spruce-dominated forests with mean age > 100 years, covers approximately 225 km 2 , thus accounting only for 1.15% of the forested land in Kokemäenjoki basin (Fig. 4A). Moreover, in all the three forest types, their overall cover decreases drastically towards older age categories indicating a severe general loss of oldest forests from the study area (S.2.2.4). Examination of the forest stand age data inside vs. outside PAs further highlights the sparsity of oldest forests outside PA network, as more than 70% of the spruce forests >160 years old, and of the pine forests >180 years old, respectively, are situated within PAs (S.2.2.4). Over >60 years old spruce forests with tall (>20 m) deciduous trees are more common and cover ca. 13% of the forested land. The total number of threatened forest-dwelling species included in the analysis was 235 (4 DD; 120 NT; 66 VU; 34 EN; 11 CR). The mean number of occurrences across all species was 49.2, ranging from one occurrence (67 species) to 7494 (Siberian flying squirrel, VU). The number of occurrences among critically endangered species (CR) varied from one (Scapania apiculate and Fuscocephaloziopsis catenulate) to 20 (Nephroma laevigatum), and in endangered (EN) from one (12 species) to 1915 occurrences (eastern pasque flower Pulsatilla patens), respectively. The summed Red List species index varied among the 96-m cells from zero with no records of the Red List species to 65 (Fig. 4 B). The spatial patterns in the summed Red List index are notably driven by the occurrences of two well-surveyed species, Siberian flying squirrel and eastern pasque flower.
The final Maxent models predicting the nesting site suitability for our four focal bird species (see S.1.2.3.) included 6-9 predictor variables depending on species (S.2.2.5). The most important variables (with >10% contribution) predominantly described forest stand characteristics at 96-m resolution, including age, height and volume of trees, with predicted suitability showing a positive relationship with these variables (S.2.2.2 and S.2.2.3). Of the landscape level variables, cover of urban areas showed a notable negative relationship with nesting suitability for common buzzard and northern goshawk, whereas stand suitability for the Eurasian three-toed woodpecker increased with increasing mean tree volume and cover of forested peatlands in the surrounding 500 m buffer. Spruce was typically the dominant tree species in stands suitable for the three-toed woodpecker.
The predictive discrimination performance of the models, as measured by Maxent's test AUC statistics, was clearly better than random (i.e. 0.5; see Phillips et al., 2009) for all four species, varying from 0.78 (common buzzard) to 0.91 (Eurasian three-toed woodpecker) (Supplementary Table BD2). The maps of predicted nesting suitability showed spatial differences between species. For the common buzzard and the northern goshawk, the most optimal sites were relatively widely and evenly spread (S.2.2.1) whereas for the European honey buzzard nesting suitability decreased towards coastal areas (Fig. 4 C). Suitable nesting sites were spatially most limited for the Eurasian three-toed woodpecker (S.2.2.1). For all four species, only 3.5-6% of their total predicted values were estimated to reside in protected areas (S.2.2.5). This suggests that the way how the optimal sites for the four bird species are managed outside PAs can have an important impact on the regional persistence of these species.

Landscape prioritization for integrated biodiversity and C evaluation
Integrating information about the regional C budget (Section 3.1) and BD values (Section 3.2), we explored how BD values and C storages and sinks could be retained within the Kokemäenjoki area under four optional conservation scenarios: i) only currently protected forests on mineral ground are retained, or; in addition to protected areas, harvesting would not be allowed in the top 10% ranked unprotected locations (~190,000 ha) identified using ii) biodiversity values only; iii) C values only, or; iv) both BD and C values (called 'balanced' here onwards) (Fig. 5 A, C-D). In total, the currently protected forests on mineral soils store approximately 2.6 Tg of C and assimilate another 0.06 Tg each year (0.22 TgCO 2 eq a −1 ), representing~2.4% of the total C pools and fluxes of mineral soil forests in the region. On average, the protected areas capture~10% of BD features' values on mineral soil forests (see Section 3.2 for more on protection).
With optimal allocation of no-harvest areas, the amount of both BD and C values retained in the study region could be significantly improved (Fig. 5B). However, using either BD or C values to identify these priority forests resulted in spatially divergent patterns (Fig. 5C- Table 3 Estimated current GHG emissions and sinks (TgCO 2 eq a −1 ) by land use in the Kokemäenjoki area (modified from Holmberg et al., 2021). Positive numbers indicate emissions and negative numbers sinks. Emissions from forests are due to cutting removals and decomposition and respiration processes. The land use class 'Other ecosystems' refers to arable land, surface waters and undrained mires.

Land use
Area (km 2 ) Area (%) Emissions TgCO 2 eq a −1 Emis. wrt total (%) Sinks TgCO 2 eq a −1 Sinks wrt total (%) Net emissions TgCO 2 eq a −1 D), indicating that best locations for conservation and climate change mitigation do not necessarily overlap. Despite spatial differences, using either C or BD values to prioritize areas made relatively little difference to how much more C pools and fluxes could be increased (Fig. 5B). On the other hand, BD values are poorly captured by areas prioritised only for C, in comparison to the BD only scenario. The balanced solution (Fig. 5A) closely mimics the BD only solution: this is because options for increasing C pools and fluxes are numerous in the region, whereas many of the BD values have restricted distributions (Fig. 5) and hence constrain the solution more (Kujala et al., 2018b). The priority areas identified in the balanced solution nevertheless contain similar C pool as the C only solution (12.1 vs 13.2 TgC, respectively), whereas C fluxes are smaller in the balanced solution (0.20 vs. 0.43 TgC a −1 in the C only solution) but still considerably higher than fluxes in just protected areas (0.06 TgC a −1 ).

EO-based tree species detection for indicator development and monitoring
For providing new datasets of biodiversity indicators, EO techniques can be utilized. UAV-based remote sensing proved to be very efficient in providing ultra-high spatial and temporal resolution RGB, multispectral or hyperspectral imagery for a detailed forest properties assessment in boreal forests at individual tree and stand-level. The best classification results were obtained with linear discriminant analysis which gave 82% accuracy for the recognition of the keystone species European aspen. Similarly, European aspen could be detected with high accuracy using high-resolution airborne hyperspectral and airborne laser scanning (ALS) data. To extend the geographical area from the UAV based, tree level aspen discrimination, we used high-resolution airborne hyperspectral imaging and ALS. The best classification accuracy of 92% for aspen from the other dominant species was reached using the support vector machine (SVM) classifier with mean spectral reflectance values together with selected vegetation indices (cf. Viinikka et al., 2020). Overall, the best method for tree species classification considering all species was 3D convolutional neural network (3D-CNN) utilizing 250 spectral bands, with overall accuracy of 86% (cf. Mäyrä et al., 2021).
Very-High-resolution EO data enabled the identification of the species at tree-level with reliable and accurate results . However, data obtained from satellites with lower spatial resolution is mainly suitable for area-based estimation of species occurrence. Therefore, upscaling of the tree species detection using satellite data to large geographical areas seems challenging. The spatial resolution affects the identification of individual trees, which is demonstrated in Fig. 6. Whereas canopies can be delineated in Lidar data and individual trees can be recognized from UAV and airborne imagery, in Sentinel-2 data, each pixel consists of several trees which makes it difficult to detect the spectral traits of a single scarcely occurring tree species. Our preliminary results on tree species detection with Sentinel-2 data showed that correlations between spectral bands and species proportions are rather weak, especially in case of species with scarce occurrence, such as aspen and larch (Larix spp.).
At current state, we have been able to produce estimates of aspen occurrence and abundance on an 83 km 2 area from fusion of airborne hyperspectral and Lidar data ( Fig. 7; Viinikka et al., 2020). This indicates the potential of the derived methodology for spatial BD mapping and identification of new areas for protection.

Discussion
Fig. 2 presents the integrated modelling and evaluation framework of our study and shows how the different sectors are linked. Our discussion is structured based on this conceptual framework and we describe how it can be used for finding optimized solutions for co-management of C and BD in the landscape.

Assessing C pools and fluxes and C neutrality
Estimation of C budgets is a key necessity in assessing whether and by which means warming can be constrained to the limits of global average temperature increase, as set out in the United Nations Paris Agreement (1.5°C and well below 2°C, Rogelj et al., 2019). Reaching this goal requires both rapid transformation of national energy, industrial and land-use sectors and large-scale deployment of negative emissions technologies (NETs) and other land-based measures (Smith et al., 2016;Roe et al., 2019). These NETs and land-based mitigation include techniques such as increased use of bioenergy combined with C capture, reduced forest deforestation and degradation, less intensive forest management, and afforestation. These measures are also central for reaching 'net zero emissions' in climate mitigation target setting. Forests thus play a major role in the global and regional climate change mitigation efforts.
In Finland, forests cover ca. 70% of the land area, and the forest industry is an important part of the national economy. Total Finnish GHGs emissions were 52.8 TgCO 2 eq in 2019 and the LULUCF sector was a sink, varying from approximately 19 to 50% of the total annual emissions from other sectors during 1990 to 2017 (Statistics Finland, 2020). Forests (trees and soil) thus absorb a significant proportion of Finland's CO 2 emissions. C budget calculations are therefore of importance for decisions on land-use policies and management at different spatial scales.
Our results show the large importance of forests in the regional C budget also of our study area (estimated present net emissions −3.7 TgCO 2 eq a −1 , Table 3; negative number indicating a sink) and, subsequently, the potential for using management measures of forests, other land-use sectors and anthropogenic emission reductions in mitigation. The possibilities for reaching policy goals of C neutrality are clearly demonstrated (Table 4), indicating that C neutrality could be achieved but would require significant reductions in the current forest cutting levels. As shown in Table 4, the estimated mitigated anthropogenic emissions in the area by year 2030 amount to 5.1 (WEM scenario, planned existing measures) and 4.8 TgCO 2 eq a −1 (WAM scenario with stringent additional measures), respectively. The net emissions of forests in the area according the Low intensity harvesting scenario by around 2030 are estimated to −6.6 TgCO 2 eq a −1 (sink) but would be significantly lower for the Policy (base) scenario for the same period (−1.2 TgCO 2 eq a −1 , Table 4). Cuttings close to a Low scenario could thus compensate for mitigated anthropogenic emissions by around 2030. Forest harvesting according to the high-intensity MaxSust scenario would convert forests into an emission source (Table 4). Our model system can therefore be used to assess the impacts of different scenario combinations for different time frames, allowing long-term planning and adaptation. Municipalities and other regional actors can use this information to make decisions that affect local GHG emissions and sinks, in particular via land use planning and management (Vanhala et al., 2016). Importantly, after clear-cutting, a boreal forest stand becomes a source for C and other GHG emissions for several decades (Grant et al., 2010;Korkiakoski et al., 2019), indicating the long time-frame involved. Our results also clearly demonstrate the potential for integrated C and BD management (see Section 3.3 and discussion below). Compiling comprehensive spatial data for BD has well-known challenges. Although forest attributes provide useful and relatively accessible proxies for BD at larger scales, conserving BD values often requires data on the exact distribution of individual species. Observation records of species are becoming rapidly available and are increasingly used in various BD assessment (GBIF Secretariat, 2019). However, the pointnature of most species data limits their use in landscape level analyses, such as integrated land-use optimizations, as they include large unsurveyed areas where species presence is unknown. One option is to combine records of species within relevant BD groups, as done here for the forest-dwelling threatened species, which allows exploration of BD patterns in broad terms (Section 3.2). Such approach is particularly suitable for including very rare, elusive or taxonomically cryptic species for which data quality is typically low (Edenius and Mikusiński, 2006;Becker and Encarnação, 2012). Combining data across species may alleviate data gaps and biases, thereby enabling identification of valuable forest stands for poorly known species (Mikusinski et al., 2001;Wiens et al., 2008;Burgas et al., 2014).
For species with enough observations (>30), Species Distribution Models (SDMs) offer a powerful tool for filling data gaps in unsurveyed parts of the landscape. There are numerous established SDM methods and their use in a range of conservation decisions is proliferating (Guisan et al., 2013). Here we modelled the nesting suitability of forest stands for four bird species. Earlier research has suggested that these species are dependent on the presence of old, mature forests with high stand volume, and their optimal nesting sites can indicate highquality mature forest in terms of BD value (Pakkala et al., 2014;Burgas et al., 2014;Björklund et al., 2015). Our results, based on nest data and fine-grained forest data accumulated across the whole Kokemäenjoki basin (Section 3.2), are in line with these earlier findings with respect to the study species ecology. Important to acknowledge is that most BD observation data are spatially biased towards easily-accessible and frequently visited locations and that it is central to account for these biases when building SDMs (Phillips et al., 2009).
In our Kokemäenjoki study, we used complementary BD measures known to be relevant for forest BD conservation, and which reflect more than just one element of forests (i.e., forest maturity, mixed tree composition, amount of dead wood and forest characteristics at local and landscape scale). From forest conservation planning perspective, more data on ecological forest attributes is needed, as e.g. many threatened forest species are epiphytes, favoring old stands or stands with large aspen trees. In managed forests these mature forest stands are systematically harvested and clear-cut (Björklund et al., 2020;Virkkala et al., 2020), causing declines in mature and old-growth forest associated species, such as our focal forest birds, three of which are already regarded as red-listed (Hyvärinen et al., 2019). In the study region, the remaining old mature forests appear to be retained predominantly in protected areas.

Improving spatial BD data compilation using EO techniques
The diverse structure of nature makes it difficult to identify and map surrogates that accurately represent the many levels of BD and ES values (Rodrigues and Brooks 2007). Combining EO and field-measured training data, some surrogates, such as stand mean age, are particularly challenging to define, whereas others (i.e., stand volume) can be mapped reliably from appropriate data sets (Vastaranta et al., 2016). Largescale data on BD relevant forest attributes can be, to some degree, sourced from national forest inventories (Section 2.3.1), which combine remote sensing data with intensive field surveys (Tomppo et al., 2008;Pukkala, 2020a). However, the maximum level of detail depends on the amount of field data and the EO material used in the process. Scarce tree species (e.g., European aspen) and stochastic phenomena (e.g., deadwood) that are not directly linked to visible canopy cannot be mapped reliably at large scale. In Finland, current approach providing forest resource data for operational forestry struggles in discriminating European aspen and other deciduous trees as the field data are typically targeted for mapping the dominant and economically significant tree species (Scots pine, Norway spruce and birch) and stands.
Our results showed that EO data with very high spatial and spectral resolution provide spatially accurate information of tree species that can help to target BD conservation efforts in boreal forests (Section 3.4). Previous studies have made similar conclusions (Säynäjoki et al., 2008;Dalponte et al., 2014;Puliti et al., 2015;Fassnacht et al., 2016;Nevalainen et al., 2017;Alonzo et al., 2018;Tuominen et al., 2018). According to our results in the Evo test site, UAV data, airborne hyperspectral reflectance and ALS data enable the tree-level detection of ecologically important tree species, such as European aspen Mäyrä et al., 2021). However, spatial coverage of UAV data is often limited, and airborne data are overly expensive to acquire for large areas. Thus, there are still obvious challenges in gathering of ecologically useful data over wide geographical regions. In the future, linking highly detailed UAV data with high resolution satellite imagery is a viable option for tackling these challenges (Kattenborn et al., 2019).
The spatial and spectral resolution of EO data greatly affects what kind of attributes may be mapped. Low-density (≈ 1 pulse/m 2 ) ALS data accompanied with aerial images and extensive field data can produce robust area-based estimates of features such as stand mean height and volume (Naesset, 2002;Hyyppä et al., 2008), whereas higher density ALS data supplemented with field measurements are often required for detecting single trees (Hyyppä et al., 2001), or to describe the vertical structure of forests (Lindberg et al., 2012;Vastaranta et al., 2016). Multispectral data from high spatial resolution commercial satellite sensors have been available for two decades, but their high cost and limited spectral resolution limit their operational use. Currently, Sentinel-2 is the most suitable freely available satellite data in terms of spatial (10 m) and spectral (13 bands) resolution. Our preliminary results indicate that the high classification accuracies obtained at tree-level cannot be easily repeated using Sentinel-2 data for area-based estimation of species occurrences (see also Stratoulias et al., 2015;Immitzer et al., 2016;Radoux et al., 2016), particularly when trying to detect scarcely occurring tree species among the dominant species. Moreover, the broader spectral bandwidths of Sentinel-2 cannot capture the subtle differences in spectral reflectance which are important in accurately discriminating between tree species. To overcome the issue, efforts of multi-temporal monitoring by Sentinel-2 data have given promising improvement (Bolyn et al., 2018;Persson et al., 2018). There are also other satellites either in the orbit or satellite missions in preparation with very high spatial (e.g. WorldView-3) or spectral (e.g. GaoFen5, PRISMA) resolution which could enable more accurate large-scale detection of rare tree species. The data from these missions, however, is not yet openly available, or requires heavy processing to be usable. Our results clearly demonstrate the potential for fusion of different EO based data for spatial BD mapping and identification of new areas for protection (Fig. 7).

Enhancing joint spatial planning and prioritization of BD and C benefits
There is a great need for integrative research and conservation planning in forest environments to identify feasible options, environmental risks and sustainable policies. A key objective of the national strategy and action plan for BD conservation is to halt BD loss in Finland. The strategy includes several measures such as consideration of climate change impacts, halting the decline in forest species and habitats, enhancing the network of protected areas and accounting for the maintenance of BD and ecosystem services (ES) in land use planning. The recent EU BD strategy has similar objectives. According to Hyvärinen et al. (2019), 2667 species in Finland (11.9% of the assessed species) are classified as threatened (CR, EN, VU), and only some 2% of the forested area in southern Finland is protected.
Considering BD conservation is essential also for GHG management strategies. A key aspect in reaching BD and C policy targets is landscape and regional level land-use planning (Thomas et al., 2013;Moilanen et al., 2014;Kujala et al., 2018a;Buotte et al., 2020;Reside et al., 2020). With optimized planning, intensive commercial forestry can be directed to locations less valuable for BD and C storage/sequestration, and BD-friendly management and local conservation actions applied in areas with notable BD and ES values. In addition, different sets of approaches, including varying intensity of forestry actions along with protection and nature management, can be applied, accommodating consideration of other land-use aspects (e.g. limited conservation budgets, private land owner interest). In the Kokemäenjoki basin, for example, the reduced forest cuttings needed to reach the C neutrality target (Table 4) could be allocated to areas with high conservation and C values (integrated solution, Fig. 5A), which would significantly improve the amount of both BD and C values retained in the study region. To achieve this outcome by the choices of private landowners requires reshaping current economic incentives, such as METSO programme (Forest Biodiversity Programme for Southern Finland), to account for both BD and C values (see below).
As our results on spatial BD mapping and EO based approaches indicate, gathering robust data for detailed landscape-level spatial optimizations is not a trivial task. Our study also demonstrates how C measurements from intensive research sites can be extrapolated to larger regions using modelling tools and national forest inventory data. Similar approach has been used, for example, in the design of research infrastructures for ecosystem research where intensive sites with heavy instrumentation are integrated into larger sites for socioecological research (so-called LTSER sites, Mirtl et al., 2018). Nevertheless, it would be valuable to improve the existing extensive measuring schemes with simple but informative additional measures such as soil C and additional BD indicators.

Using economic incentives for enhanced BD conservation and C sequestration
Implementing an integrated approach for the enhancement of BD conservation and C sequestration requires support from economic incentives, and active voluntary participation of forest owners is key for the success when forest land is privately owned. In Finland forest owners accept economic incentives for voluntary BD conservation (e.g. Horne, 2006;Juutinen and Ollikainen, 2010). Inclusion of climate targets in the incentives provides, however, a special challenge requiring co-creating through forest owners' participation to ensure social acceptance. An ideal participatory payment mechanism provides an extension of the existing METSO compensation scheme for BD conservation to C storage/sequestration (Kangas and Ollikainen, 2021;Kosenius, 2021). At the heart of the process, there is a need to develop economic mechanisms of joint ES provision (Lankoski et al., 2015;Ollikainen, 2016) and the associated trading rules, as well as means to assess the coherence of various environmental policies and their acceptance among the stakeholders. It has previously been shown that subsidizing forest C sequestration by 50 € t −1 (metric ton) would increase the C sequestration of Finnish forestry by 50%, and with higher subsidies even more (Pukkala, 2020b).

Potential for fully integrated multi-objective spatial prioritization
There are also apparent challenges in developing and successfully applying multi-objective plans for regional land use and conservation. Prioritization methods can be tailored to find optimal "win-win" solutions, accounting for several factors such as expected climate change, land-use change, land cost, area size, ecological connectivity, and other ESs, and willingness of land owners to participate in conservation schemes. However, the diversity of elements and how they are valued and traded-off against each other make it a highly complex issue (Moilanen et al., 2014;Kujala et al., 2018aKujala et al., , 2018b. Our study showed that it is possible to optimise for both C and BD values in a landscape, but also that the best locations for conservation and climate change mitigation do not necessarily overlap (Section 3.3, Fig. 5). Similar results have been obtained in other studies (Buotte et al., 2020;Repo et al., 2020), indicating the need for diverse strategies. These should accommodate a range of actions that can promote multiple goals, including continuous-cover silviculture, lengthening of harvest cycles and conservation of managed landscapes Buotte et al., 2020;Eyvindson et al., 2020;Repo et al., 2020). Indeed, the next step in developing the multi-objective spatial prioritization approach used here is to include different types of BD-friendly forest management measures as well as land owner-based socio-economic elements into one analysis. Despite their complexity, such analyses have been demonstrated to be very useful for uses such as zoning, conservation planning, or impact avoidance in development (Whitehead et al., 2017;Jalkanen et al., 2020).

Conclusions
Climate change and loss of biodiversity are massive global challenges and deeply interconnected problems. Efforts to develop common solutions at different scales and in different ecosystems are therefore required. As demonstrated in this paper, developing a detailed integrated multidisciplinary modelling and evaluation framework, where different carbon processes and biodiversity aspects can be analysed simultaneously and where optimized solutions in the landscape are sought for, can assist in these efforts. We link these analyses to formal environmental policy goals of the society. The next step will be to expand the tested methodologies to the national scale, thereby further increasing the policy-relevance of our approach. Still, it is likely that much of the real progress aiming at finding integrated solutions to these problems will be achieved at the local and regional scale where the actual landuse decisions are made, and anthropogenic emission mitigation efforts are planned and implemented. For the data needs of these actions, the detailed spatial scale of our approach is well suited.
Our results show that it would be possible to achieve the stated policy goal of carbon neutrality in our study area by mitigating anthropogenic greenhouse gas emissions and substantially decreasing the level of forest harvesting, and that one can further utilise this information to optimise for retaining both biodiversity and carbon values in the landscape. Such a policy can be implemented by a combination of anthropogenic emission mitigation actions, municipal and regional landuse planning, and voluntary actions by land-owners, incentivised by economic instruments for joint biodiversity protection and carbon storage/sequestration. The impacts of these actions in the landscape can be studied and monitored by utilizing advanced Earth Observation-based techniques and indicators.
However, our study also shows the complexity of the issues and the large data needs and heavy computing resources required to carry out the detailed analyses. Moreover, we recognise the large uncertainties involved in the ecosystem processes studied, in the model simplifications made, and in the data extrapolation efforts carried out. Continued intensive research, data collection and mapping efforts are therefore needed. The long time-scales involved also necessitate frequent updating of the different impact and management scenarios to minimize uncertainties and keep the results relevant for the end-users.

Role of authors
MF planned and organised the study and wrote/edited the first version of the paper. HK, NM, SKu, RM, AM, RV, MPek and RKK collected biodiversity data and site information, carried out Maxent and Zonation based modelling and mapping work and participated in writing of the paper. FM, MH, AA, JB, AM, KR, TR and PVa collected data on carbon processes, carried out modelling work based on the PREBAS model and landscape level analyses on carbon processes, and participated in writing of the paper. NL and IA carried out GIS-based data aggregation and analyses. V-VP and NK carried out modelling work using the FRES model on anthropogenic greenhouse gas emissions. TT, PH, JM, SKi, SK-S, AV, PVi, TK, AK, and LP carried out the analyses using Earth Observation techniques and participated in writing of the paper. A-KK and MO provided information on using economic incentives in land-use planning and provided comments on the paper. MPel and ST provided data on forest segmentation used in carbon modelling.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.