Identifying Nutrient Export Hotspots Using a Spatially Distributed Model in Boreal-Forested Catchments

Abstract

The implementation of the Water Framework Directive (WFD) aimed to reduce nutrient export from catchments to water courses. Forest operations cause diffuse loading, which challenges the efficient targeting of water protection measures. We formed 100 equally probable clear-cut scenarios, to investigate how the location of the clear-cuts influenced the total nitrogen (TN) and phosphorous (TP) export on different scales. The nutrient export was calculated by using a distributed nutrient export model (NutSpaFHy). The clear-cut-induced excess TN and TP exports varied by 4.2%–5.5% and 5.0%–6.5%, respectively, between the clear-cut scenarios. We analyzed how the sub-catchment characteristics regulated the background export. The results also suggested that there was no single sub-catchment feature, which explained the variation in the TN and TP exports. There were clear differences in the background export and in the clear-cut-induced export between the sub-catchments. We also found that only 5% of the forest area could contribute up to half of the total nutrient export. Based on our results, we presented a conceptual planning framework, which applied the model results to finding areas where the nutrient export was high. Application of this information could improve the overall effectiveness of the water protection measures used in forestry.

1. Introduction

The European Union’s Water Framework Directive (WFD) forbids deterioration of the chemical and ecological status of water bodies [1]. If the status of the water body is worse than ’Good’ in the WFD classification scheme, then active measures, both in the water bodies and in their catchments, are required, to reclaim the ’Good’ status [2]. The aim of water protection is to reduce the export of nutrients and suspended solids to a level where the recovery of the water body is possible. In Finland, boreal forests cover approximately ¾ of the total land area [3], and the majority of this area is privately owned and commercially managed; therefore, forestry plays a key role in implementing the WFD obligations.

Forestry operations, such as final fellings, cause diffuse nutrient load [4,5], meaning that the elevated nutrient export emerges in large areas simultaneously. A typical clear-cut size in southern Finland is less than 2 hectares [6]. Because the rotation period from regeneration to the final felling is typically tens of decades, less than 2% of the forestry area becomes clear-cut annually [3].

The clear-cut and the following site preparation alter the water and nutrient fluxes [4,7], and induce increased nutrient leaching. The nutrient cycle in boreal forests is tight [8], and excess nutrients are efficiently taken up by trees, ground vegetation and soil microbes [9,10]; nevertheless, elevated nutrient concentrations in runoff water are observed regularly for several years after the clear-cut [4,11,12,13,14,15,16,17,18,19,20], because the removal of the tree biomass decreases interception and transpiration, raises the water table and increases the runoff [21,22]. At the same time, litter decomposition [23] and nitrification rates are accelerated [24,25], due to the increased soil temperature [26]. The nutrient supply exceeds the nutrient demand, as the uptake by trees collapses [27,28] and the remaining logging residues release nutrients [29]. Site preparation can further decrease the ground vegetation biomass, and add decomposing litter to the soil, thus increasing the release of nutrients [30].

Recent studies suggest that total nitrogen (TN) and phosphorus (TP) exports are not decreasing, either on the national scale [31,32] or on the headwater catchment scale, in Finland [5,33,34], which may jeopardize the attainment of the WFD goals [31]. The main drivers behind this increasing trend in nutrient export are likely the gradual changes in the climate, and atmospheric deposition [5]. The trend is pronounced in peatland-dominated catchments [32], which are very common in Finland, as the peatlands cover a fifth of the total forest area [3]. The peatlands contain a vast storage of organic matter, which is susceptible to enhanced biological decomposition by increased soil temperatures [35]. The effect is further amplified if the water table lowers and exposes peat to aerobic decomposition [36]. With decomposition, nutrients are released [37] to the soil solution, and can be exported further to the water bodies, thus increasing nutrient exports to water courses.

The above, together with increasing harvesting, calls for efficient methods of planning and implementing water protection measures [2,38]. Diverting water protection investments equally over the whole area or applying similar water protection methods everywhere decreases the overall cost efficiency of the water protection [39]. Thus, there is an urgent need for the development and implementation of more efficient water protection that will allow specific targeting of the water protection measures to places where the nutrient export risk is highest [40], considering the catchment characteristics regulating the export [41]. In this study, we applied a distributed nutrient model, to identify the high nutrient export risk areas, and we discuss whether the identification can be used to direct water protection measures efficiently.

The magnitude and duration of the nutrient export caused by clear-cutting depend on several location-specific factors, such as catchment topography, soil properties, site fertility, atmospheric deposition, vegetation recovery, the timing of management practices, weather conditions and distance to water courses [4,7,19,30,41,42,43,44,45]. Catchment characteristics contribute to the timing of elevated nitrogen (N) concentrations in stream water after the clear-cutting [46]. Water N concentrations can increase quickly in areas with high atmospheric deposition and/or fertile soils. By contrast, the response can be delayed and low in low-fertility and low-atmospheric deposition areas, such as Finland and central and northern Sweden [41,43,47], because the growth of boreal upland forests is nitrogen-limited [8]. Any surplus N from the atmospheric deposition, or any other source, will efficiently be retained in the stand, ground vegetation, soil and microbes; consequently, the leaching losses will remain very low [48], especially in low-fertility sites [49]. In sparsely populated areas in central and northern Finland and Sweden, the atmospheric deposition is low [50].

