Application of SWAT to assess peatland forestry impacts on water quality

Soil material that contains at least 30% content by mass of dead organic matter (Joosten & Clarke 2002) or organic matter with less than 20–35% mineral content can be considered peat (Turetsky et al. 2014). Peat deposition occurs when peat formation is higher than peat decay and is generally slow under waterlogged conditions (Holden et al. 2004). Favourable wet conditions lead to the accumulation of organic matter in specific areas, resulting in the formation of peatlands. Of the 3% peatland cover on Earth (Greenup et al. 2000), Finland has more than 30% peatland in its territory. In northern Finland, the largely drained peatlands usually cover more than 20,000 ha of the area and are over 30 cm in depth (Turunen 2008).

Peatlands are hotspots for biodiversity and ecosystem services (Kadlec & Wallace 2008). They are substantial carbon reservoirs in boreal and subarctic regions, constituting one-third of the global soil carbon pool (Yu et al. 2010). Peatland degradation occurs worldwide because of the drainage of peatlands (Ramchunder et al. 2012). In boreal and temperate zones, approximately 15 million hectares of peatland have been drained for forestry since the 1950s. Drainage activity peaked during the 1970s. In the past two decades, almost half of the total peatland area in Finland has been affected by drainage (Paavilainen & Päivänen 1995).

Peatland drainage modifies the water-holding capacity of peat and the amount of water that leaves it (Landry & Rochefort 2012). Subsurface flow dominates the runoff generation mechanism in drained peatlands, mainly because of the lower water table conditions compared with surface flow-dominated undisturbed peatlands (Holden et al. 2006). Previously drained areas may alter the natural state of hydrological conditions and nutrient loading (Paavilainen & Päivänen 1995; Nieminen et al. 2017b; Finér et al. 2021). Different forest management activities in peatlands can also affect nutrient leaching (Kreutzweiser et al. 2008).

Small catchment research on peatlands (Marttila et al. 2018) clearly shows the short-term (Finér et al. 2021) and long-term (Nieminen et al. 2017a) effects of forestry on water quality. Although long-term national monitoring programs are active at river outlets, little information on peatland drainage in larger catchments is available. Several models have been developed, ranging from the conceptual level to describing the processes of water movement in catchments with different soil types, land use, and management (Singh 2018) through the simulation of hydrology and water quality. A few modelling studies have focussed on peatland-dominated catchments to analyse the impacts of forestry (Lepistö & Syri 2001). However, the effect of the drainage network as the primary driver of change has not been accounted for directly. Often, monitoring large catchments fails to pinpoint the influence of specific land use operations. In recent years, there has been an increase in organic nitrogen exports from peatland-dominated catchments (Lepistö et al. 2021) in the Baltic Sea region. These require attention but are rarely included in catchment-scale simulations. This has led to an unclear understanding of how drainage of peatlands affects water quality at the catchment level.

This study was motivated by a lack of knowledge and the need for better modelling approaches to assess spatiotemporal changes in nutrient and suspended solids (SS) status in a peatland-dominated catchment. The soil and water assessment tool (SWAT) has been successfully implemented under various conditions (Molina-Navarro et al. 2017), including in peatlands (Melaku et al. 2020). However, to our knowledge, it has not yet been used in drained peatlands, especially in areas with intensive peatland forestry.

In this study, we aimed to observe the major land uses in peatlands and their specific properties, using the SWAT model to compare nutrient and SS predictions against observations in peatland-dominated catchments. The Simojoki catchment and one of its sub-catchments, Hosionkoski, in northern Finland, were used to assess the SWAT model performance in a peatland-dominated catchment. The specific objectives were to (i) identify the specific parameters of the SWAT model that can represent the effect of drainage at the catchment scale; (ii) determine annual loads and annual mean concentrations of SS, organic phosphorus (Org-P), total phosphorus (TP), organic nitrogen (Org-N), and total nitrogen (TN); (iii) understand how nutrients and SS vary across different geographical units across the catchment; and (iv) investigate the uncertainty ranges for the estimates.

