NutSpaFHy—A Distributed Nutrient Balance Model to Predict Nutrient Export from Managed Boreal Headwater Catchments

: Responsible forest management requires accounting for adverse environmental effects, such as increased nutrient export to water courses. We constructed a spatially-distributed nutrient balance model NutSpaFHy that extends the hydrological model SpaFHy by introducing a grid-based nutrient balance sub-model and a conceptual solute transport routine to approximate total nitrogen (N) and phosphorus (P) export to streams. NutSpaFHy uses openly-available Multi-Source National Forest Inventory data, soil maps, topographic databases, location of water bodies, and meteorological variables as input, and computes nutrient processes in monthly time-steps. NutSpaFHy contains two calibrated parameters both for N and P, which were optimized against measured N and P concentrations in runoff from twelve forested catchments distributed across Finland. NutSpaFHy was independently tested against six catchments. The model produced realistic nutrient exports. For one catchment, we simulated 25 scenarios, where clear-cuts were located differently with respect to distance to water body, location on mineral or peat soil, and on sites with different fertility. Results indicate that NutSpaFHy can be used to identify current and future nutrient export hot spots, allowing comparison of logging scenarios with variable harvesting area, location and harvest techniques, and to identify acceptable scenarios that preserve the wood supply whilst maintaining acceptable level of nutrient export.


Introduction
Responsible forest management requires balancing between the economic gains and adverse environmental effects of the production [1]. In Nordic and Baltic countries, forestry has a significant role in national and regional economies but also risks adverse environmental effects. Forest management operations increase nutrient exports to water courses deteriorating water quality especially in headwater streams and lakes [2][3][4][5], but forestry is also a significant source of nutrient load to the Baltic sea [2]. Forest harvesting increases nutrient export through combined hydrological and biogeochemical responses. After the removal of trees, runoff increases [6] and nutrient release exceeds the uptake resulting in increased nutrient transport to water courses [7].
Numerous experimental studies have quantified the effects of forest management on the water quality and nutrient export in headwater catchments [2,4,7,8]. The results are, however, difficult to generalize because nutrient export is affected by several compounding factors. These include, for instance, site fertility and soil properties [9][10][11], topography [12], tree species and ground vegetation [13,14], varying climate and weather conditions [2] and atmospheric deposition [10,15]. Moreover, the size of the harvested area in relation to the size of the catchment, the distance from the receiving water body [4,9,16], dimensions of the buffer zones [8,9], and the intensity and timing of harvesting affect the magnitude and dynamics of nutrient release and export from the terrestrial system to streams.
Novel planning tools are needed to detect areas where forest management may cause increased nutrient leaching, and to choose suitable management practices and water protection methods, such as buffer zones. Such planning tools are currently not available, although different alternatives to estimate nutrient loads from forest management exist. The simplest one is the specific export approach [17][18][19] that computes nutrient loading as a product of managed areas and static export coefficients, which are based on a small number of plot or catchment-scale experiments [17]. As the method lacks spatial context and ignores hydrometeorological variablity altogether, its use is limited to coarse estimates of regional loading and its attribution to different sources. On the other end of the spectrum, there are sophisticated process models on coupled hydrological and nutrient cycle and nutrient transport within stands and catchments [8,20,21]. These models describe water and solute transport and nutrient dynamics at different levels of complexity; there is, however, often no satisfactory means to describe and independently parameterize different processes. As a result, a large number (from tens to hundred) of parameters are commonly calibrated against only a few datasets such as stream runoff and nutrient concentrations, which leads to problem of model equifinality [22] and limits the use of complex models to intensively monitored research catchments.
In the boreal forests, considerable differences in water fluxes, soil moisture and site fertility emerge at scales of kilometres to meters [23,24] following the topographic relief [25]. The water and nutrient availability are the decisive factors for determining forest productivity and the scale of response under forest management operations [26]. Based on this rationale, topographic wetness indices are used to identify wet areas in the landscape, improve multi-criteria forest planning [26] and optimize location and dimensions of buffer zones [27]. While already being implemented in forest planning systems in Finland and Sweden [28,29] the approach remains as a static and qualitative rank of moisture conditions within the catchment.
In forest management planning, the traditional scale has been the individual stand or set of stands managed by a single forest owner. However, from a water quality perspective, a catchment forms the basic unit to quantify the impacts of forestry as water and nutrients released from a site are transported to the downstream watercourse. To bridge these two scales, we developed a distributed nutrient balance and transport model, NutSpaFHy. It extends the recently-proposed Spatial Forest Hydrology model (SpaFHy, [30]), and is designed as a robust and practically applicable tool to quantify nutrient export from managed forested catchments. We test the model's predictive capability against measured streamwater nitrogen (N) and phosphorus (P) concentrations and export loads at 18 headwater catchments across Finland. After testing, we demonstrate the model in a practical forest planning at a single catchment. We analyzed how location, extent and intensity of forest harvesting affects catchment-level nutrient export and discuss the practical implications of our approach. As NutSpaFHy uses only openly available data as input, and contains only four calibrated parameters, it it is designed to be applicable for whole Finland and extendable to other Nordic and Baltic countries.