Statistically significant effects of clear-cutting on water quality can be difficult to separate by using a paired catchment method [4,51], because the treatment effect can become masked behind the large natural variation in the nutrient export [51]: this highlights the need to account for the processes affecting the nutrient retention along the flow paths. Based on paired catchment studies, Palviainen et al. [4] found that clear-cutting significantly increased N and phosphorus (P) exports only when the clear-cut area exceeded 30% of the total catchment area.

Blackburn et al. [45] concluded that N export to stream is driven by biogeochemical processes and hydrological routing in riparian zones, emphasizing the role of the near-stream areas as nutrient export hotspots or nutrient-retaining buffer zones. Undisturbed riparian buffer zones can considerably decrease the nutrient export [9] because, during the transport, nutrients can be retained in soil, vegetation and microbes before reaching the water body. However, if near-stream areas are disturbed, the nutrient retention may decrease and, consequently, the nutrient transport to the stream increases [45].

The specific load method [4,52] is currently used in Finland to estimate forestry-induced nutrient export to watercourses. The specific load method is based on coefficients that are determined separately for each management measure at mineral and peatland sites (kg ha⁻¹year⁻¹) using paired catchments experiments [4,19,51,53,54,55]. The specific load coefficients represent the average excess load caused by the management measure. The specific load method is easy to apply, but the experimental quantification of the coefficients takes years, is very expensive and the number of paired catchment studies for a management method remains small. Furthermore, the method is not responsive to changing weather conditions, location of treatments, or interconnected biogeochemical processes changing along the flow path of the water [9,40,41,56]. Therefore, it is necessary to develop and use methods that include the most important biogeochemical processes, catchment properties and the explicit locations of clear-cuts.

We used the NutSpaFHy-model [57] to estimate how the positioning of clear-cuts affected TN and TP exports in an area consisting of 33 sub-catchments in southeastern Finland. NutSpaFHy is a spatially distributed nutrient export model incorporating grid-(16 m × 16 m)-based nutrient balance calculation into a hydrological model SpaFHy [58], and including a conceptual solute transport routine to approximate N and P exports to the sub-catchment outlets [57]. The model was parameterized using up-to-date forest inventory data, which allowed for realistic positioning of potential clear-cut sites. The model results were explored, to identify nutrient export hotspots, and to consider catchment characteristics contributing to nutrient export.

By comparing the modeled nutrient exports between the sub-catchments in different clear-cut locations, we aimed to answer the following questions:

i

How much does the nutrient export vary according to the location of the clear-cuts, on catchment and sub-catchment scales, and by how much can the nutrient export be reduced by changing the locations of clear-cuts, while keeping the total clear-cut area constant?

ii

What is the importance of the background export in the total nutrient export, and what are the most influential sub-catchment characteristics regulating the background load?

iii

How can nutrient export hotspots be identified, mapped and used in the water protection planning process?

We hypothesized that the locations of clear-cuts, and catchment characteristics, have a remarkable impact on nutrient export. In a mature forest, nutrient uptake and release are in balance, and only a small fraction of these rather large nutrient fluxes can be transported to water courses. Clear-cuts change the balance, and more nutrients are at risk of transporting to the water course if the site-specific characteristics support the transport to water courses. We wanted to test this assumption, in order to understand the significance, to the responsible management of boreal forests, of identifying high nutrient export areas (i.e., hotspots).

2. Materials and Methods

2.1. Study Area

The Vehka–Kuonanjärvi study area is a 7328 ha catchment situated in the boreal zone in southeastern Finland, South Savo region (61°59${}^{\prime }$ N, 29°12${}^{\prime }$ E), which drains to Lake Puruvesi. The catchment consists of 33 sub-catchments (Figure 1) whose forest areas vary between 9 and 1036 ha, the average being 158 ha. Water covers 22% of the total study area. The study area included several larger water bodies and numerous smaller water elements, which are connected by stream networks.

Forestry land covers 90% of the total land area (5723 ha). Agricultural fields cover 6% of the total land area: these were not included in the nutrient export modeling. Roads, lots and other constructed areas (4% of the study area) were also omitted from the simulations.

The dominant tree species are Scots pine (Pinus sylvestris L.) and Norway spruce (Picea abies Karsten), mixed with deciduous species, downy and silver birch (Betula pubescens Ehrh. and P. pendula Roth) together with European aspen (Populus tremula L.). Due to a long tradition of even-aged forestry, the forests consist of compartments (stands), where the age and stand structure are rather homogeneous within the compartment, but may vary between the compartments.