Granites and gneisses are the most common bedrock types and are usually covered with till. There was sand over the bedrock near the catchment outlet. Forests on mineral soils (39%) and peatland forests (53%) dominate (Lepistö et al. 2014), whereas agriculture covers only approximately 3% of the catchment area, mainly near the catchment outlet (Perkkiö et al. 1995). The catchment contains 41.5% mixed soil types, whereas peat covers only 32.87%. A combination of a thin peat layer (THP) and mixed soil type (MST) covered 13.61% of the catchment area. In contrast, other combinations of multiple soil layers cover a very small portion of the catchment (<1%). We found that 54.26% of the catchment area had a slope of less than 2.5%, whereas 25.64% had a slope of 2.5–5% (Figure 1).

SWAT is one of the most globally used hydrological models for simulating the influence of hydrology and land use on nutrient concentration and loading. Catchment-scale modelling in SWAT includes the dynamics of quantitative and qualitative information and has the potential to simulate drained peatlands. However, the implementation of drained peatlands in a semi-distributed SWAT model is challenging. In addition, model conceptualization, parameter specifications, and management activities are aspects of uncertainty that can arise when applying SWAT (Gupta et al. 2005), particularly for peatland and organic soils.

The SWAT uses Manning's equation to define the flow rate and velocity. Water is routed through the river network using the variable storage routing method or the Muskingum River routing method. To compute sediment yield, SWAT uses a modified universal soil loss equation (MUSLE; Williams & Berndt 1977). The sediment yield of a specific catchment depends on surface runoff, peak flow rate, soil erosion, slope length, and crop management factors. Multiple processes are followed in SWAT to estimate the number of different forms of nitrogen (N). NO₃-N in SWAT represents the product of the water volume and average concentration in runoff, lateral flow, and percolation. SWAT uses a loading function for Org-N that depends on the concentration of Org-N in the top soil layer, sediment yield, and enrichment ratio (Williams & Hann 1978). As phosphorus (P) is associated with the sediment phase, soluble P estimation depends on the labile P concentration and runoff volume (Knisel 1982).

The Simojoki catchment contains many ditches (drainage networks and individual drainage lines) in peatlands, particularly in the forest area (Figure 1(b)). For different land uses, the area drained with open ditches varied substantially throughout the catchment (Figure 1). To implement existing drained peatlands in the SWAT model, we:

• - compiled a ditch database identified from aerial images of the catchment based on a method developed in an earlier study (Bhattacharjee et al. 2021a);

• - collected a high-resolution digital elevation model (DEM) (1 m × 1 m) from the National Land Survey of Finland (NLSF [Paituli] 2019);

• - the drainage (ditch) database contains the ditch length in a shapefile. Thus, to integrate the ditch database into the DEM, we burned the ditch database shapefile using the delineation option in ArcSWAT (Version-2012.10_7.24; SWAT2012.exe Revision 681);

• - the burned DEM was masked with the existing catchment boundary obtained from NLSF, particularly when the land was flat;

• - created a flow accumulation raster (1,000 ha) to identify potential stream outlets;

• - streams and outlets were created based on the flow accumulation raster;

• - verified streams and outlets using the existing database from NLSF [Paituli] (2019); and

• - delineated the watershed based on the final outlet of Bothnian Bay, which resulted in 111 sub-catchments.

This process essentially sets up the necessary spatial and topographic data for the SWAT model, ensuring that the drainage networks are accurately represented. Integrating a ditch database into the DEM is a critical step for modelling the hydrology of peatlands in the catchment. The simulation period was also considered so that the model could represent the existing ditch condition in the catchment as closely as possible. The verification steps with existing databases add a quality control layer to the process. Therefore, the final streams were generated based on the hydrological network of the drainage lines of the burned DEM.