Field Data
We used data from 18 forest-dominated headwater catchments covering a large geographical area across Finland (Figure 1). Twelve of the catchments were used for model calibration and six as independent test catchments. The area of the catchments ranged from 31 to 1966 ha ( Table 1). Five of the catchments were pristine or outside forest management, while 13 were under normal forest management, which typically includes growing of even aged semi-natural stands treated with thinnings and having rotation time from 70 to 120 years. The dominant tree species in the catchment were Scots pine (Pinus sylvestris L.) and Norway spruce (Picea abies Karst.) and the catchment mean stand volume varied between 25 and 166 m 3 ha −1 . The distribution of site fertility classes in peatland and upland sites are presented in Table A1. Runoff at catchment outlets was continuously measured from 2014 to 2016 using the V-notch weirs equipped with automatic water level loggers. On average, 20 water samples (data available: https://metsainfo.luke.fi/fi/vesistokuormitukset (accessed on 6 June 2021)) were collected per year and total nitrogen (N) and total phosphorus (P) were determined at certified laboratories using accredited methods (SFS-EN 45001:1990, SFS-EN ISO/IEC 17025:2000, SFS-EN ISO/IEC 17025:2005) [31].
Daily runoff and linearly interpolated concentration time-series were used to calculate the monthly and annual export loads. Table 1. Characteristics of calibration and test catchments: id denotes catchment identification number, p is pristine and m is managed catchment, Area is catchment area in ha, T sum is temperature sum in degree days, P r is mean annual precipitation in mm, V is the mean stand volume in m 3 ha −1 , f coni f is the fraction of coniferous trees from the total stand volume in the catchment, Slope is the mean slope in the catchment in %, D water is the mean distance to water body in m, and Peat is the share of peatland sites from total area of the catchment.

Model Description
NutSpaFHy is a spatially-distributed nutrient balance model, which extends the hydrological model SpaFHy [30] by introducing grid-based nutrient balance sub-model and a conceptual solute transport routine to approximate total N and P export load to streams ( Figure 2). A detailed description of the model is provided in the Appendices.
The model domain consists of a rectangular grid (cell size 16 m) that is initialized using open-source GIS data available throughout the Finland [30]. These include the Multi-source National Forest Inventory (MNFI) data for initial forest structure, soil maps and topographic databases for determining soil type (fine, medium or coarse mineral soils and peatlands) and permanent water elements to locate streams and water bodies. The digital elevation model (DEM) is used to compute local slope, topographic wetness index (TWI) [30,32] and flowpath distances from a grid cell to the nearest stream. Basic daily meteorological data (global radiation, air temperature, humidity and precipitation rate) are used to drive the SpaFHy-model, which produces daily water table (WT), root layer water content (W), water flux from root layer to deeper soil (Q dr ), water flux to root layer from below (Q ex ), and runoff (Q runo f f ). For NutSpaFHy, the water fluxes and state variables are further aggregated to monthly values. Prior to the NutSpaFHy simulation, a database of development parameters for stand height and yield are computed using Motti-simulator [33]. We treat separately growth and yield of Scots pine (Pinus sylvestris L.) and Norway spruce (Picea abies Karst.) growing on peatland and mineral soil sites in site fertility classes (s f c) 2, 3 and 4 (growth decreases with increasing site fertility class, [34]). The initial state of the forest stand in each grid cell is defined by the MNFI data that contains information for the main tree species (spe), s f c and site main class (smc, peat/mineral) and the mean stand height (h g ), volume (V) and age (A). Thereafter, the stand growth follows the development parameters and Equations (A1) and (A2). The growth rate is adjusted according to the annual temperature sum, which incorporates a simple dependency between the growth performance and weather conditions (A13). The growth of stand determines the net nutrient uptake as described by Palviainen and Finér [13] (Equation (A7)). To obtain the gross nutrient uptake we have to add the nutrient uptake by ground vegetation, and the nutrients that are lost in the litter fall as described in Laurén et al. [35]. The total monthly uptake of N and P (N up , P up ) were given as input to soil nutrient balance module.
Decomposition of organic matter in mineral soils and peatlands is described using empirical equations presented by Pumpanen et al. [36] and Ojanen et al. [37], respectively. We use monthly mean air temperature, W m and WT m as drivers of decomposition rate (Appendix B.2.3). During decomposition, a fraction of the released N and P is immobilized to microbe biomass (Equations (A16) and (A20)) as in [35]. Immobilization parameter (Imm peat NP and Imm miner NP ) values in NutSpaFHy were first calibrated using 12 catchments (see Section 2.3), and then a regression model was built to predict these parameters from the catchment properties. The net release of N and P (r min NP , r peat NP ) were given as input to nutrient balance module.
The nutrient balance module keeps an account of the total nutrient storage and nutrient concentration in the root layer. At each timestep, the nutrient storage is updated accounting for the nutrient release and uptake and nutrient concentration in the liquid water (W m ) in the root layer is computed. Using the monthly mean nutrient concentrations and monthly cumulative drainage (Q dr m ) and surface runoff (Q ex m ), we calculate the respective nutrient fluxes to groundwater (Q dr NP ) and as a form of surface runoff (Q ex NP ).
The nutrient transport via groundwater flow from a grid cell to the receiving water body is delayed by the residence time (t delay , in months), and reduced by an empirical retention factor that depends on the grid-point distance to the receiving water body [38] (Appendix B.3). The t delay is computed for each grid cell using the slope and distance to the receiving water body. The resulting variable, Q dr NP delayed , represents the input to the groundwater nutrient storage. The output from the groundwater nutrient storage is nutrient flux with runoff. Nutrients in the surface runoff are transported to the receiving water body without a time delay and no retention processes are considered. The total nutrient flux at the catchment outlet is the sum of the groundwater and surface runoff transport from all grid cells. NutSpaFHy extends the hydrological model SpaFHy [30] (yellow box). Blue boxes refer to model inputs, orange box to a priori computation of regional forest growth parameters, green boxes are computation modules, grey boxes are variables or data transferred between the modules, and purple box is the model output. Outlines of the boxes describe the frequency of the computation: The box with double solid outline is computed a priori, boxes with a wide dashed outline are calculated only in the beginning of the simulation, and boxes with a square dotted outline are computed in monthly time-step. SpaFHy (yellow box) is calculated in daily time-step. Variable symbols are provided in Abbreviations. Symbols B.1. . . B.3 refer to Appendix B, where the submodules are described in detail.

Calibration
Nutrient export load is strongly dominated by runoff, which is modeled using SpaFHy. The hydrological predictions (runoff, evapotranspiration, snow and soil water storage) of SpaFHy have been thoroughly tested against the observed data in [30]. To calibrate the new processes in NutSpaFHy, the monthly mean observed N (c N obs ) and P concentrations (c P obs ) in the stream water from the twelve catchments were used ( Figure 1, Table 1). Optimization was done separately for N and P as there is no interaction between them in the model. For both N and P we had two calibrated parameters, the immobilization factors for mineral and peat soils (imm N peat , imm Nminer and imm Ppeat , imm Pminer ). Catchments were calibrated individually in a process, where, in each iteration of the optimization, the observed monthly concentrations were compared to the predicted monthly concentrations using regression through the origin (Equation (1)).
where c obs m was the observed concentration of N and P in month m, c pred m was the predicted concentration of N and P for month m, and s is a slope parameter, and ε m is the residual term for month m with expectation value of zero. When s equals to one, NutSpaFHy predicts the nutrient concentration with no bias. Thus, the minimized target variable in the optimization was (s − 1) 2 . Initial values for all calibrated immobilization parameters was set to 0.9, and the valid range was for the parameter values was set to [0.5, 1.0]. The optimization was conducted with minimize function in scipy.optimize-package in Python 3.7. To support the model use for ungauged catchments, we investigated whether the immobilization parameters ( imm N peat , imm Nminer and imm Ppeat , imm Pminer ) can be predicted from catchment characteristics. Using the calibration catchments, we predicted the optimized immobilization parameter values and the catchment characteristics in linear stepwise regression. The backward step-wise regression drops the start from the full model and drops the weakest explaining variables until the model is not improved anymore. Forward stepwise regression starts with a single explanatory variable and adds variables until no improvement is achieved. The used olsrr package in R applies a combination of backward and forward regression analysis. (Tables 1 and A1, Appendix A).

Model Testing
The NutSpaFHy model was tested against independent data from six catchments ( Figure 1, Table 1). The immobilization parameters were predicted from the catchment characteristics and the root mean square error of the modeled immobilization parameter was used to create 95% confidence intervals for the nutrient concentration and export load predictions (Appendix A). The simulated monthly and annual concentrations and export loads with their confidence intervals were compared to observations.

Application to Clear-Cut Scenario
To test the NutSpaFHy in realistic forest planning situations, we simulated alternative harvesting operations at catchment 2. We compare N and P export loads from 25 alternative clear-cut scenarios to a reference scenario (uncut) where no management was done ( Table 2). To follow common forest management practice, we applied clear-cut only at grid cells where the stand volume exceeded 100 m 3 ha −1 . The scenarios contained alternative harvest regimes, where similar size clear-cut area was located close and far away from the water (<35 m, >100 m), where similar total tree volume was harvested from peat and mineral soils (Peat < 100 m, Min < 100 m), and where total volume was harvested from low fertility (s f c 4,5,6) vs. high fertility sites (s f c 1,2,3). In addition, we located clear-cuts extending from stream to 10 . . . 190 m distance (<10 m . . . <190 m), and finally we calculated 10 scenarios, where the clear-cuts with same harvested total volume were randomly located around the catchments (Rand 1 . . . 10).
The results were expressed both as mean annual export loads (per unit catchment area) during the first 10 years after harvesting, and as specific export loads, where the increase of export load compared to the reference scenario was divided by the fraction of the harvested area from the total catchment area. We selected Catchment 2 for this analysis because it allowed the compilation of versatile management scenarios due to the high mean stand volume distributed quite evenly on both mineral soils and peatlands. These kinds of catchments are prone to high nutrient export loads under intensive forest management.

Model Calibration and Immobilization Parameters
The mean (and range) values for imm N peat , imm N min , imm P peat , imm N min were 0.88 (0.79-0.94), 0.92 (0.84-0.96), 0.92 (0.78-0.97), 0.92 (0.83-0.98), respectively. The N immobilization parameters were found to depend on catchment characteristics, whereas P immobilization was best explained using the mere average over the calibration catchments (Table 3). N immobilization in peat soils increased with an increasing proportion of coniferous tree volume from the total tree volume in the catchment, and decreased with an increasing share of bogs (Appendix A). N immobilization in the mineral soils was increased by the share of low-fertility mineral soils in the catchment. Table 3. Nitrogen and phosphorus immobilization parameters as a function of catchment characteristics calculated using stepwise regression analysis. Catchment characteristics are shown in Tables 1 and A1.

N and P Concentration and Export Load
The highest monthly export loads of N and P typically emerge during the spring floods, which were well captured by NutSpaFHy (Figures 3 and 4). For the calibration catchments, both the mean concentration and the range of seasonal concentration variation were rather well-captured, especially when the immobilization parameters were calibrated for the specific catchment ( Figure 3). However, the timing of the concentration fluctuation was rather poorly predicted. As expected, the model performance decreased when parameters were predicted from catchment characteristics without calibration (Figures 3 and 4). However, in most cases the range of seasonal concentration variation was still rather well as the observed concentrations fit within the model confidence intervals. According to R 2 and RMSE the export loads were better predicted than the concentrations.
The observed annual N export varied from 0.5 to 5 kg ha −1 year −1 , and the P export from 0.1 to 0.25 kg ha −1 year −1 ( Figure 5). For the test catchments, NutSpaFHy slightly overestimated the annual export loads both for N and P, whereas the annual N concentrations were remarkably well predicted ( Figure 5).

NutSpaFHy Application
To test the NutSpaFHy in realistic forest planning situations, we compared N and P export loads from alternative clear-cut scenarios to a reference scenario (uncut) where no management was done ( Table 2). In the uncut scenario, the N and P export were 2.0 and 0.08 kg ha −1 year −1 , respectively ( Figure 6). Locating a similar size of clearcut area close to receiving water body (scenario < 35 m) resulted in considerably higher specific N and P export loads than clear-cuts located further away (>100 m) (Figure 6b,d).
Harvesting similar tree volume close to the waterbody from peatlands (scenario Peat < 100 m) resulted in higher specific N export, but lower P load than harvesting from mineral soils (Min < 100 m). Clear-cuts in fertile sites (s f c 1,2,3) led to clearly higher specific exports than in lower-fertility sites (s f c 3,4,5). Scenarios <10 m . . . <190 m describe a situation where the clear-cut area is extended gradually from close to water body towards further-away locations. The total export load (Figure 6a,c) increased with increasing clear-cut area and harvested volume ( Table 2), but the specific export decreased when including harvests from further-away grid cells (Figure 6b,d). Results from Rand 1 . . . Rand 10 indicate low variability in export loads when half of the catchment area was clear-cut from randomly selected grid-cells.  Table 2.

Model Requirements
We built a distributed nutrient balance model NutSpaFHy to simulate forest harvesting impacts on N and P export to watercourses at forested headwater catchments. NutSpaFHy is aimed as a robust and practically applicable tool to assist forest management planning, and thus the model development was guided to meet the following requirements: (1) Modular structure to allow future development; (2) Use only open data for model calibration and applications; (3) Produce realistic nutrient exports in a versatile stand, site and soil combinations; (4) Describe how nutrient export depends on spatial harvesting patterns in the catchment; (5) Respond to climate gradient across Finland, as well as on seasonal and inter-annual variability of meteorological conditions; (6) Have a low number of calibrated parameters to be applicable outside specific research catchments. Moreover, the model code must be distributed under open-source license to stimulate its use and future development.
Meeting the above constraints is challenging, as they mean prioritizing data availability and model applicability over detailed descriptions of biogeochemical and hydrological processes. Thus, the structure of NutSpaFHy reflects balancing between the practical demands and constraints, and the scientific rigor to present the processes regulating nutrient release and transport in boreal ecosystems. Most importantly, we included only the main nutrient in-and outflows, storages, and transport mechanisms in the model, and focused especially on describing those processes and fluxes that are directly affected by forest management practices.

Evaluation of Model Structure
NutSpaFHy consists of several modules ( Figure 2). The hydrological model SpaFHy has been rigorously tested and successfully applied for a large range of different catchments [30]. Based on open data, SpaFHy is particularly good in describing daily grid-cell water fluxes and catchment runoff in heterogeneous catchments. In SpaFHy, grid-cell hydrology responds to local stand characteristics and can well describe the response of these fluxes to leaf-area and species composition changes due to forest management [30]. Drainage from the root zone or return flow from the groundwater link the grid-cell water balance to catchment-level groundwater module. The latter is, however, described using the classical Topmodel-approach [39] that does not explicitly account for water movement from a grid cell to another. Therefore, we created a separate nutrient export module to account for the delay and retention during the export.
Nutrient release in organic matter decomposition and atmospheric bulk deposition are the largest inputs of N and P to soil water [8], whereas the largest outputs are nutrient uptake by stand [13], ground vegetation [14] and soil microbes [40]. In contrast, nitrification, denitrification and N 2 O emissions are considerably smaller fluxes and omitted from NutSpaFHy as they are difficult to describe and parameterize at the relevant spatiotemporal scales using available data. As organic fractions of N and P dominate in the export loads in boreal forested catchments [18,41], we accounted only for total dissolved nutrients that include both inorganic and organic forms of N and P. This considerably simplifies the model structure, and it is the total nutrient export which determines the water quality because organic N and P are eventually mineralized or photochemically degraded in the water courses [42,43].
Microbial immobilization is an important biochemical nutrient flux in soils [40], and it has been considered to be the primary process affecting especially N leaching [44]. We calibrated N and P immobilization parameters against the observed N and P concentrations in the runoff water. The mean values for all the four immobilization parameters were close to 0.9, which equals the value used in the decompostion-the nutrient dynamics model Romul [45] for early-stage decomposition. A high immobilization rate is expected because organic matter C/N ratio is typically high in boreal soils [46]. Modeling of stream N concentrations was slightly more consistent than P, and the imm N values could be reasonably well predicted from the catchment properties (Table 3). It is widely known that accurate modeling of P dynamics in terrestrial ecosystems is difficult [47,48] and requires data that is not available at spatial scales relevant for distributed catchment models. This is especially because processes controlling P kinetics and transport are closely related to biogeochemical and hydrological processes such as soil redox conditions and the presence of iron and aluminum in soil [49][50][51].
The growth of trees and understory vegetation is the primary driver for nutrient uptake. The initial state of the forest stand in NutSpaFHy is described based on reasonably accurate gridded forest data [52]. The stand growth onwards from the initial state is described using statistical growth curve parameterized for each site type a priori, and the monthly nutrient uptake follows the prevailing weather conditions in a simple way. Furthermore, the development of understory vegetation follows the stand development. It is clear that the tree stand development is too simplistic to make NutSpaFHy a growth and yield simulator; however, the current stand description is flexible enough to produce reasonable estimates for nutrient uptake for versatile sites, stands and weather conditions. Forest harvesting, such as thinning or clear-cutting affect the nutrient balance through decreasing the nutrient uptake and increasing runoff, and consequently more nutrients are leached. Decomposition models in NutSpaFHy [36,37] are robust empirical models which are responsive to prevailing temperature, soil moisture, fertility and forest stand volume. However, the current versions of the decomposition models do not include the input of logging residues into the soil and the consequent increase in the organic matter decay. A future development of NutSpaFHy should include a decomposition model that accounts for the organic matter mass balance in the soil.
The time delay and the nutrient retention during the transport are the most important catchment-level factors that regulate the dynamics and magnitude of nutrient load. The time delay and the nutrient retention can be computed in physically-based solute transport models [8,53,54], where diffusion-advection-dispersion equation is solved in two or three dimensions. These models are computationally very demanding and the parameterization requires detailed soil data that is not available outside of research areas. In NutSpaFHy, we used pre-calculated residence time to describe the time delay and an empirical model to account for the nutrient retention. This approach is simplistic, but computationally efficient. It is unlikely that this approach is capable of producing daily nutrient export dynamics; however, in forest planning the export loads become meaningful only in much longer (i.e., annual and decadal) time-scales.

Model Performance at Study Catchments
In the experimental datasets, we had catchments with different soil types and fertility distributions, climatic conditions, management history and a large range of stand properties. NutSpaFHy was able to produce realistic estimates for seasonal dynamics of N and P export (Figures 3 and 4) and mean annual concentrations and export loads ( Figure 5). However, the timing in monthly concentrations was more inaccurate, probably because of the simplified computation of residence time and nutrient retention (Appendix B.3). The time scale in forest management planning is, however from years to decades, and therefore the mismatch in monthly concentrations is not detrimental.

Nutrient Export from Different Harvesting Scenarios
The NutSpaFHy was tested for different harvesting scenarios at a single catchment. The model predicted higher export loads from harvests that were located close to streams than from those located further away. This is consistent with several field studies showing that uncut buffer zones between the clear-cut and the stream can efficiently reduce nutrient loading to water courses [8,9,55]. Moreover, NutSpaFHy predicts higher specific N export loads from peatland than from mineral soil clear-cuts, which has been documented in several field studies [2,5,56]. In contrast, higher specific P loads were detected from mineral soil than peatland clear-cuts ( Figure 6). This most likely reflects the relation between the calibrated imm P peat and imm P min raising from the calibration data, which contained both undrained and drained peatlands. Finér et al. [2] used partly same dataset than us, and found that P export load increased for soil types in the following order: undrained peatlands < mineral soils < drained peatlands. Both drained and undrained peatland grid-cells were treated similarly in our calibration procedure, which may have resulted in lower modeled specific P export loads for peatlands than mineral soils.
According to the clear-cut scenarios, the nutrient export from harvests on fertile sites was higher than that from low fertility site clear-cuts. This reflects primarily differences in the quality and nutrient concentrations in the organic matter, which enables faster decomposition and higher nutrient release from the fertile sites [57][58][59][60].

Potential for Forest Management Planning
Implementation of the Water Framework Directive (WFD) in commercial forestry has raised a need to quantify the impact of various harvesting scenarios on nutrient exports from catchments. Accounting for the landscape heterogeneity is essential in the harvest planning. Using the NutSpaFHy model with up-to-date forest inventory data, and including spatial information of timing and type of planned harvests, it is possible to estimate the future increase in the nutrient export caused by the logging. The model can be used to identify current and future hot spots of nutrient export in the catchments. This allows testing and comparison of logging scenarios with variable harvesting area, location and harvest techniques to identify acceptable scenarios which preserve the wood supply whilst maintaining an acceptable level of nutrient export.
We argue that the current state of nutrient loading (and its predicted future dynamics) at the actual catchment of interest, not the background load from an imaginary pristine forest, should set the baseline for comparing the impacts of any future forest operations. Calculating export loads using up-to-date forest data with no new loggings planned, we can estimate how the past logging history affects the present and future 'background' load at a given catchment. Thus, the NutSpaFHy also facilitates setting more realistic goals for what can be achieved by optimising future commercial forestry operations.
In the Anthropocene, the meteorological conditions regulating biogeochemical processes and nutrient cycle in forests are expected to change. This is commonly expected to increase nutrient leaching from forested areas compared to the present climate; however, increase in nutrient loads may differ between catchments and regions. Running NutSpaFHy using climate model predictions can help to identify catchments and within-catchment locations particularly vulnerable to climate change.

Conflicts of Interest:
The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Catchment Properties
Various catchment properties derived from the open GIS data were tested in predicting the calibrated immobilization parameters using stepwise regression. Bogs, fens and open peatlands are identified in the MNFI data as site main classes (smc) 2, 3, and 4, respectively (Table A1). For more detailed description of smc see [34]. Stand productivity in peatland and upland mineral soil decreases when the site fertility class (s f c, [34]) in MNFI increases from sc f 1 to s f c6. We combined s f c1 and s f c2 as rich peatlands (p rich ) and mineral soils (m rich ), s f c3 and s f c4 to medium fertile peatlands (p medium ) and mineral soils (m medium ), and s f c5 and s f c6 to poor peatlands (p poor ) and mineral soils (m poor ). The areal share of grid cells with these characteristics were used in the stepwise regression analysis.

Appendix B. Description of NutSpaFHy
NutSpaFHy is a grid-based, semi-distributed model for predicting the impact of forest management on total N and P loads to surface waters. It uses variety of datasets and contains sub-modules that represent regional-scale, catchment-scale and grid-cell scale. The sub-modules are described below.
Appendix B.1. Regional Scale Growth and Yield Forest growth from the initial state to the end of the simulation is required to simulate the nutrient uptake. Forest growth pattern, i.e., the shape of height and stand volume development curves, varies between tree species and geographical regions, and is distinctively different for mineral and peat sites. Growth can conveniently be described with scalable anamorphic growth curves [61] (Equations (A1) and (A2)) that describe stand development at a grid-cell provided that the shape parameters are known. Therefore we applied, a priori, Motti-simulator [62] to model growth and yield patterns of forest along a regional grid over Finland (12 grid points). The simulation was done for medium fertility upland mineral soil and peatland sites with Norway spruce (Picea abies Karst.) or Scots pine (Pinus sylvestris L.) as a dominant species (4 simulations/grid point). A whole rotation from regeneration to harvesting was computed using the recommended thinning regime in practical forestry [63]. We recorded stand age (A, years), mean height (h, m), and yield (y, m 3 ha −1 ) from each simulation. The development of h and y follows Schumacher equation [64,65], whose shape parameters were defined through a curve-fitting procedure to the Motti-simulated results (Equation (A1)).
where h (m) is stand mean height at age A (years), h i (m) is the mean stand height at index age A i , p 0 and p 1 are fitted shape parameters. Index age was chosen to be the median age in the simulation time series. The curve defined by Equation (A1) passes through point (x = A i , y = h i ), therefore changing h(A) scales the curve and produces a family of anamorphic (same shaped) curves, which allow a convenient method for describing forest height development for different site productivity using the same shape parameters and one observed A, h point [66].
Similarly y is where y(A) (m 3 ha −1 ) is yield at age A (years), y i (m 3 ha −1 ) is the yield at index age A i , p 2 and p 3 are fitted parameters. Parameters p 0 , p 1 , p 2 and p 3 were indexed by grid point coordinates, site main class (upland, peatland) and tree species and saved into a database. When a grid-cell is clear-cut in the model, an immediate regeneration is assumed, and the age of the stand is set to one year and the stand development follows the same initial growth curve. Hydrology is calculated using spatially-distributed SpaFHy model [30], which produces daily rooting zone water content (W, m 3 m −3 ) and local soil moisture deficit at grid scale to be used as an input for the nutrient balance modelling. SpaFHy also accounts for daily water flux from the rooting zone down to groundwater (Q dr , m day −1 ), daily return flow from groundwater to rooting zone (Q ex ) in grid scale, and the daily volume of runoff (Q runo f f , m 3 day −1 ) in catchment scale. For NutSpaFHy, the water fluxes and state variables were aggregated to monthly timescale.

Appendix B.2.2. Nutrient Uptake
Nutrient uptake at each grid cell was estimated based on predicted stand yield. First, the site main class and tree species were obtained from the MNFI data for each grid point. Based on this information, the shape parameters p 0 and p 1 were retrieved from the Mottisimulation database. Then the observed stand age (A obs , years) and the stand mean height (h obs , m) in the grid point were retrieved from the MNFI data. Next, A obs and h obs were substituted to index age A i and height at index age h i and in Equation (A1) resulting in an anamorphic growth curve for the grid point. A relative growth performance (g rel ) describing local difference in productivity: where h obs (m) is the observed stand mean height at age A obs (years), and h Motti is the height at age A obs using a priori computed Motti parameters. The stand yield (y(t sim ), m 3 ha −1 ) during the simulation period (t sim ) was computed as: where Equation (A2) is substituted by observed age (A obs , years), t sim is the duration of the simulation (years), p 2 , p 3 , A i , y i are a priori computed parameters from Motti-simulations. The stand volume (V end , m 3 ha −1 ) in the end of the simulation period is: where V obs is the observed initial stand volume from MNFI data (m 3 ha −1 ), and y(t sim ) is the projected yield during the simulation period (m 3 ha −1 ). Nutrient storage of the tree stand can be computed from the dominant tree species (s p ) and the stand volume according to Palviainen and Finér [13]: where C N,P is content of N and P in stand (kg ha −1 ), V is stand volume (m 3 ha −1 ), lna, b, and k are species-specific parameters provided by Palviainen and Finér [13]. The net N and P uptake (U net N,P , kg ha −1 ) resulting from stand volume growth during the simulation period can now be calculated as: U net N,P = e (lna+b * ln(V end )+k) − e (lna+b * ln(V obs )+k) . (A7) The amount of nutrient m that the tree stand loses in the litterfall must be compensated with uptake from the soil (U comp N,P , kg ha −1 ). It is calculated accounting for the nutrient retranslocation from senescing leaf mass: where LAI obs is the observed leaf area index (m 2 m −2 ), c lea f N,P is concentration of N, P in leaf (kg kg −1 ), ret N,P is share of nutrient re-translocation in litterfall, sla is specific leaf area (kg m −2 ), llo is leaf longevity (yrs), and t sim is length of simulation (years). Ground vegetation nutrient uptake was estimated from the net change in the ground vegetation biomass (Equation (A9)) and the uptake needed to compensate the nutrient loss in the litterfall (Equation (A10)) as presented in [35]. Grid cells where the site main class (smc) was 1 were classified as upland sites and grid cells where smc>1 were classified as peatland sites. Thereafter, we used the dominant tree species to match the cell with a correct empirical ground vegetation biomass model presented by Muukkonen and Mäkipää ([34]; Tables 6-9) and solved the ground vegetation biomass in component i (gv i , kg ha −1 ). Ground vegetation was considered in i components, where i = [dwarf shrubs, herbs and sedges, mosses]. Field layer consists of dwarf shrubs and herbs and sedges, whereas the ground layer consists of mosses only. In upland mineral soil sites, the share of dwarf shrubs and herbs from the total field layer biomass was assumed to be 91%, 9%; 71%, 29%; 38%, 62% in Scots pine, Norway spruce and broad leaved stands, respectively ( Figure 1 in [34]). Parameterization of peatland sites is presented in details in [35]. Ground vegetation nutrient contents (cgv N,P , mg g −1 ), turnover rate (τ, years −1 ), and nutrient retranslocation fractions (rgv N,P , kg kg −1 ) (Table A2) were derived from [14,67,68]. Table A2. Ground vegetation nutrient contents (cgv N,P , mg g −1 ), turnover time (τ, years −1 ), and nutrient retranslocation fractions (rgv N,P , kg kg −1 ).
where Ugvnet i is the total net N, P uptake of ground vegetation component i (kg ha −1 ) during the whole simulation period, ∆gv i is the ground vegetation biomass change from the beginning to the end of the simulation (kg ha −1 ), cgv i N,P is the nutrient concentration in the ground vegetation component i (mg g −1 ). Change in gv i occurs when the stand volume, stem number, basal area and stand age change from the beginning and the end of the simulation. In some cases with increasing stand volume ∆gv i may become negative, and then Ugvnet i N,P is set to zero.
where Ugvlitter i N,P is the annual nutrient uptake needed to compensate the nutrients lost in the litterfall (kg ha −1 year −1 ), τ i is the turnover rate of ground vegetation component i (year −1 ), and rgv i N,P is the fraction of N, P retranslocated back to living tissues before the litterfall (kg kg −1 ). Total ground vegetation nutrient uptake (U gr N,P kg ha −1 ) was obtained by: where Ugvlitter i N,P is evaluated at the beginning of the simulation (t 0 ) and at the end of the simulation (t end ), and t sim is the length of simulation (years). The combined N,P uptake (U N,P , kg ha −1 ) of the stand and ground vegetation was calculated as a sum of stand and ground vegetation uptake: U N,P = U net N,P + U comp N,P + U gr N,P . (A12) Uptake of N and P (U N,P m , kg ha −1 month −1 ) for month m was derived from U N,P tot using the temperature sum for month m (T sum m , degree-days), and the long-time temperature sum (T sum , degree-days): This approach allows a simple dependency between the nutrient uptake and the weather conditions throughout the simulation period: a period with a higher temperature sum than the long-time mean results in a higher nutrient uptake than a colder period.
When a grid-cell is clear-cut in the model, the ground vegetation is adjusted according to the new stand properties.

Appendix B.2.3. Nutrient Release
Nutrient release is computed separately for upland mineral soils and for peatland soils, which is directly affected by heterotrophic soil respiration, describing the rate of organic matter decomposition. For mineral soil sites, soil respiration was calculated as [36]: r min CO2 = r min 0 * f moist * Q T soil −10 10 10 n days in month , where r min CO2 is heterotrophic soil respiration in CO 2 (kg ha −1 month −1 ), r min 0 is reference rate of heterotrophic respiration in CO 2 (60.82 kg ha −1 day −1 ), Q 10 is temperature sensitivity (value 2.3), and T soil is monthly mean soil temperature ( • C, here T soil is assumed to be minimum of monthly mean air temperature and 16.0), n days in month is number of days in month (days month −1 ) and f moist is moisture restriction function [69]: where W m is liquid water content in root zone as simulated by SpaFHy (m 3 m −3 ), poro is soil porosity in the root zone (m 3 m −3 ). To obtain release of N,P (r min N,P , kg ha −1 month −1 ) in the organic matter decomposition, we computed: where conv CO2toC is conversion factor from CO 2 to elemental C (12 kg mol −1 /44 kg mol −1 ), c om N,P is N,P content in soil organic matter in mineral soils (kg kg −1 ). According to Tamminen [60], N content in soil organic matter of s f c 1, 2, 3, 4 and 5 is 0.024, 0.022, 0.018, 0.016 and 0.014 kg kg −1 , respectively. P content of s f c 1, 2, 3, 4 and 5 is 0.0017, 0.0015, 0.0013, 0.0011 and 0.0010 kg kg −1 , respectively. c om C is C content in organic matter (0.55 kg kg −1 ), and imm miner N,P is fraction of the released N and P that is immobilized into microbial biomass (kg kg −1 ). imm miner N,P was calibrated against measured N and P concentrations in the runoff water (see Section 2.3. Calibration). Nutrient release in peat soils followed the same principle. First the CO 2 efflux was calculated as [37]: where r peat CO2 is heterotrophic peat respiration (kg ha −1 month −1 ), r peat 0 peat heterotrophic respiration (kg ha −1 day −1 ) in reference soil temperature T soil re f (10 • C), T soil is monthly mean soil temperature at depth of 0.05 m ( • C), here represented with monthly mean air temperature T soil = min (T, 16.0), T soil 0 is soil temperature at which soil respiration is zero (−41.02 • C), B is a parameter, and n days in month is number of days in month. The peat respiration in the reference temperature depends on stand and site properties, and water table [37]: where V is forest stand volume in the grid point (m 3 ha −1 ), ρ b is peat bulk density (kg m −3 ), and WT g s is mean water table during the growing season (here we use monthly mean water table WT m ), and conv 1 is conversion factor from g m −2 h −1 to kg ha −1 day −1 . Parameter B is obtained as: where T air gs is mean growing season air temperature ( • C), d peat is peat depth (here set to 0.99 m). MNFI data on peatland site fertility class (s f c) was used together with reported field data [59] to assign the grid-point peat characteristics: ρ b for s f c 1. . . 5 was: 140, 140, 110, 100, 80 (kg m −3 ); c peat N was: 0.019, 0.019, 0.016, 0.014, 00012 (kg kg −1 organic matter); and c peat P : 0.0010, 0.0010, 0.0008, 0.0006, 0.0005)(kg kg −1 organic matter). These parameters were used to compute the release of N and P from peat (r peat N,P , kg ha −1 day −1 ): where c om C is content of C in organic matter (0.55 kg kg −1 ), imm peat N,P is fraction of the released N and P that is immobilized into microbial biomass (kg kg −1 ). imm peat N,P was calibrated against measured N and P concentrations in the runoff water (see Section 2.3. Calibration).