Peatlands cover 26% of the total forest area. The majority of the peatland forests (85%) were drained for the first time between 1950 and 1980. As a result of drainage, originally sparsely wooded fens or open bogs changed to productive forests, with growth and biomass comparable to upland forests. Only 2% of the forest area is classified as non-productive land (volume growth < 1 m³ ha⁻¹ year⁻¹), typically located in non-drained or infertile peatlands [37]. The site fertility class (sfc, Finnish classification [59]) reflects wood production potential and soil nutrient content and, thus, the risk of nutrient export after clear-cutting. Sfc varies in the study area from highly productive (sfc 1) herb-rich mineral and similar eutrofic peatlands to non-productive barren mineral soils, alluvial lands or non-productive Spaghnum fuscum-dominated oligotrofic peatlands (sfc 6), as shown in Figure 2.

2.2. NutSpaFHy Model

First, a daily time step was used to calculate the hydrological fluxes and storages, which were further aggregated to support the calculation of nutrient balance and export by monthly time step [57]. The change in the local nutrient storage was calculated individually for each of the 16 m × 16 m grid cells of the forest land. The nutrient balance calculation included organic matter decomposition, atmospheric deposition and nutrient uptake by trees and ground vegetation. The excess nutrients from each grid cell were transported to the catchment outlet, using a conceptual model, which accounted for nutrient retention and delay as a function of the distance to the nearest water body (for a detailed process description, see [57]).

We used the parameterization introduced in [57], where N immobilization parameters were based on the site type and site fertility class, the mineral and peat soil distribution, the total tree volume and the coniferous tree volume in the sub-catchment. The P immobilization parameters were based on the averages calculated from the calibrated parameters in [57]. The effect of the clear-cut on the local nutrient storage was performed by updating tree age, volume, height, leaf area and ground vegetation biomass, as described in [57]. The simulation period extended from 2002 to 2012, and the first year was used to initialize the hydrology and nutrient storage calculation, and was therefore excluded from the analysis of the simulation results.

Catchments typically contain very different sites, in terms of nutrient export generation. In sites located on mineral soils, the N load in water percolating under the rooting zone can be very small (0.1…0.2 kg ha⁻¹ year⁻¹ [48]). On the other hand, N export to water courses can be 10 times higher in pristine peatlands [20]. Because a catchment is a mosaic of different sites, temporal and spatial variation in hydrological and biogeochemical processes introduce uncertainty to any estimation of forest management effects on nutrient load, as shown by Laurén et al. [51]. Therefore, studying treatment effects using a process-based model is a particularly suitable approach, because the harvesting scenario can be changed while keeping all other parameters and input variables unchanged. Consequently, the differences between the harvested and unharvested scenarios reflect the direct and indirect treatment effects, and there is no risk that the background export will mask the treatment effect in the modeling approach.

The NutSpaFHy model was tested against experimental nutrient export datasets measured from the catchment outlets in Laurén et al. [57]. In a situation where a model is evaluated against one measured dataset, and the model output can result from a combination of multiple processes and their spatial arrangements, there is always a risk of obtaining the right results for the wrong reasons [60]. To decrease this risk, Laurén et al. [57] relied on conceptualization and construction of the model structure, keeping the number of calibrated parameters at a minimum. In NutSpaFHy, high export emerges in locations where the nutrient release exceeds the nutrient uptake for a prolonged time, and there is a water flux that may transport the nutrient load onwards.

2.3. Input Data

NutSpaFHy input contains spatially explicit data on forest resources, soil, topography, water courses and daily weather data. The forest resources input contains the soil type separating the mineral and peat soils, the site fertility class describing the site productivity (sfc, class variable from 1 to 6, highest productivity being 1), and the forest stand characteristics, such as total and species-specific tree volume (m³ ha⁻¹), biomass (kg ha⁻¹), canopy height (m), mean breast height diameter (cm), basal area (m² ha⁻¹) and stand age (years). The input was derived from openly available forest inventory data (FID) maintained by the Finnish Forest Centre (FFC) (available at [6]). The data were organized by forest compartments (stands) that had uniform site characteristics, allowing for accurate simulation of stand development, thinnings and economically feasible regeneration timing. The average forest compartment area in the FID is 1.4 ha. The data also contained the type and timing of the next harvesting operation for each forest compartment. The regeneration-mature stands created a realistic population, from which future clear-cut sites will be selected.

As the FID were mainly available for privately owned forest, the data were back-filled by Multi-Source National Forest Inventory data (MS-NFI) provided by the Natural Resources Institute of Finland. MS-NFI data are based on satellite images, field inventory data and other data sources [61], and are published in 16m×16m resolution. MS-NFI data are comparable to FID, but they are not updated yearly with growth simulation, and do not include the timing and location of the recommended harvesting operations. In our data, 87% of the forest resources data originated from FID, and 13% from MS-NFI data.

The topographic wetness index (TWI) [62] and the terrain slope were calculated using a 16 m × 16 m-resolution digital elevation model (DEM). The DEM, peatland cover, bare rock areas, roads and water elements were obtained from a public database maintained by the National Land Survey of Finland [63]. The Euclidean distance to the nearest water body was calculated for each grid cell. Soil data were derived from the Geological Survey of Finland [64]. All spatial data, if not so originally, were converted into 16 m-resolution with uniform extent.