The hydrologic response unit (HRU) is the smallest spatial unit in SWAT and consists of a unique combination of land use, soil type, and the slope band of a sub-catchment based on user-defined thresholds. Detailed spatial and temporal land use analysis of the catchment was provided in an earlier study (Bhattacharjee et al. 2021b). Based on land use analysis from Landsat (20 m × 20 m) (Bhattacharjee et al. 2021b), six classes were considered from 2015 to represent the catchment (NLSF [Paituli] 2019; Figure 1(c)). Soil data (25 m × 25 m) were collected from the Geological Survey of Finland [GTK] (2023). Ten soil types dominated the catchment (Figure 1(d)) and based on earlier studies (Menberu et al. 2021), we matched them to soils already in the SWAT soil database to be integrated with the model. Five slope bands (Figure 1(e)) were generated using the DEM. These slope bands were the most representative elevation variations in the catchment.

We used threshold values of 15, 20, and 20% for the land use, soil type, and slope, respectively, to define the final HRUs. These threshold values indicate that any land use, soil type, or slope band representing a percentage below these threshold values in each sub-catchment was deleted to simplify the model, subsequently proportionally increasing the remaining land use, soil type, and slope band. As mentioned in an earlier study, commonly used thresholds might lead to substantial information loss in catchment landscapes of the catchment (Her et al. 2015). Thus, the thresholds were chosen carefully, resulting in 634 HRUs representing the drained peatland areas within each HRU and sub-catchment.

We used open-access gridded weather data (10 km × 10 km) provided by the Finnish Meteorological Institute (FMI) from the Paituli server (NLSF [Paituli] 2019) for all meteorological inputs (precipitation, temperature, solar radiation, and relative humidity) except wind speed. We used simulated wind speed data based on a modified exponential equation in SWAT (Neitsch et al. 2011). These meteorological stations cover the entire Simojoki catchment at 26 different locations. The long-term (1996–2015) average annual temperature in the area varied between −3.83 and 4.94 °C, and the long-term annual average precipitation ranged from 580 to 746 mm. The mean monthly temperature (1996–2015) varied between 20.5 °C (July) and −16.3 °C (January). Average monthly precipitation ranged from approximately 30 mm (January) to 83 mm (July).

We collected daily streamflow and concentration data of Org-N, TN, Org-P, TP, and SS from the Hertta database (Suomen ympäristökeskus [SYKE] 2020). These monitoring data were used as the observed time series in the SWAT model to account for temporal changes during calibration and validation in two different periods. Since 1996, the sampling frequency at both outlets (Figure 1(a)) of the Simojoki catchment has been 11–25 samples per year (Lepistö et al. 2014). We used two monitoring stations to calibrate the model (i) Hosionkoski sub-catchment in the central part of the catchment and (ii) Simojoki at the catchment outlet (average summary is shown in Table 1). To estimate the load from the measured streamflow and concentration data, we applied a regression method (Walker 1996) and followed the steps described in a previous study (Bhattacharjee et al. 2021b).

We used the Penman–Monteith equation to estimate potential evapotranspiration and the variable storage routing method to define flow routing. First, we focussed on calibrating the SWAT model for flow and then for SS, Org-P, TP, Org-N, and TN sequentially. The full SWAT calibration period was from 1985 to 2005, with the first 10 years as the warm-up period.

We chose the SUFI2 algorithm in the SWAT-CUP (V-5.1.6.2) to perform semi-automatic calibration (Abbaspour 2015). First, we analysed the sensitivities of these parameters. As mentioned in earlier studies (Arnold et al. 2012), there are some calibration parameters in SWAT-CUP for which the model result can vary to a certain extent. Sensitivity analysis is an effective method for understanding the effects of the parameters in the model. Initially, we identified multiple parameters for use in the sensitivity analysis. We used parameter values from published literature (Menberu et al. 2021) to set up the initial parameter ranges of the physical and hydrological properties of peatlands. We used a Latin hypercube (LH) sampling scheme to calculate the sensitivity index of the parameters for global sensitivity analysis. The most sensitive parameters, based on t-stat and p-value for the calibration of flow, nitrogen, phosphorus, and SS, are shown in Supplementary Appendix Figure A1 where t-stat is the coefficient of a parameter divided by its standard error. If a coefficient is large compared with its standard error, then it is different from 0 and the parameter is sensitive. The p-value tests the null hypothesis (the coefficient is equal to zero, with no effect). A low p-value (<0.05) indicates rejection of the null hypothesis which represents the parameter associated with the response variable or affected by it (Abbaspour 2015).

