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

Abstract

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.

1. 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.

2. Materials and Methods

2.1. 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. 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³ ha⁻¹, $f_{conif}$ 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.

$\mathit{id}$	$\mathit{Area}$	${\mathit{T}}_{\mathit{sum}}$	P	V	${\mathit{f}}_{\mathit{conif}}$	$\mathit{Slope}$	${\mathit{D}}_{\mathit{water}}$	$\mathit{Peat}$
Calibration catchments
2 p	167	1118	674	156	0.88	1.94	70	0.42
10 p	74	1145	668	88	0.94	1.65	123	0.48
13 m	436	1454	661	148	0.83	5.50	105	0.05
14 m	154	1283	606	165	0.82	2.82	86	0.09
21 m	1053	1043	807	96	0.82	3.37	86	0.28
22 m	1560	942	600	61	0.89	5.82	218	0.25
24 m	1719	968	623	51	0.84	1.49	89	0.53
25 m	1072	1261	626	142	0.81	3.00	104	0.23
27 m	1373	942	600	25	0.89	3.91	235	0.04
31 m	31	1425	574	137	0.85	3.43	32	0.05
32 m	37	1439	583	145	0.90	2.45	51	0.08
33 m	51	1118	674	130	0.84	2.18	215	0.30
Test catchments
3 p	72	1106	731	166	0.89	4.36	226	0.17
6 p	49	878	697	74	0.88	3.96	351	0.37
9 p	75	898	731	62	0.85	3.02	403	0.31
15 m	1455	1287	667	108	0.88	1.53	126	0.34
16 m	505	1366	646	152	0.88	2.05	49	0.42
17 m	1966	1275	665	140	0.83	1.97	52	0.33

). 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³ ha⁻¹. The distribution of site fertility classes in peatland and upland sites are presented in

Table A1. Volume fraction of coniferous trees ($f_{conif}$), areal share of fens, bogs and open peatlands ($fen$, $bog$, $open$), and peatlands in different fertility classes ($p_{rich}$, $p_{medium}$, $p_{poor}$), and upland mineral soil sites in different fertility classes ($m_{rich}$, $m_{medium}$, $m_{poor}$) in the catchments.

$\mathit{id}$	${\mathit{f}}_{\mathit{conif}}$	$\mathit{fen}$	$\mathit{bog}$	$\mathit{open}$	${\mathit{p}}_{\mathit{rich}}$	${\mathit{p}}_{\mathit{medium}}$	${\mathit{p}}_{\mathit{poor}}$	${\mathit{m}}_{\mathit{rich}}$	${\mathit{m}}_{\mathit{medium}}$	${\mathit{m}}_{\mathit{poor}}$
Calibration
2	0.88	0.17	0.22	0.03	0.03	0.46	0.09	0.03	0.55	0.00
10	0.94	0.02	0.35	0.11	0.00	0.23	0.25	0.00	0.48	0.02
13	0.83	0.02	0.03	0.00	0.01	0.03	0.02	0.19	0.55	0.01
14	0.82	0.05	0.04	0.00	0.01	0.12	0.01	0.25	0.63	0.00
21	0.82	0.06	0.22	0.00	0.02	0.24	0.06	0.02	0.67	0.00
22	0.89	0.02	0.17	0.06	0.02	0.16	0.08	0.01	0.73	0.00
24	0.84	0.05	0.39	0.09	0.04	0.30	0.16	0.00	0.45	0.00
25	0.81	0.10	0.13	0.00	0.02	0.19	0.04	0.25	0.43	0.00
27	0.89	0.00	0.02	0.02	0.02	0.03	0.00	0.00	0.09	0.17
31	0.85	0.00	0.05	0.00	0.00	0.03	0.02	0.12	0.77	0.03
32	0.90	0.00	0.08	0.00	0.00	0.03	0.05	0.06	0.65	0.06
33	0.84	0.06	0.24	0.00	0.01	0.28	0.05	0.06	0.64	0.00
Test
3	0.89	0.03	0.14	0.00	0.02	0.15	0.02	0.06	0.76	0.01
6	0.88	0.01	0.33	0.04	0.01	0.27	0.09	0.00	0.62	0.00
9	0.85	0.01	0.15	0.15	0.00	0.22	0.08	0.00	0.69	0.00
15	0.88	0.01	0.32	0.02	0.00	0.17	0.19	0.02	0.48	0.11
16	0.88	0.03	0.39	0.00	0.02	0.33	0.09	0.05	0.51	0.00
17	0.83	0.13	0.19	0.00	0.03	0.31	0.02	0.19	0.45	0.01