The weather input was obtained from interpolated weather data in a 10 km×10 km grid between 2002 and 2012 [65]. The data included the daily temperature and precipitation, the global radiation, the relative humidity and the cumulative temperature sum of the growing season. The simulation period represented typical weather in the region, as the mean air temperature, precipitation and reference evapotranspiration (ET0) were close to the long-term median of the period 1981–2018 in all seasons [66]. More detailed information about the input weather is presented in [66].

2.4. Clear-Cut Scenarios

Predicting future nutrient loads is challenging, because the exact clear-cut locations depend on the individual landowners’ decisions, and thus remain largely unknown. To tackle this uncertainty, we formed 100 clear-cut scenarios with equal total harvested areas and unique harvesting locations, using the FID containing the regeneration-mature forest compartments. All the clear-cuts occurred on March 1st of the second simulation year. The total clear-cut area was matched to the realized average annual area from the previous 10-year period. Clear-cuts in mineral and peat soils were tracked separately. The average clear-cut area was 192 ha year⁻¹, of which 83% was performed on mineral soils, and the rest on peatland. The total area of the mature forest compartments, from which the selection of the clear-cuts was performed, was 732 ha, which represented 14% of the productive forest area. The compartments were selected randomly from the potential compartments, and some compartments ended up appearing more frequently in the scenarios than others (Figure 3). Regeneration-mature forest compartments were not evenly located in sub-catchments, because of the heterogeneous forest development class distribution resulting from the long tradition of even-aged forestry.

2.5. Data Analysis

The effect of the location of the clear-cuts on the whole catchment and sub-catchment scales was evaluated by calculating the average annual export (kg year⁻¹) and the relative increase (%) compared to the background export. In this study, ’background export’ refers to the nutrient export from forest areas where no clear-cuts were done during the simulation period. This meant that the background export might include the effects of past forest management. The specific export was calculated by subtracting the background export from the total nutrient export, and dividing it by the area of the clear-cuts. We used these results to quantify how much the clear-cut locations affected the nutrient export, while keeping the total clear-cut area constant in the catchment scale.

When identifying the hotspots in catchments with a long forest management history, we separated the background export from the excess nutrient export caused by the forest management: this was important, because in some cases the nutrient export can be naturally high, because of the catchment characteristics, such as topography, soil types and hydrological and biogeochemical processes [40,57,67]. In other cases, high nutrient export is caused primarily by forest management. Focusing solely on the clear-cut-induced absolute (kg year⁻¹) and relative (%) export increase would neglect the background export hotspots, where the reduction of nutrient export could also be possible; therefore, the role of background export was also studied, to identify the nutrient export reduction priorities. The variation in the background export was presented by exploring the absolute (kg year⁻¹) and relative (kg ha⁻¹ year⁻¹) exports among the sub-catchments.

We used correlation analysis (Pearson) to reveal the influence of the catchment properties on the background export. The variables for the correlation analysis were extracted by calculating the mean or median of the input raster cells (excluding open-area and deciduous tree volume) located inside each sub-catchment, or by using the model outputs. The open area (%) was calculated as a fraction of raster cells where the stand volume was below a threshold value (either 2 m³ ha⁻¹ or 10 m³ ha⁻¹), from the forest area in the sub-catchment. The proportion of deciduous trees of the total stand volume (%) was calculated for each grid cell.

Based on preliminary analysis, we selected two sub-catchments according to different criteria. Sub-catchment 21 was selected because it contained two (out of 100) clear-cut scenarios (A and B) with identical clear-cut areas (17% of the area), but substantially different TN exports. Sub-catchment 26 was selected because it had two clear-cut scenarios with nearly identical TN export, but with different clear-cut areas. These two clear-cut scenarios, in addition to the no-clear-cut scenario (background export) from the two selected sub-catchments, were used to examine the nutrient export within the sub-catchment: this was done by calculating the nutrient balance for each grid cell, to discover the specific locations inside the sub-catchment that produced the highest nutrient export (the top 5% of the area). These locations were named as a within-catchment hotspot.

3. Results

3.1. Clear-Cut-Induced Nutrient Export

Based on the simulated 100 scenarios, the clear-cuts increased the catchment scale TN export on average by 602 kg yr⁻¹ (range 522–677 kg year⁻¹) (Figure 4a). Background export accounted for most of the total export, and the clear-cuts increased this by 4.9% (range 4.2%–5.5%, Figure 4b). The specific TN export varied between 3.0 and 4.1 kg ha⁻¹year⁻¹ in the catchment scale (Figure 4c).