We found some sensitive parameters (i.e., LAT_ORGP) that should be used only for calibration when the observed data are available to justify the selected ranges. As we had data from earlier studies that focussed on parameters related to peatland forestry (Huttunen et al. 2016; Räty et al. 2020), we used these parameters for calibration.

Based on the propagation of the parameter uncertainties, two prominent factors were considered when calibrating the model response variables, as recommended in a previous study (Abbaspour 2015):

• At least 70% of the measured data (observations) were between the 2.5th and 97.5th prediction percentiles (95PPU) of the model results. The measurement factor by which the uncertainties are addressed is known as the p-factor.

• The standard deviation of the measured data was greater than the difference between the 2.5th and 97.5th percentile predictions. The gap is known as the thickness of the 95PPU band or the r-factor.

Next, the ranges of parameters that fulfilled the above criteria were considered to check the reasonable agreement between the modelled and measured outputs. The Kling–Gupta efficiency (KGE) (Gupta et al. 2005) was selected as the objective function in the SWAT-CUP as recent studies have reported that KGE is suitable for evaluating the variability, peak, and mean of observations (Pechlivanidis et al. 2011).

Next, we started flow calibration based on the initial values of the sensitive parameters. We used 1,000 simulations in the SWAT-CUP for the first iteration. Because the initial model output did not match the observations well, we ran another iteration with the same number of simulations. After verifying the results of the second iteration, we run another 500 simulations. In each iteration, the parameter ranges were altered to obtain the most plausible model output based on the objective function. Finally, 23 parameters with reasonable ranges that affect the simulated flow (Table 3) were selected.

After calibrating the flow parameters, we started with the SS calibration and then moved to phosphorus and nitrogen. In addition to the flow parameters, we found nine parameters for calibrating the SS, nine for phosphorus, and three for nitrogen. We used three iterations with 500 simulations at each calibration stage to calibrate the nutrients and SS.

Finally, a split-sample validation was performed based on calibrated parameter values (Supplementary Appendix Table A1) in an independent dataset. The applicability of the calibrated model was tested from 2006 to 2015.

In addition to KGE, the Nash–Sutcliffe efficiency (NSE) was also used to evaluate the model. Recent studies have argued that using NSE could lead to a poor representation of flows, whereas KGE does not have this issue (Gupta et al. 2005; Pechlivanidis et al. 2011). Thus, we used two metrics for model evaluation to add robustness to the calibrated model. NSE varied between −∞ and 1.0. A value of 1 represents a perfect fit between the simulated and observed values, whereas values less than zero indicate unacceptable performance (Moriasi et al. 2015). We also evaluated model performance using R² and PBIAS. The R² evaluates the dynamics between the observed and simulated values, whereas the PBIAS quantifies the difference between the observed and simulated values. Thus, a PBIAS value of zero indicates the best result. Based on the output of these statistical indices, we selected the best simulations to finalize our calibration for the studied catchment.

For the calibration period (1996–2005), the ranges of the fully calibrated parameters were considered to generate the uncertainty (95PPU) boundary, whereas, for the validation (2006–2015), we used a single parameter value for the model evaluation (Supplementary Appendix Table A1). We also added the results of another calibrated model without peatland drainage. We integrated the results of both models with and without peatland drainage only in Table 2 in the Results section because it is important to know how the modelling approaches represent the drained peatland.