Appendix B.2.4. Nutrient Balance
Nutrient balance module takes the nutrient demand and the nutrient release as input, adds the atmospheric deposition to the nutrient supply, keeps track on the nutrient storage in the rooting zone and provides the nutrient outfluxes to the groundwater and surface runoff. These are inputs to nutrient export module. Rooting zone depth (s depth , m) is parameter that is inherited from the SpaFHy [30]. First, the mean N and P concentration in the rainwater is multiplied with the monthly infiltration calculated with SpaFHy to obtain the nutrient input with atmospheric deposition (depo N,P , kg ha −1 ). N and P storage in the rooting zone was obtained as: sto N,P (t) = sto N,P (t − 1) + r N,P + depo N,P − U N,P m , where sto N,P is storage of N and P in the rooting zone (kg ha −1 ) at time-steps t and t − 1, r N,P is N and P release from organic matter decomposition (r min N,P or r peat N,P depending on the site type), and U N,P m is the monthly nutrient uptake by the stand and ground vegetation. Only positive values for sto N,P are allowed. Thereafter, N and P concentration (c water N,P , kg m −3 ) in the rooting zone was calculated as: c water N,P = sto N,P 10 4 * W m * s depth , where W m is the monthly mean water content (m 3 m −3 ), s depth is the rooting layer depth (m). Thereafter, N and P flux to groundwater (Q dr N,P , kg ha −1 month −1 ) was calculated as: Q dr N,P = 10 4 * Q dr * c water N,P , where Q dr N,P , kg ha −1 month −1 , Q dr is water flux down from root layer (m month −1 ). The nutrient flux from soil to surface runoff (Q ex N,P , kg ha −1 month −1 ) was obtained as: Q ex N,P = 10 4 * Q ex * c water N,P , where Q ex is water flux from soil to surface runoff (m month −1 ). Then, Q dr N,P and Q ex N,P were subtracted from sto N,P .