The average clear-cut-induced increase in TP export was 26.0 kg year⁻¹ in the catchment scale (range 22.4–29.5 kg year⁻¹) (Figure 5a): this increase accounted for 5.7% of the TP export (range 5.0%–6.5% (Figure 5b). The specific export varied between 0.12 kg ha⁻¹year⁻¹ and 0.19 kg ha⁻¹ year⁻¹, the average being 0.15 kg ha⁻¹ year⁻¹ (Figure 5c).

There was a substantial variation between the sub-catchments, in terms of how the clear-cuts increased the TN export (Figure 6). The greatest variation occurred in sub-catchment 32, where the total clear-cut-induced excess load varied between 75 and 179 kg year⁻¹ (Figure 6a). However, in relative terms, the variation in this sub-catchment was rather small. For example, in sub-catchment 21, the relative increase in TN export varied from 5% to 69%. Sub-catchments 15, 18 and 21 emerged when the average export increase per hectare was examined (Figure 6c). When examining the range of the specific export, other sub-catchments stood out, particularly catchments 6, 7 and 14 (Figure 6d).

The TP export in the 100 clear-cut scenarios followed the same pattern as that of the TN in the absolute export increase (Figure 7a), the relative export increase (Figure 7b), the hectare-based export increase (Figure 7c) and the specific export (Figure 7d), meaning that the same sub-catchments stood out for the TP export as for the TN export. However, the TP export was not identical to the TN export (e.g., sub-catchment 18, in Figure 6c and Figure 7c) ).

3.2. Background Export

The average annual background export for the TN and TP varied between the sub-catchments considerably (Figure 8). A major part of the variation in the absolute exports can be explained by the forest area of the sub-catchments (9–1036 hectares) (bars in Figure 8). On the other hand, there was a considerable fluctuation in the exports expressed relative to the sub-catchments forest area for both TN (line in Figure 8a) and TP (line in Figure 8b). The average annual background TN export for the whole study area was 2.4 kg ha⁻¹year⁻¹. The difference in the background TN export between the sub-catchment with the lowest (1.4 kg ha⁻¹year⁻¹) and highest export (5.6 kg ha⁻¹year⁻¹) was threefold. The average annual TP export in the study area was 0.08 kg ha⁻¹year⁻¹, whereas the lowest sub-catchment TP export was 0.01 ha⁻¹year⁻¹, and the highest was 0.21 ha⁻¹year⁻¹. The highest relative TN and TP exports did not occur in the same sub-catchments.

We found that the variables that correlated with the TN differed, to some extent, from those that correlated with the TP (Figure 9 and Figure A1). The highest positive correlation with the sub-catchment TN export was the peatland percentage (0.61). Catchment characteristics that were connected to lower TN export included mean distance to the water body (−0.36) and stand biomass (−0.31). For TP export, the highest positive correlations were with the share of the open area (0.62, 0.57) and with the runoff (0.57). Open areas here represented recent clear-cut sites, where wood volume was still low. Stand age (−0.55), fertility class (median sfc) in peat soil (−0.35) and stand biomass (−0.31) correlated negatively with TP export.

3.3. Nutrient Export Hotspots

The background export dominated the TN and TP export (Figure 10). In the majority of the clear-cut scenarios, the clear-cut area did not exceed 10% of the forest area. However, there were considerable differences between the sub-catchments as to how much the clear-cuts increased TN and TP export (Figure 10, dots). In most of the sub-catchments the clear-cut-induced increase was modest, but there were some exceptions, where nutrient export distinguished them from others (e.g., catchments 15, 18 and 21): these sub-catchments could, therefore, be defined as clear-cut-induced hotspots. However, the export increase alone would not necessarily justify the hotspots, if they were reviewed on the absolute scale: consequently, the sub-catchments 15 and 21 would stand out any more, as their background export was low compared to some other sub-catchments (e.g., 27 and 33).

When shifting the focus, in identifying within-catchment hotspots, the results showed clearly, in the selected sub-catchments 21 and 26 (Figure 11), that the nutrient export was very unevenly distributed over the grid cells. The results showed that within-catchment hotspots (top 5% of the area) produced 21%–25% of the TN export and 27%–57% of the TP export, both for background and clear-cut scenarios (

Table 1. Contribution of the within-catchment hotspots (top 5 % of the catchment area) for the TN_kg and TP_kg export sum. TN_% and TP_% are the share of nutrient export of the hotspot areas from the total nutrient export during the 10 year period.

Scenario	Clear-Cut%	TNkg catch.	TNkg hotspot	TN%	TPkg catch.	TPkg hotspot	TP%
Background 21	0	41,236	10,304	25	756	433	57
Clear-cut 21A	17	51,820	11,222	22	803	461	57
Clear-cut 21B	17	62,404	14,810	24	1574	674	43
Background 26	0	45,830	11,540	25	2391	671	28
Clear-cut 26A	4	51,940	10,556	21	2738	739	27
Clear-cut 26B	2	51,864	11,089	21	2811	770	27

).

Figure 11 shows that the increase caused by the clear-cuts was highly location-dependent. In sub-catchment 21 (Figure 11a,b), the cumulative TN export was 20% higher in clear-cut scenario 21B than in 21A (Figure 11b), even though the clear-cut area was the same 17% (Figure 11a). In sub-catchment 26, the nutrient export was similar in clear-cut scenarios 26A and 26B (Figure 11c), even though the clear-cut areas were different.