We further developed a Python script to analyse the calibration and validation results (Supplementary Appendix C) and extract loading information for any reach, sub-catchment, or HRU from the SWAT output database. The output results of the script can be used to map information regarding the nutrients and SS at different geographical boundaries.

The performance metrics for flow, nutrients, and SS in the SWAT-CUP varied from very good to unsatisfactory (Moriasi et al. 2015) for the Hosionkoski sub-catchment and the entire Simojoki catchment. We also added indices from the SWAT-CUP model without peatland drainage (in parentheses, Table 2). For all indices, the model performance was better when peatland drainage was implemented within the model. For the flow and nitrogen, the performance metrics were satisfactory or higher for calibration and validation. However, the estimations for phosphorus and SS varied for different indices (Table 2). As the observed SS concentrations were low and mainly contained organic fractions, the SS budget was not closed, which affected the modelling accuracy of the SS and phosphorus.

Based on the input and output components of the SWAT model, GW flow (shallow aquifer) had the highest share of water balance in the Simojoki catchment, followed by ET and the other components (Table 3).

We also compared the simulated potential ET with the potential ET available from gridded data (1 km × 1 km) from the FMI (Pirinen et al. 2022). From 1996 to 2005 (April to September of each year), we compared the average monthly potential ET (Supplementary Appendix Figure E1) and found that the difference was always less than 10%.

During the calibration period (1996–2005), based on the model results with peatland drainage, the uncertainty boundary (95PPU) captured 74% of the observed flow at Simojoki station, whereas the percentage was 70% for the upstream Hosionkoski catchment. The percentage was almost the same (69%) during the validation period in the Simojoki catchment area. The baseflow was separated from the simulated flow. After separation, we found that the simulated baseflow ratio was within 20% of the observed baseflow ratio. This can be considered as a benchmark for satisfactory calibration, as suggested in an earlier study (Bracmort et al. 2006). Gridded meteorological inputs were used to achieve better spatial coverage. Thus, these data minimized the uncertainty issues that could affect the simulated flow. ArcSWAT considers a meteorological station near the centroid of each sub-catchment. For low-flow conditions, the modelling results followed the pattern of the observations, but not for high-flow conditions caused by snowmelt (Figure 2).

For Org-N and TN, the uncertainty band estimated by SWAT-CUP included more than 82% of the observations in both catchments for calibration. The calibrated Org-N differed from the observations by 0.6% in the Hosionkoski sub-catchment and 11.2% in the Simojoki catchment. For validation, the differences were 9.4 and 16.1% for Hosionkoski and Simojoki, respectively. The difference between the simulated and observed TN values varied from 20 to 3.2% from upstream to downstream (Figure 3).

The ratio between the average Org-N concentration of the observations and simulations was 1.04 for the Simojoki catchment, whereas it was 0.9 for the Hosionkoski catchment (Supplementary Appendix Figure B1). The average TN concentration ratio between the simulations and observations for both catchments is 1.1.

The difference in TP load between the simulations and observations in the Hosionkoski catchment was approximately 20%, whereas it was higher (38%) in the Simojoki catchment (Figure 4). We found that 60% of the observations fell within the 95% prediction uncertainty (95PPU) boundary generated from the calibration results.

The simulated average TP concentration was 21.8 μg/l in the Simojoki catchment, which was 4.6% higher than the observations. For the Hosionkoski catchment, the simulated TP values differed by 25% from observations (Supplementary Appendix Figure B2). The simulation and observation analysis showed a higher concentration in Hosionkoski than in the Simojoki catchment.

For SS, the simulated model captured only 37 and 38% of the observations from 1996 to 2005 within the uncertainty boundary (95PPU) in the Hosionkoski sub-catchment and the Simojoki catchment, respectively (Figure 4). Although the prediction of the average SS load differed from the observations by 1.1% in the Simojoki catchment, the difference was high (∼50%) in the Hosionkoski catchment, as imbalances might also exist in the nutrient pool of the SWAT at the sub-catchment level.