.

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.

2.2. 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_{runoff}$). 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 ($sfc$) 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$), $sfc$ 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 $W{T}_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 ($Im{m}_{peatNP}$ and $Im{m}_{minerNP}$) 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_{minNP}$, $r_{peatNP}$) 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_{drm}$) and surface runoff ($Q_{exm}$), we calculate the respective nutrient fluxes to groundwater ($Q_{drNP}$) and as a form of surface runoff ($Q_{exNP}$).

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_{drNPdelayed}$, 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.

2.3. 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_{Nobs}$) and P concentrations ($c_{Pobs}$) 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 ($im{m}_{Npeat}$, $im{m}_{Nminer}$ and $im{m}_{Ppeat}$, $im{m}_{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)).

$$c_{obsm}=s\ast c_{predm}+{\epsilon }_m,$$

(1)

where $c_{obsm}$ was the observed concentration of N and P in month m, $c_{predm}$ was the predicted concentration of N and P for month m, and s is a slope parameter, and ${\epsilon }_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)². 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 ( $im{m}_{Npeat}$, $im{m}_{Nminer}$ and $im{m}_{Ppeat}$, $im{m}_{Pminer}$ ) can be predicted from catchment characteristics. Using the calibration catchments, we predicted the optimized immobilization parameter values and the catchment characteristics in linear step-wise 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. (Table 1 and Table A1, Appendix A).

2.4. 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.

2.5. 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. Clear-cut scenarios calculated for Catchment 2 (area 166 ha). In all scenarios, clear-cut was only conducted to grid where stand volume exceeded 100 m³ ha ⁻¹.

Scenario	Distance, m	Clear-Cut Area, ha	Harvested V, m3	Mean Harvested V, m3 ha −1
Uncut	-	-	-	-
>100 m	>100	41	9066	211
<35 m	<35	43	6955	168
Peat <100 m	<100	30	4977	168
Min <100 m	<100	27	4979	184
sfc 4,5,6	no limit	21	3032	144
sfc 1,2,3	no limit	15	3035	196
<10 m	<10	34	5582	166
<50 m	<50	60	10,298	174
<70 m	<70	73	12,775	174
<90 m	<90	82	14,445	175
<110 m	<110	91	16,118	177
<130 m	<130	99	17,802	179
<150 m	<150	105	19,154	181
<170 m	<170	112	20,471	183
<190 m	<190	118	21,692	184
Rand 1 ...10	no limit	80	15,000	188

). To follow common forest management practice, we applied clear-cut only at grid cells where the stand volume exceeded 100 m³ ha⁻¹. 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 ($sfc$ 4,5,6) vs. high fertility sites ($sfc$ 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.

3. Results

3.1. Model Calibration and Immobilization Parameters

The mean (and range) values for $im{m}_{Npeat}$, $im{m}_{Nmin}$, $im{m}_{Ppeat}$, $im{m}_{Nmin}$ 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. Nitrogen and phosphorus immobilization parameters as a function of catchment characteristics calculated using stepwise regression analysis. Catchment characteristics are shown in Table 1 and Table A1.

Variable	${\mathit{imm}}_{\mathit{N}\mathit{peat}}$	${\mathit{imm}}_{\mathit{N}\mathit{miner}}$	${\mathit{imm}}_{\mathit{P}\mathit{peat}}$	${\mathit{imm}}_{\mathit{P}\mathit{miner}}$
$Intercept$	0.652 ($p<0.001$)	0.894 ($p<0.001$)	0.846 ($p<0.001$)	0.882 ($p<0.001$)
$f_{conif}$	0.282 (p 0.013)	-	-	-
$bog$	−0.150 (p 0.009)	-	-	-
$m_{poor}$	-	0.284 (p 0.038)	-	-
$Adjusted{R}^2$	0.607	0.301	0.031	0.0
$RMSE$	0.019	0.019	0.070	0.054

). 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.

3.2. 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 (Figure 3 and Figure 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 (Figure 3 and Figure 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⁻¹ year⁻¹, and the P export from 0.1 to 0.25 kg ha⁻¹ year⁻¹ (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).

3.3. 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⁻¹ year⁻¹, respectively (Figure 6). Locating a similar size of clear-cut 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 ($sfc$ 1,2,3) led to clearly higher specific exports than in lower-fertility sites ($sfc$ 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.

4. Discussion

4.1. 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.

4.2. 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₂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 $im{m}_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.

4.3. 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 (Figure 3 and Figure 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.

4.4. 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 $im{m}_{Ppeat}$ and $im{m}_{Pmin}$ 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].

4.5. 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.