The range of the TN export at the grid-cell level was similar for both of the clear-cut scenarios in both sub-catchments. In sub-catchment 21, the range of the cumulative TN export was 5–261 kg, while in sub-catchment 26 it was 6–131 kg, during the 10 year simulation period. In both sub-catchments, there were areas with very high TN export, both inside and outside the clear-cut areas. These high nutrient export areas could have arisen from permanent location-specific properties, such as soil type and topography, or recent clear-cut, where nutrient uptake had not yet recovered. In sub-catchment 21, there was considerable difference between the clear-cut scenarios in the location of the clear-cuts within the sub-catchment and, in this case, this clearly affected the total nutrient export (see also Figure 10). Opposite to sub-catchment 21, in sub-catchment 26, even the doubling of the clear-cut area hardly increased the TN export.

4. Discussion

4.1. Impact of Clear-Cut Location on Nutrient Export, on Catchment and Sub-Catchment Scales

The results supported our hypothesis that the location of the clear-cut area, and catchment characteristics, affect the nutrient export. Previous studies have indicated that the increasing proportion of clear-cuts in the catchment increases the nutrient export load [15,41,51]. We demonstrated that the location of the clear-cut affects nutrient export both on the catchment scale and the sub-catchment scale. In the studied catchment, the simulated excess TN and TP export varied between the clear-cut scenarios, although all scenarios had the same total clear-cut area. In our study area, the clear-cuts increased the export load by 4.2%–5.5% for TN and 5.0%–6.5% for TP, with respect to the background load. At first glance, the difference between the lowest and highest export load may seem to be quite modest, but it should be noted that the clear-cut area was less than 4% of the productive forest area, and all the harvesting was done at the beginning of the simulation period. The simulations showed that the difference between the clear-cut scenarios was up to 155 kg for TN export and 7 kg for TP export in the whole study area, whilst keeping the total harvesting area constant. The previously set TP export reduction target to achieve improvement on water quality in Lake Kuonanjärvi (sub-catchment 22) was approximately 110 kg [68]. Compared to the target, locating clear-cuts covered only a small part of the reduction target and, therefore, additional water protection efforts were needed. The results indicated that there was a large variation between the sub-catchments in the clear-cut-induced nutrient export (Figure 6 and Figure 7). Thus, locating efficient water protection measures required shifting focus from the larger catchment to the sub-catchment scale.

The range in the exports between the clear-cut scenarios could be wide, both in absolute (Figure 6a and Figure 7a) and in relative terms (Figure 6b–d and Figure 7b–d). Figure 6a and Figure 7a show that in the two largest catchments (17 and 32) the absolute increase in the export was highest, but the relative increase remained modest.

The export increase in relative terms (Figure 6b and Figure 7b) could be useful to pinpoint those catchments where new clear-cuts might increase the nutrient export considerably in comparison to the background export, and where, therefore, mitigating future clear-cut-induced nutrient export is also an effective strategy on the sub-catchment scale (15, 18 and 21).