Appendix B.3. Nutrient Export
The nutrient export module moves the N and P in groundwater (Q dr N,P ) and in surface runoff (Q ex N,P ) to the receiving waterbody, keeps track of N and P storage and concentration in the groundwater, accounts for the nutrient retention and calculates the time delay connected to the groundwater transport. Hydrological processes close to soil surface have a shorter characteristic timescale than those deeper in soil [70]. The residence time of surface runoff was assumed to be shorter than the monthly time step applied in NutSpaFHy; therefore Q ex N,P was transferred to catchment outflux without a time delay and without retention.
In the module initialization, the euclidian distance (dist, m) and the elevation difference between each grid-point and the closest receiving water-body grid-point were calculated, and the mean slope (slope, m m −1 ) of the water flow path was solved accordingly. Using the a macro-scale hydraulic conductivity of K sat = 10 −4 m s −1 for both mineral soil slopes [53], and for peat soils [71], the time delay (t delay , months) from the generation of Q dr N,P to the visibility in Q runo f f N,P was obtained as: where conv to month is conversion of K sat from m s −1 to m month −1 , and poro is the soil porosity in the grid-point (m 3 m −3 ). In the initialization of the Nutrient export module, the three dimensional Q dr N,P matrix (dimensions: time,x,y) was shifted with t delay along the time axis for each grid-point to account for the effect of distance and slope to the residence time of groundwater. Thin till soils and peatlands are the dominant soil deposits if Finland [72]. Water movement in both these soils is dominated by the top layer with macroporosity and preferential flowpaths [53,71], and where the hydraulic conductivity can be orders of magnitude higher than in the deeper layers. Thus, we assumed the lower boundary of the active groundwater storage (where the transportation and storage occurs) reach two times the depth of the rooting zone. Current groundwater storage of the catchment (GW sto , m 3 ) was determined as the water located between the WT and the lower boundary of the active groundwater storage. N and P storage in the catchment groundwater storage GW N,P (kg) was calculated by considering the outflow with the runoff, and the inflow with the delayed input matrix Q dr N,P delayed (kg). The inflow was further decreased by N and P retention according to distance to the receiving water body using empirical function [38]. The N and P outflow with the groundwater Q gwout N,P (kg) was calculated as: where A catchment is the catchment area in m 2 , and Q runo f f is the monthly runoff (m month −1 ). N and P surface runoff flux Q sruno f f N,P (kg month −1 ) to the receiving water body were calculated as a sum of grid-cell surface runoffs: Q sruno f f N,P = ∑ x,y Q ex N,P * conv to grid−cell , where x and y are the grid-cell coordinates, Q ex N,P is N,P surface runoff flux from each gridcell (kg ha −1 month −1 ) and conv to grid−cell converts the outflux to kg grid-cell −1 month −1 .
The new N,P storage in the groundwater was obtained as: GW N,P = GW N,P (t − 1) − Q gwout N,P + ∑ x,y Q dr N,P delayed * (1 − f ret N,P ) * conv to grid−cell , where GW N,P (t − 1) is the groundwater N,P storage in the previous time-step (kg), ∑ x,y Q dr N,P delayed is the sum over the catchment area of the delayed N,P input to groundwater (kg), and f ret N,P is N and P retention (kg kg −1 ) factor. The retention factor was empirically connected to distance to receiving water body (dist, m) following the data presented in [38]: f ret N = (15.4 * ln(dist) − 52.7) * 10 −2 and f ret P = (19.1 * ln(dist) − 61.9) * 10 −2 . Finally, the total N and P outflux (Q out N,P , kg month −1 ) from the catchment was obtained as a sum of N P fluxes in groundwater outflow and surface runoff: Q out N,P = Q gwout N,P + Q sruno f f N,P .
The concentration of N and P in the runoff was obtained by dividing Q out N,P with the runoff volume.