The average simulated concentration of SS was 7.6 mg/l, whereas it was 3.8 mg/l for the observed measurements in the Simojoki catchment (see Supplementary Appendix Figure B2). The simulated concentration in the Hosionkoski sub-catchment was 1.5 times higher than the observed concentrations. Although the simulated suspended sediment (SS) concentrations exhibited significant variations compared with the observations in both catchments, the ratio of the average SS concentration between the observed and simulated data remained nearly constant.

The Org-N load was 80% higher in 2015 than that in 1996 at the Simojoki outlet, and this increase was even more pronounced for the TN load. However, the load ratio between the most downstream and upstream reaches was almost the same in 1996 and 2015 for both Org-N and TN. We found the highest load of 748 ton/year Org-N and 907 ton/year TN at the Simojoki outlet in the catchment.

Similar to Org-N and TN, Org-P and TP were higher in the middle and downstream parts of the Simojoki catchment. Between the most upstream and downstream reaches, the ratio of Org-P and TP loads was higher in 2015 than that in 1996. A total of 159 ton/year of Org-P and 169 ton/year of TP, the highest total load in a year, were found at the outlet of the Simojoki catchment.

The SS load at the furthest downstream reach in the Simojoki catchment was four times the SS load found at the furthest upstream reach in 2015 and almost five times the simulated load in 1996. In contrast to the nutrients, the highest total average SS load was found at a reach (reach 98) which was located above the most downstream outlet of the Simojoki catchment (Supplementary Appendix Figure F1).

The loading of nutrients in the peat soil in the clear-cut areas was the highest compared with all other land use categories. However, forests in the peat soil cover a large portion (31%) of the catchment. The loading of the mineral soil was also higher in the clear-cut areas than in the forest. From 1996 to 2015, the differences in loading between forests on peat soil and forests on mineral soil were 0.002 ton/ha/year for Org-N and TN, 0.0002 ton/ha/year for Org-P and TP, and 0.022 ton/ha/year for SS. Although peat extraction had the lowest area coverage within the Simojoki catchment, we observed the highest loading (1.03 ton/ha/year) of SS in that land use category.

The SWAT model has been used in various catchments worldwide and under northern conditions (Grizzetti et al. 2003). However, previous studies have not evaluated the SWAT performance in peatland-dominated catchments with intensive drainage activities, especially for various land uses in organic soils (forestry). The approach of burning the drainage network (collected from Bhattacharjee et al. (2021a)) into the DEM in a forestry catchment and integrating the process within the SWAT model resulted in good to satisfactory model performance for flow and nitrogen and slightly lower performance for phosphorus and SS. Specifically,

• I. The sensitive parameters identified in this study can be used to simulate the nutrients and SS in similar catchments. The gap that remains in this study is answering the question of why certain ranges were selected for specific parameters based on the sensitivity analysis. We relied on previous studies (e.g., Huttunen et al. 2016; Räty et al. 2020; Menberu et al. 2021) to select parameter ranges related to peat properties. Thus, the uncertainties in the current SWAT model can be reduced if there is an opportunity to obtain physical measurements of these parameters to support selected ranges.

• II. The nutrient and SS movement processes in multi-pools in SWAT depend heavily on the chosen parameters and their ranges (Supplementary Appendix Table A1), especially for forest land and organic soils. This process is more complex when the catchment has a drained peatland area of more than 30% (Rankinen et al. 2006). For example, owing to the organic fraction and the low concentration of SS and phosphorus, there are mismatches between the modelled output and the estimated loads from the monitored observations. Many parameters were dropped during the calibration process because of their insensitivity to peatland forestry. The uncertainty (95PPU) in estimates captured 35–81% of the observations that involved all plausible sources of uncertainties related to the model structure and the trade-off between parameter values and the objective function (Gupta et al. 2005).