In sub-catchments 6 and 7, the additional export caused by clear-cuts (Figure 6d and Figure 7d) can be explained by the fact that the clear-cuts were located in the areas with properties increasing nutrient export. In these sub-catchments, while the relative increase (Figure 6b and Figure 7b) in export may have been modest, there may have been individual clear-cut sites that had a high local impact on the condition of small water bodies (Figure 6d and Figure 7d). On the other hand, the increase in clear-cut-induced export could remain small in such sub-catchments, where the number of stands of regeneration age (sub-catchments 5, 24 and 28, see Figure 2) were low, or where the specific export was low (e.g., sub-catchment 3, Figure 6d and Figure 7d). The TP export was not identical to the TN export (e.g., sub-catchment 18, in (Figure 6c and Figure 7c): this was logical, because the local catchment characteristics had a different impact on the TN and TP exports (see Figure 9, Figure A1 and Figure A2).

Studying the clear-cut-induced increase in the nutrient export on the sub-catchment level revealed clearly those sub-catchments with high export increase potential, but it also showed equally clearly those sub-catchments where the nutrient export reductions by methods associated with clear-cuts (e.g., additional buffer zones, transfer to continuous cover forestry) were not efficient. This information is highly valuable when forming an effective water protection strategy concentrating on using methods applicable to reducing background export (e.g., constructed wetlands.)

When reviewing the results, it should be taken into account that the clear-cut scenarios were formed by randomly selecting the regeneration-mature forest compartments, until the area threshold was reached on the catchment scale, not on the sub-catchment scale: this meant that the clear-cut area in each sub-catchment may have differed between the clear-cut scenarios, even though the total clear-cut area remained constant.

The range in the TN-specific export was 3.0–4.1 kg ha⁻¹year⁻¹ in the whole study area (Figure 4c), which was in the same order of magnitude as those reported from paired catchment studies in Fennoscandia [4] (range 0.08–4.94 kg ha⁻¹year⁻¹). The TP-specific export from the simulations (0.12–0.19 kg ha⁻¹year⁻¹) was also in line with the field measurements (0.01–0.53 kg ha⁻¹year⁻¹) [4] (Figure 5c).

4.2. Background Export

In large forested catchments, the background load is typically the major export component [4,9,33]; therefore, it is important to separate the background export from the excess nutrient export caused by the forest management [33,40,51,69,70]. Here, the background export referred to the export from the observed continuously managed forest area, which was not subject to any forest management during the simulation period.

The results showed that the background load varied greatly between the sub-catchments (Figure 8): this reflected differences in the runoff between adjacent sub-catchments, which has been generally observed to play a role in nutrient export in experimental studies [71]. Our simulated mean annual nutrient exports (TN 2.4 kg ha⁻¹year⁻¹, TP 0.08 kg ha⁻¹year⁻¹) were coherent with the measured values from managed (20 pcs.) and unmanaged (10 pcs.) forested catchments in Finland [33] (TN 2.31 kg ha⁻¹year⁻¹, TP 0.095 kg ha⁻¹year⁻¹). Furthermore, the simulated variation in TN and TP export in this study was well within the measured range [33].

The results suggest that there is no single sub-catchment feature that could explain the variation in TN or TP export (Figure 9): this was expected, as the biogeochemical and hydrological processes regulating nutrient export are highly dependent on local conditions, which can vary between the 16 m×16 m cells even more than between the catchments. For example, tree volume may be a significant factor contributing to nutrient balance on the grid-cell level, but not on the sub-catchment level, as the effect becomes masked with increasing spatial extent. Peatland area had the highest correlation (0.61) with the elevated TN export: this was in line with the results from empirical studies from managed [33] and unmanaged catchments [33,72]. The high TN export from peatland-dominated areas results from the large TN stocks in the peat [73]. The share of the open area (0.62, 0,57) and runoff (0.57) correlated positively with TP export. The strongest negative correlation with TP was with stand age. These TP correlations suggest that the appearance of the TP export hotspot was more dependent on forest management than TN export.

4.3. Hotspots

In this study, the hotspot concept was used in two different contexts. From the sub-catchment populations, we searched for hotspot sub-catchments where the TN or TP exports were exceptionally high: this was justified by the need to prioritize the water protection investments in those sub-catchments with high nutrient export. However, when defining hotspots, it is also important to take into account the condition of the receiving water body, and its ability to receive nutrient export without a considerable negative effect on water quality [74]: this means that there is no precise threshold value for hotspots based on absolute or area-based nutrient export. Because the study area consisted of 33 sub-catchments with different receiving water bodies that had different properties, we neglected the properties of the receiving water body, and determined the hotspot sub-catchment based on the simulated export only (see Figure 6, Figure 7 and Figure 10). Because the WFD also prohibits the deteriorating of the present condition of the water courses, the other premise, besides using the absolute export for hotspot identification, was comparing sub-catchments by their relative (%, see Figure 6b and Figure 7b) or area-based (kg ha⁻¹yr⁻¹, see Figure 6c, Figure 7c and Figure 10, dots) increase in nutrient export caused by forest management.

When identifying hotspots in catchments with a long forest management history, it is important to separate the background export from the excess nutrient export caused by the forest management. In some cases, the nutrient export may be naturally high, because of the catchment characteristics, such as topography, soil types, and hydrological and biogeochemical processes [40,57,67]. In other cases, high nutrient export is caused primarily by forest management. In both of these cases, a small share of the sub-catchment area can significantly contribute to nutrient export (Table 1). Therefore, hotspot analysis could help to improve the efficiency of the water protection measures, when the measures can be located where the nutrient export is actually realized or along their transport paths to water bodies. A topic for future research would be to find methods to assess the costs and effectiveness of the different water protection strategies.

Enhanced water protection with wider buffer zones [40,75,76,77,78,79,80] or the favoring of continuous cover forestry [81,82,83] can be applied to future forest management only and, therefore, separating background export hotspots from management-induced hotspots (Figure 10) is essential. Cost-efficient water protection calls for targeting the investments to areas where the background export is high or will become high as a result of forest management. Therefore, hotspot analysis supports decision making concerning where water protection improvements should be carried out.

To preserve the spatial resolution of the characteristics affecting nutrient export, and thus to maintain the usefulness of sub-catchment-level hotspot analysis, the sub-catchments delineated should be small enough. Based on this study, the sub-catchment area should not exceed a few hundred hectares. The two largest sub-catchments (17 and 32) in this study had a forest area of ca. thousand hectares, which probably should have been divided further into smaller sub-catchments.

The results also indicated that it is useful to extend the hotspot examination on both the sub-catchment and the within-sub-catchment scales. On the within-sub-catchment scale, the water protection planning should be done in conjunction with planning the harvesting operations: this level of planning requires sufficient resolution, to identify areas where clear-cut may mobilize a considerable amount of nutrients. The resolution of 16 m×16 m meets this requirement, and can easily be used together with forest compartment data to guide the width of the buffer-zone (Figure 11,84,85,86,87]), locate water protection structures (e.g., [88,89]) or apply continuous cover forestry [82,90].