• III. The current SWAT model showed a PBIAS of 9.6% between the simulated flow and the observations. TIMP (snowpack temperature lag factor) was the most sensitive parameter for calibrating flow. The amount of snowmelt is affected by the lagging factor because it accounts for the snowpack density, snowpack depth, and other exposures (Neitsch et al. 2011). As the range of TIMP was between 0.64 and 0.84 in the current SWAT model, we found that the mean air temperature for each day during the studied period had a larger impact on the snowpack temperature on that specific day than the snowpack temperature for the preceding day.

• IV. The simulation of SS loading is challenging as transported SS in peatland-dominated catchments typically consists mostly of organic particles (Marttila & Kløve 2015) with very fine particle sizes (Marttila et al. 2016). These physical processes, including the slow settling velocity and wash load nature, cannot be fully simulated using the SWAT. Compared with the observations (9,128 ton/year), the SWAT model generated approximately 11,292 ton/year SS load at the Simojoki outlet, whereas the simulated load was 1.4 times the observed load at the Hosionkoski outlet.

• V. Within the soil profile, the soil surface contained the highest amount of nutrients and organic matter. As organic soils have a distinct process from mineral soils (Menberu et al. 2021), the movement of nitrogen from the fresh organic pool to the inorganic mineral pool in the HRUs was insignificant compared with the concentration of Org-N in the baseflow (parameters of nitrogen calibration in Supplementary Appendix Table A1). The differences between the simulated load and the observations for Org-N and TN varied by 11.2 and 3.2%, respectively. In a previous study (Huttunen et al. 2016), the authors found a 1,012 ton/year gross TN load between 1991 and 2011, whereas the current SWAT model resulted in a 907 ton/year TN load. A minimal change in the Org-N concentration in the baseflow substantially affected model calibration. Saturated hydraulic conductivity and bulk density in the peat soil also affected nitrogen processes in the catchment, as nitrogen was highly reactive in the peat soil.

Some inorganic processes might have occurred in the mineral soil layer of the catchment, but during 1996–2015, denitrification was the dominant process. When the temperature was 5 °C or higher on a specific day, nitrification and volatilization could also occur in the catchment (Neitsch et al. 2011). Most plant-essential nutrients are cations that are attracted to and sorbed onto negatively charged soil particles. As plants extract these cations from the soil solution and the peat soil particles are positively charged, the leached nitrogen in the current SWAT model was not noticeable compared with nitrogen in the surface runoff. Although the Simojoki catchment has lakes covering 5.7% of its total area, nutrient retention within the catchment is not sensitive to lake effects. The parameters responsible for calibrating Org-N also fit well with TN calibration. However, peat properties are site specific and may vary owing to the different peat decomposition levels in different peatlands.

In response to a concentration gradient, phosphorus usually moves through the migration of ions over small distances on the soil surface in SWAT (Neitsch et al. 2011). Solubility is limited to phosphorus in most environments compared with nitrogen, which is highly mobile (Chaubey et al. 2006). Peat soils contained more than 50% phosphorus, whereas mineral soils contained approximately 40%. Runoff in the drained peat soils was higher than that in the undisturbed soils, which also affected the phosphorus concentration. Similarly, the concentration of Org-P in the baseflow also influenced phosphorus calibration, as it was the most sensitive parameter. Because peat soils have high water-holding capacity, they usually have little surface runoff. After the snow melting period, TP was higher than in the other seasons, during which stagnant water in peat soils might affect phosphorus loading (Huttunen et al. 2016).

Previous studies have already mentioned (Joensuu et al. 2002; Nieminen et al. 2018) that drainage activities and ditch maintenance are the main reasons for the increased nutrient and SS loads and concentrations in peatland-dominated catchments. In a previous study (Lepistö et al. 2008), the authors noticed an increasing trend in Org-N and found that changes in hydrological dynamics and soil temperature during low flows might affect its concentration. Temporal estimates using the current SWAT model for the Simojoki catchment showed an increasing nitrogen trend.