In the model description, the export load occurs as a combined effect of biological and hydrological processes. Nutrients become soluble as a result of organic matter decomposition and input from atmospheric deposition. The nutrient uptake is explicitly accounted to the forest stand and to the ground vegetation, and implicitly to microbes, by using the immobilization term in connection with the decomposition. Water fluxes transport the dissolved nutrients along groundwater or surface water flow. Nutrient export emerges in locations where the nutrient release exceeds the nutrient uptake for a prolonged time, and where there is a water flux that may transport the nutrient load onwards. One major source of uncertainty is the actual model description, i.e., does it include the relevant processes regulating nutrient export in reality. In our case, the model included a description of the nutrient storage change calculating nutrient uptake and release by the vegetation in monthly time steps. The transport phase was then calculated, based on the Euclidean distance to the nearest water body or channel. Even though the model was evaluated against empirical data gathered from comparable small catchments, and the use of the calibrated parameters was kept to a minimum (see Section 2.2, [57]), there was still uncertainty as to whether differences between the sub-catchments and within-sub-catchments in the nutrient export were on a scale with the field conditions: this uncertainty could be diminished by actual measurements from the study area. However, we can argue that the results were probably plausible in finding differences in export generation, i.e., hotspots, because the range of the calculated export was in line with independent measurements [4,33,72] and, at the same time, the most influential catchment characteristics impacting the N and P exports were consistent with the published literature (e.g., [33]). This does not eliminate the need to verify the hotspots in the field, especially if within-sub-catchment hotspot maps are included in operational use.

The further development of the NutSpaFHy-model, to directly guide a spatial selection of the water protection measures, would require re-thinking of the model description, so as to also include sediments and their transport, which is important, especially for the calculation of P export.

4.4. Implementing Results in the Water Protection Planning Process

We present a framework for multi-scale water protection planning in two phases (Figure 12), by extending the procedure introduced by Santos et al. [91]. In Phase 1, we define the nutrient export reduction target required for the receiving water body, and we identify sub-catchments where nutrient export reduction is needed. Background and management-induced export are then be compared to the threshold nutrient export defined individually for each receiving water body (e.g., [74]). If the expected nutrient export is smaller than the threshold, the hotspot analysis is not necessary, because no additional water protection investments are needed (Figure 12, Phase 1a). If the export exceeds the threshold (Figure 12, Phase 1b^1–3), the planning advances to Phase 2, where hotspots are identified (hotspot sub-catchments, see Figure 6, Figure 7 and Figure 10). Then, we detect the reason behind the excess nutrient export: if it is caused by forest management (Figure 12, Phase 1b¹), water protection should be directed to sub-catchments with forthcoming forest management operations (e.g., sub-catchment 5 in Figure 12).

If the high nutrient export is due to the background export, more focus should be diverted to background export hotspots (Figure 12, Phase 1b²). In this case, water protection can be enhanced by constructed wetland and surface flow areas situated along the flow paths transporting the nutrients from the hotspots (e.g., sub-catchment 6 in Figure 12). If the forest management has a large impact on nutrient export and, at the same time, the background export is high, a combined water protection strategy is needed (Figure 12, Phase 1b³): this includes both investing in water protection structures along flowpaths, and directing additional investments to water protection of the forest operations (e.g., sub-catchment 2 in Figure 12).

From hotspot sub-catchments (e.g., sub-catchments 2, 4, 5, 6 in Figure 12), actual areas contributing to the nutrient export should be mapped at grid-cell level. The grid-level nutrient export maps are useful for practical forest management planning, e.g., delineating clear-cut areas or finding suitable places for water protection wetlands (see Figure 11). The next step in the development of this approach is to compose methods to assess the nutrient reduction efficiency of different water protection strategies.

5. Conclusions

The results revealed distinctive spatial differences in the background and in the clear-cut-induced nutrient exports. The differences occurred between the sub-catchments, and at grid-cell level within the sub-catchments. Ranking the locations according to the export load allowed for identification of the nutrient export hotspots, from where the nutrient concentration was expected to be high. A very limited area (here, 5% of the sub-catchment area—see Table 1) produced up to half of the total nutrient export. Our state-of-the-art simulation model allowed for composing a new planning framework for water protection in boreal forested catchments. The water protection efficiency of methods based on adsorption and biological uptake of nutrients (e.g., buffer zones, overland flow fields and biochar reactors) typically improves when the concentration of the incoming water is increasing. Better purification efficiency can lead to improved cost efficiency of water protection. If a simulation indicates a particularly high clear-cut-induced export, we could—instead of a clear-cut—apply continuous cover harvesting methods, to maintain a part of the stand nutrient uptake, thus alleviating the formation of the export.