Clear-cutting in peat soils has not been implemented on a large scale in the Simojoki catchment area. However, the loading is predicted to increase in the future (Lepistö et al. 2008). From the spatial distribution of loading in different land uses (Figure 7), the clear-cut peat soil contained the maximum loading of all nutrients. In a previous study based on an export coefficient database, the authors also found the maximum loading for the same land use (Launiainen et al. 2014; Bhattacharjee et al. 2021b). However, the loading estimates from the export coefficient differ from those of the current SWAT model. For example, the authors mentioned 34 kg/ha/year TN, 3.9 kg/ha/year TP, and 2,400 kg/ha/year SS, whereas we noticed 50 kg/ha/year TN, 6 kg/ha/year TP, and 1,200 kg/ha/year SS from the forest in the peat soil. A plausible explanation for the difference in loads could be the nutrient decay processes in large rivers, which were not covered by the export coefficients. The export values were monitored in different empirical studies that included seasonal loading information from different land use activities but did not include decay processes from land to sea in the surface water system. In addition, the SWAT model is needed in this current study to provide a more detailed, process-based understanding of nutrient export and its spatial and temporal variability, which the export coefficient model cannot fully capture. While both methods offer valuable insights, the SWAT model's ability to simulate hydrological processes, nutrient decay and transformation, and the effects of management practices under changing conditions make it essential for addressing the research questions in this study.

Historical land use activities may show impacts for years or even decades. For example, the long-term prediction of Org-N in this study resulted in increasing trends from drained peatlands, which was in line with another recent study (Lepistö et al. 2021). Therefore, from a land manager's perspective, it is beneficial to know the spatial and temporal distribution of the drainage network, nutrients, and SS loading. Especially in the catchments containing a high percentage of peat soil, a proper management plan is essential.

Therefore, the modelling study can be used to understand the comprehensive status of a catchment, starting from knowing peatland drainage conditions to considering future projections. The study can also be used to monitor or achieve the good ecological status of surface water bodies based on the guidelines of the Water Framework Directive (WFD 2022–2027) (Suomen ympäristökeskus [SYKE] 2021).

For instance, the catchment manager can use the modelling approach to focus on the catchment-scale consequences of land use changes for nutrient loads and SS exports. However, including more detailed spatial information about drainage activities and the effects of hydrometeorological inputs would make the proposed approach more practical. Thus, the manager can apply the peatland drainage activity in the SWAT model in future hydrological modelling efforts by considering the impact of climate change.

The application of the SWAT tool in a drained peatland catchment demonstrated satisfactory calibration for flow and nitrogen, although there were differences in modelling SS and phosphorus. Despite these uncertainties, the selected parameters and ranges provide a good foundation for simulating physical processes in peatland-dominated catchments.

The observed increasing trend in Org-N during long-term predictions of drained peatlands is particularly important. Despite showing a PBIAS of 9.6% for simulated flow, the current SWAT model exhibited variations in loading between simulations and observations for SS, Org-N, TN, Org-P, and TP. Spatially, the clear-cut area's peat soil displayed the highest loading compared with other land uses, with the downstream outlet of the Simojoki catchment experiencing the maximum annual loading of nutrients, excluding SS and phosphorus.

While forestry practices are inherently dependent on local interpretations, our approach to integrating peatland drainage into SWAT holds promise for future hydrological modelling endeavours. The spatial estimation of nutrient distribution and SS loading has substantial environmental management implications, offering valuable insights for formulating catchment management strategies and implementing water protection measures.

This work was part of the Nordic Centre of Excellence BIOWATER, funded by Nordforsk under project number 82263.

J.B. and H.M. conceived the idea and developed the methodology. J.B. and H.M. performed formal analysis and wrote the first draft of the manuscript. E.M.N. reviewed and edited the manuscript. B.K. supervised and edited the manuscript.

The authors declare that they have the consent to participate.

The authors declare that they have the consent to publish.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.