AGRICULTURAL AND FOOD SCIENCE Agricultural and Food Science (2023) 32: 112–127 112 https://doi.org/10.23986/afsci.125385 National-scale nitrogen loading from the Finnish agricultural fields has decreased since the 1990s Inese Huttunen1, Markus Huttunen1, Tapio Salo2, Pasi Mattila2, Liisa Maanavilja2,3 and Tarja Silfver4 1 Finnish Environment Institute, Latokartanonkaari 11, 00790 Helsinki, Finland 2 Natural Resources Institute Finland, Latokartanonkaari 9, 00790 Helsinki, Finland 3 Geological Survey of Finland, P.O. Box 96, 02150 Espoo, Finland 4 Natural Resources Institute Finland, Halolantie 31 A, 71750 Maaninka, Finland e-mail address: Inese.Huttunen@syke.fi The national scale nutrient load modelling system VEMALA-ICECREAM was used to simulate agricultural total nitrogen (TN) loading and its trends for all Finnish watersheds for the period from 1990–2019. Across Finland, agricultural TN loading (ATNL) has decreased from 17.4 kg ha-1 a-1 to 14.4 kg ha-1 a-1 (moving 10-year averages) since the 1990s. The main driver of the decrease in simulated ATNL is a reduction in mineral fertilizer use, which has decreased the N surplus in the soils. The TN leached fraction, however, did not show a trend but did have high annual variability due to variations in runoff; this corresponds to an average of 14.4% of the TN applied. The ATNL was considerably higher in the Archipelago Sea catchment compared to other Finnish Baltic Sea sub-catchments, with the lowest ATNL found in the Vuoksi catchment in Eastern Finland. The highest decrease of ATNL was simulated for Vuoksi and Gulf of Finland catchments. In the Bothnian Sea, Bothnian Bay and Archipelago Sea catchments, the decreasing trend of ATNL was smaller but still significant, with the exception of the Quark catchment, where there was no significant change. The differences in decreasing trends between regions can be explained by the heterogeneity of catchment characteristics, hydrology and agricultural practices in different regions. Key words: agriculture, nutrient load, modelling, VEMALA model, leached fraction Introduction Intensive agricultural production and the use of mineral and organic fertilizers in agriculture has led to increased nitrogen (N) loading to inland and coastal water bodies, both in Finland and globally. The use of N fertilizers started to increase during the 1960s and reached its maximum in the 1990s (Luke 2023a, Luke 2023b). The proportion of mineral fertilizer N from total N input was close to 60% in the 1990s, with the proportion of manure N was at 35%. In recent years, the proportions are 55% and 40%, respectively. Elevated N loading is particularly detrimental in the Baltic Sea, which suffers from eutrophication due to its physical properties and anthropogenic drivers (HELCOM 2018). There is a consensus across Baltic Sea countries about the need to reduce nutrient loading in order to retain the sea’s positive environmental state (Baltic Sea Action Plan, BSAP 2007). The Water Framework Directive (WFD) is a policy instrument used to regulate water protection at the river basin scale (Rekolainen et al. 2003). The Nitrate Decree (Government Decree 2014) sets regulations for fertilizer and manure application methods, timing and amount of available N for different crops, as well as setting an annual maximum rate of 170 kg ha-1 total N in manure. For the purposes of BSAP and WFD implementation, there is a need to quantify the sources of apportionment of N loading to water bodies and for the evaluation of the trends of N loading over recent decades. As agricultural loading is diffuse, it is difficult to measure on site. Therefore, models are used to assess nutrient loading from the agricultural sector. There is a range of models available with varying levels of complexity that can be used for this purpose, from export coefficient models to dynamic process-based biogeochemical models. The scales at which the models are applied also vary to a great extent, usually in reverse correlation with their complexity. The most commonly used field scale ATNL models in Nordic countries are ICECREAM (Tattari et al. 2001), COUP (COUP manual 2022) and DAISY (Abrahamsen and Hansen 2000). Field scale models are often incorporated into the catchment scale models to assess agricultural loading at the catchment scale, e.g., VEMALA-ICECREAM (Huttunen et al. 2016) and DAISY-MIKE-SHE (Trolle et al 2019). Catchment scale models also integrate sub-models to simulate other non- agricultural sources of nutrient loading. ATNL estimates are also needed for the national GHG inventory. Emissions reported in the agricultural sector involve accounting for indirect nitrous oxide (N2O) emissions from managed soils which are directly related to N leaching from agricultural lands. Currently, the Finnish inventory uses the IPCC default value 0.3 kg N leached / kg N input, which is much higher than the estimates used in other Nordic countries, such as those used in Sweden (Swedish Environmental Protection Agency 2021) and Norway (Bechmann et al. 2012). Received 13 December 2022 / Accepted 16 August 2023 The Scientific Agricultural Society of Finland ©This is an open access article under the CC BY 4.012 I. Huttunen et al. 113 The aim of this study is to apply the VEMALA-ICECREAM modelling system for ATNL simulation to: 1) improve the temporal coverage and spatial resolution of the national scale N simulation related inputs to the VEMALA model, 2) improve volatilization and denitrification process description in the ICECREAM N model, and 3) provide national scale N loading from agricultural fields and analyze trends for the period from 1990–2019. Materials and methods Model description VEMALA-ICECREAM modelling system description The national-scale modelling system VEMALA used for nutrient (phosphorus (P) and N) loading simulates runoff processes, nutrient processes, leaching and transport on land, in rivers and in lakes in Finland (Huttunen et al. 2016, Korppoo et al. 2017). The VEMALA model provides an estimate of the input, output and retention of nutrients in all of the circa 40 000 lakes in Finland larger than 1 ha, as well as giving the source apportionment of loading: agriculture, forests, scattered settlements and point sources. The ICECREAM model (Tattari et al. 2001) is used to simulate P and N storages and fluxes in the soil (e.g., fertilisation, uptake, transport, annual balance, soil test P (STP) value changes) for all individual field plots in Finnish catchments and to calculate the estimates of mean catchment-scale P and N fluxes. The VEMALA model simulates the total nutrient (P and N) loading generated within the catchment, including diffuse and point sources, and retention in the streams and lakes, as well as simulating loading that enters the Baltic Sea from each watershed. A comparison of total P (TP) and total N (TN) concentrations and loads in the streams and lakes is one way to vali- date the model performance on a catchment scale. The scheme of the modelling chain for agricultural and total TN loading simulation using the VEMALA system is shown in Figure 1. The ICECREAM model is used to simulate the N fluxes for all, the approximately 970 000 field plots in Finnish river catchments, which cover an area of 23 359 km2. Loading from forested areas consists of natural background loading and forestry loading; in VEMALA, these are estimated using a combination of methods. Long term annual forest loading is based on Metsävesi project equations (Finér et al. 2020). The daily distribution of TN concentrations is simulated using a VEMALA-N model (Huttunen et al. 2016). The total TN concentrations and loads in the rivers are then calibrated against in-stream observations. In this study, ATNL results are summarized for Finnish catchments in the Baltic Sea sub-basins: Archipelago Sea (AS), Gulf of Finland (GF), Vuoksi in Eastern Finland (VUO), Bothnian Sea (BS), the Quark (QUA), and Bothnian Bay (BB) (see Figure S4). The Vuoksi catchment is also included in the study, since it covers a considerable part of Finland and, through the Neva River, contributes to the Baltic Sea. In this study, we report the gross ATNL generated on the national and catchment scale for the Baltic Sea sub-basins. ATNL retention in the inland water bodies must be taken into account when estimating the net ATNL entering the Baltic Sea. Fig. 1. The modelling chain for agricultural and total N loading simulation using the VEMALA modelling system Agricultural and Food Science (2023) 32: 112–127 114 Agricultural N loading simulation with ICECREAM The field-scale nutrient (P and N) loading model ICECREAM includes a process-based description for hydrology and N and P cycles in the soil, which are essential for the estimation of the effect of a changing climate and varying weather on nutrient (P and N) processes and loading. The ICECREAM model is based on the CREAMS and GLEAMS models (Knisel 1993). It has been developed further in order to be applicable to winter conditions and soil types in Finland. In the ICECREAM model, water flow in the soil is divided into three pathways: surface run- off and macropore flow (in clay soils) are subtracted from precipitation / snow melt and the rest of the water infiltrates through the soil profile as a matrix flow. ICECREAM calculates the daily balance of organic matter, organic N, ammonium-N (NH4 +-N) and nitrate-N (NO3 - -N) pools in soil from a 1 m deep soil layer. The input fluxes are plant residues, organic and mineral fertiliser, at- mospheric deposition, fixation by plants (for N) and the decay of organic matter. The processes that reduce N and P in the soil are plant uptake, denitrification (for N), transport via runoff and leaching, and volatilisation from mineral fertilizer (for N). More detailed description of hydrology simulation in the ICECREAM model is given by Sundholm (2021), and a description of N cycle simulation has been provided by Kämäri et al. (2019). Basic calibration of ICECREAM has been done in previous projects using field scale hydrological and TN leaching measurements. The processes of plant N uptake, denitrification and volatilization are described here in more detail. The plant N uptake is the highest outward flux in the N balance of the soil; it is important to simulate it as accurately as possible. The increase in plant biomass is related to air temperature sums over the vegetative season (Rekolainen and Posch 1993). The plant is divided into below- and above-ground biomass and yield, which all con- tain N in a user defined C:N ratio. Plant N uptake is simulated using a daily N demand, taking into account daily plant biomass increase, soil moisture stress, and N and P availability as limiting factor. (1) Where Bmat is biomass at maturity, input data, kg m-2, w is the plant dependent growth parameter, 2.0 for cereals, 1.0 for grass, and 3.0 for root crops. Denitrification is simulated according to equations (2), (3) and (4), and also depends on the NO3 --N content in the soil, according to the Michaelis-Menten formula (4). The denitrification potential, KdenN,max (g C m-2 d-1), is calculated using the mineralization rate of organic matter (Chatskikh et al. 2005): (2) where CFOM, CSMB and CSOM are the carbon pools of fresh organic matter, soil microbial biomass and native organic matter in g m-2, and λFOM, λSMB and λSOM are the efficient decay rates for the pools, respectively. fNO3 is the cali- brated parameter for adjusting denitrification. The denitrification rate, KdenN (g N m-2 d-1), is calculated as (3) (4) where KdenN,max (g N m-2 d-1) is the denitrification potential, calculated from the mineralization rate of organic matter; Fq is the dependency of denitrification rate on water-filled porosity; θ is volumetric water content; θs is the water content at saturation; Ft is the dependency rate on soil temperature; and Fcn is the dependency of the denitrifi- cation rate on NO3 --N concentration: (5) where NO3 -N is the soil NO3 --N concentration in mg kg-1. In the model, Fcn can obtain values between 0 and 1. 𝐵𝐵𝐵𝐵 = ( ∑𝑇𝑇𝑇𝑇𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 𝑇𝑇𝑇𝑇𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠,𝑠𝑠𝑠𝑠𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚 )𝑤𝑤𝑤𝑤 × 𝐵𝐵𝐵𝐵𝑠𝑠𝑠𝑠𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚 𝐾𝐾𝐾𝐾𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑,𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚 = 𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓3 ∙ (𝜆𝜆𝜆𝜆𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹 ∙ 𝐶𝐶𝐶𝐶𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹 + 𝜆𝜆𝜆𝜆𝑆𝑆𝑆𝑆𝐹𝐹𝐹𝐹𝑆𝑆𝑆𝑆 ∙ 𝐶𝐶𝐶𝐶𝑆𝑆𝑆𝑆𝐹𝐹𝐹𝐹𝑆𝑆𝑆𝑆 + 𝜆𝜆𝜆𝜆𝑆𝑆𝑆𝑆𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹 ∙ 𝐶𝐶𝐶𝐶𝑆𝑆𝑆𝑆𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹) 𝐾𝐾𝐾𝐾𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 = 𝐾𝐾𝐾𝐾𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑,𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚 ∙ 𝐹𝐹𝐹𝐹𝑡𝑡𝑡𝑡 ∙ 𝐹𝐹𝐹𝐹𝑞𝑞𝑞𝑞 ∙ 𝐹𝐹𝐹𝐹𝑐𝑐𝑐𝑐𝑑𝑑𝑑𝑑 ∙ 𝐹𝐹𝐹𝐹𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑑𝑑𝑑𝑑 𝐹𝐹𝐹𝐹𝑞𝑞𝑞𝑞 = 0.0116 + 1.36 1 + 𝑒𝑒𝑒𝑒−(𝜃𝜃𝜃𝜃 𝜃𝜃𝜃𝜃𝑠𝑠𝑠𝑠−0.815)/0.0896⁄ 𝐹𝐹𝐹𝐹𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐 = 1.17 × 𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁3−𝑁𝑁𝑁𝑁 32.7 + 𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁3−𝑁𝑁𝑁𝑁 I. Huttunen et al. 115 In earlier studies, simulated TN concentrations were overestimated throughout the 1990s because fertilization levels were higher than in the 2000s. Two parameters limiting the rate of vmax and concentration of Km in the Michaelis -Menten formula (5) were changed to achieve better responses between denitrification and NO3-N concentration in the soil. Changing the parameters (from 1.17 to 1.80 for vmax and from 32.7 to 10 for Km) did not provide any major increase in the simulated denitrification at the catchment scale. The additional linear denitrification coeffi- cient Flin in 1990s has been introduced in order to increase simulated denitrification and to give differing conditions for before and after EU regulations began in 1995. Linear coefficient had a value of 1.5 in 1990, and 1.0 in 1999. The reasoning for the new coefficient could be due to changes in the use of mineral fertilizers and manure N. In the period 1990–1994, before Finland’s entry into the European Union and it’s Common Agricultural Policy in 1995, mineral N fertilizer rates were aimed at providing the highest possible yields due to high yield values. The lower intensity of N fertilization during the latter half of 1990’s is indicated by a lower application rate of fertilizer N per hectare (Luke 2023b). Furthermore, the timing and rate of manure application was less regulated, which may also have led to higher N losses before 1995. However, it remains unclear how the extra fertilizer N has been distributed to various pathways of N losses from the soil and water system. The most likely loss pathways were denitrification as N2 and ammonia volatilization. During this project, an additional process is added to the model – volatilization from manure N. A suitable equa- tion has been found in the original GLEAMS manual (Knisel 1993). (6) (7) VOLNi, is ammonia volatilization, kg ha-1, AWNHi is ammonia in animal manure, kg ha-1, Kv,I is the volatilization rate constant, and ATPi is the mean daily air temperature. It is assumed that volatilization happens only from surface applied slurry. Volatilization continues over a 7-day period after application. For the implementation of the volatilization simulation, the percentage of the surface spread or injected slurry has also been added to the catchment scale ICECREAM model application. The volatilization rates suggested by Grönroos et al. (2009) vary from 30 to 60% of ammoniacal manure N. In the model, it is assumed that 59% of N in manure is in the form of NH4 +-N (Grönroos et al. (2009) give a range from 32% to 72% for cattle). The simulated volatilization for one test field was 67% of the surface applied NH4 +-N in manure (Fig.S1a). National scale input data for the model A comprehensive national-scale database has been collected; it is updated regularly and is maintained as input data for the VEMALA model (Huttunen et al. 2016). This database includes the following: daily air temperature and precipitation observations from the Finnish Meteorological Institute; water quality monitoring data in streams and lakes from SYKE for model validation; point loading data (including peat production loads) gathered from the Compliance Monitoring Data System (YLVA); and loading from scattered settlements based on a built environment information system and specific loading per person. The temporal coverage and spatial resolution of mineral fer- tilizer, manure, crop distribution and placement data has been improved within this project. Mineral fertilizer data Mineral fertilizer data estimates are based on the amount of mineral fertiliser sold at the level of Centres for Eco- nomic Development, Transport and the Environment (ELY Centres); the data is received from the Natural Resources Institute Finland (Luke). Luke calculates the amounts of P and N in sold fertilizers based on fertilizer companies’ annual declarations to the Finnish Food Authority and on information obtained from fertilizer companies directly. Updated mineral fertilizer data is implemented in the model for the period 1990–2019 (Fig. 2). The mineral N fertilizer usage has decreased by circa 30%, from 220 Gg a-1 to 150 Gg a-1, since the beginning of 1990s. The highest fertilizer usage is in Southwest Finland’s ELY region and is the lowest in the Ostrobothnia ELY region (Fig. 2a). The field-specific application rate of fertilizer N is estimated. The mineral fertilizer (and manure) usage has been distributed to each field depending on the cultivated crop and its P and N demand, according to the recommended 𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑖𝑖𝑖𝑖 = (𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝑉𝑉𝑉𝑉𝐴𝐴𝐴𝐴𝑖𝑖𝑖𝑖) × [1 − 𝑒𝑒𝑒𝑒−𝑘𝑘𝑘𝑘𝑣𝑣𝑣𝑣,𝑖𝑖𝑖𝑖] 𝑘𝑘𝑘𝑘𝑣𝑣𝑣𝑣,𝑖𝑖𝑖𝑖= 0.409 × (1.08)𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴,𝑖𝑖𝑖𝑖−20 Agricultural and Food Science (2023) 32: 112–127 116 fertilizer amounts for different crops required to receive agri-environmental support payments (according to MAVI 2009 guidelines). Then the field scale fertilizer inputs are further adjusted to match the annual areal ELY centre values of fertilizer sale. Manure data In the VEMALA model, manure data has to date been based on one year of data (2017). Animal number data for each farm and year for the period 1995–2019 is now available through the Finnish Food Authority, allowing the manure input data to be improved. The amount of N entering the soil from the manure of one animal of each animal type (cattle, swine, poultry, sheep and goat) was calculated in Luke for each year in the period from 1990– 2019; this was based on the N mass flow model version used in the national greenhouse gas inventory of Finland (Grönroos et al. 2009, Grönroos 2014, Statistics Finland 2021) and the national-scale shares of animal subcate- gories (dairy cows, heifers etc.) within each animal type (Luke 2021, Statistics Finland 2021). The information for the shares of each animal waste type (solid, slurry, urine, etc.) and the slurry incorporation method (injection or surface spreading) (Grönroos 2014) was included in the data. The share of slurry from total manure has increased from 50% at the beginning of the 1990s to 70% by the year 2012 (Grönroos et al. 2009, Grönroos 2014). The total amount of injected N in the slurry has also increased from about 5% at the beginning of the 1990s to 30% in the last decade (Grönroos et al. 2009, Grönroos 2014). The in- corporation method for the slurry is implemented in the ICECREAM for each field so that the amount of injected N matches LUKE estimates on the Finnish scale. 0 20 40 60 80 100 120 140 19 90 19 92 19 94 19 96 19 98 20 00 20 02 20 04 20 06 20 08 20 10 20 12 20 14 20 16 20 18N m in er al fe rt ili ze r, kg h a-1 a-1 Whole country Southwest Finland Ostrobothnia a Fig. 2. (a) Statistics for N mineral fertilizer usage in Finland 1990–2019 for the whole country, Southwest Finland and Ostrobothnia; (b) amount of mineral N fertilizer used in Finland – statistics and model inputs 0 50 100 150 200 250 19 90 19 92 19 94 19 96 19 98 20 00 20 02 20 04 20 06 20 08 20 10 20 12 20 14 20 16 20 18 N m in er al fe rt ili ze r, Gg a -1 Statistics Model b I. Huttunen et al. 117 The total amount of N in manure as input to the soil (Gg a-1) in Finland (Fig. 3a), which is used as input to the model, varies annually from 80 Gg a-1 to 93 Gg a-1. The reason for the increase of N in manure is the change in animal size and excretion. The increase is not caused by increased animal numbers as these numbers have actually decreased since the 1990s. The share of organic N in the total mineral and organic N fertilizer has increased from 31% in the 1990s to 39% in the 2010s. The estimation of N in manure returned to the soil is a challenging task, and there is a limited amount of national scale data available. VEMALA-ICECREAM estimates have been compared with N mass flow model estimate for the year 2014, which was 73 N Gg a-1 (Luostarinen et al. 2017) in manure returned to soil plus 11 N Gg a-1 left in pastures. The difference between the two methods is 10.5%, which is a reasonably close estimate. There are quite large differences between different areas of Finland in terms of N in manure usage (Fig. 3b). The areal distribution is based on the real animal type number on each farm (Finnish Food Authority data) and manure N input per animal in each animal type (Grönroos et al. 2009, Grönroos 2014). According to our simulations, the highest manure N usage was in the Pohjois-Savo area and the lowest was in the Uusimaa area. Detailed field-specific application of manure data is not known, and therefore is simulated using the VEMALA-ICECREAM model. The assumptions are that: manure is spread within a 20 km radius around the farm; manure is spread on fields with all crop covers (except fallow fields); and field-specific manure and mineral fertilizer usage depends on the cultivated crop and its P and N demand (according to MAVI 2009 guidelines). In areas of high livestock production these limits are often exceeded. Fig. 3. (a) The total amount of N in manure as an input to the soil (Gg a-1) in Finland, which is used as input to the model; and (b) N in manure used at the ELY centre level (kg ha-1 a-1) 30 40 50 60 70 80 90 19 90 19 92 19 94 19 96 19 98 20 00 20 02 20 04 20 06 20 08 20 10 20 12 20 14 20 16 20 18 N in m an ur e, G g a-1 a 100 0 10 20 30 40 50 60 70 19 90 19 92 19 94 19 96 19 98 20 00 20 02 20 04 20 06 20 08 20 10 20 12 20 14 20 16 20 18 N in m an ur e, k g ha a -1 Southwest Finland Ostrobothnia North Ostrobothnia North Savo Uusimaa Mean b Agricultural and Food Science (2023) 32: 112–127 118 Database of agricultural field characteristics The VEMALA-ICECREAM modelling system contains a detailed agricultural soil description for all field plots in Finland. Field plot register data, which includes field ID, field borders and cultivated crops, was provided by the Finnish Food Authority. Standard soil fertility analysis for the tillage layer was provided by Eurofins Viljavuuspalvelu Oy (or other companies) for around 40% of field plots. The analysis includes the sensory estimation of soil type for the field plots. For the rest, the soil texture was estimated based on data from the Finnish soil database (Lilja et al. 2017). The VEMALA-ICECREAM database contains the following characteristics for each field: farm number, field ID, area (ha), coordinates, 3rd level sub-catchment number, soil texture, STP level, class of soil organic matter con- tent, and slope. There are four soil types in the model: clay (main particle size <0.002 mm), silt (0.002–0.02 mm), sand (>0.02 mm), and organic soils (organic matter content in plough layer >20%). In the text, we refer to any soil with an organic matter content of >40% as peat soil. The proportion of the fields covered by mineral and organic soils are 2.08 × 106 ha of fields on mineral soils and 0.24 × 106 ha on organic soils. Crop distribution and location data was collected from the Finnish Food Authority for the period from 2000–2020. Data from 2000 is used for 1990–1999, because data for earlier years is not available. The field borders and total area of the fields do not change in the model, they are based on the 2012 data. During the period from 2000–2020, the grass crop area slightly increased, from 29% to 34% of the total area, and the cereal area has slightly decreased from 46% to 42% of the total area (Fig. 4). The grassland area increased from 2016–2019, with a possible reason being a shift towards organic production, with silage added to crop rotations, and the low economical profit in cereal production. This may also create a slight reduction in TN loading because mean specific loading from grass fields is usually lower than from cereal fields. Mann-Kendall trend tests, or, in the presence of covariates, partial Mann-Kendall trend tests (Libiseller and Grim- vall 2002), were used for testing the N input, ATNL and leached N fraction trends from 1990–2019 using R. The trends of TN river concentrations are estimated using seasonal Mann-Kendall tests. The moving average of previous 10-year values was used to analyse ATNL trends: (8) MA10 (kg ha-1a-1) – is the moving average of the previous 10 annual loading values Ln (kg ha-1a-1). Results Calibration of the denitrification parameter In this project, only one ICECREAM parameter, denitrification parameter fNO3, was used to calibrate deni- trification in both mineral and organic soils. Denitrification potential is multiplied by fNO3. The higher fNO3, 0 500 1000 1500 2000 20 00 20 02 20 04 20 06 20 08 20 10 20 12 20 14 20 16 20 18 20 20 Cr op a re a, 1 03 ha other crops green fallow grass spring wheat barley oats Fig. 4. Distribution of crops in Finland (ha) used in the VEMALA-ICECREAM system for 2000–2020 (based on Finnish Food Authority data) 𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀10 = 𝐿𝐿𝐿𝐿𝑛𝑛𝑛𝑛−10+1 + 𝐿𝐿𝐿𝐿𝑛𝑛𝑛𝑛−10+2 … . +𝐿𝐿𝐿𝐿𝑛𝑛𝑛𝑛 10 I. Huttunen et al. 119 the more denitrification and the less N load. In previous model versions, the value of fNO3 was 0.80, but this decreased when volatilization was added to the model. fNO3 values now vary from 0.48 to 0.64 and were manually calibrated against TN concentrations and loads in 9 Finnish rivers located in different regions of the country (river names are given in Table 1). Calibration of the VEMALA-ICECREAM modelling system is quite challenging because ICECREAM model field scale parameters should be calibrated against output of the other model (VEMALA), which simulates river concentrations. N sub-processes in agricultural soils The results for N sub-processes in agricultural soils for all the fields in two contrasting catchments are shown as a mean in Figure 5. The Aurajoki catchment (Fig. 5a) uses intensive agriculture, a high proportion of which is cereal cropping, with high mineral fertilizer rates applied and high manure use. The Vantaanjoki catchment (Fig. 5b) has lower mineral N input and much lower manure use. Despite the lower total N input to the soils in Vantaanjoki, there is only a small difference in the N in crop yields in these two catchments. N in yield slightly increased after the 2000s by about 10 kg ha-1 a-1 due to the increasing length of the growing season. However, dry summers have resulted in lower yields at the end of the 2010s. The N surplus shown in Table 2 is calculated as the difference between the total N fertilizer applied and the N in crop yield. A full N balance calculation is quite complex and is done in the ICECREAM model; it includes more fluxes such as mineralization, leaching, denitrification and volatilization. The N surplus in catchments similar to Vantaanjoki is lower (Table 2), meaning there is a greater decrease in ATNL. In catchments similar to Aurajoki, N surplus in soils has remained high (Table 2), meaning there is only a slight decrease in ATNL leaching. In other words, N fertilizer decrease leads to a decrease in N surplus in the soils, which in turn causes a decrease in ATNL. Denitrification from soils plays an important role in the N cycle in soils, both in nature and in the model. In the latter, denitrification is a correcting outward flux, where the excess N from the system is allocated. Simulated denitrification varies between different catchments, from a mean of around 40 kg ha-1 a-1 in Aurajoki catchment with higher N inputs to around 20 kg ha-1 a-1 in the catchments with lower N inputs. In the 1990s, there is elevated simulated denitrification due to the additional denitrification coefficient in the model. Table 1. Trends of observed daily TN concentrations, loads and simulated annual loads between 1990–2019 using seasonal Mann- Kendall tests River ELY-centre Agricultural loading, % Peat soils*, % Simulated monthly average TN concentration Observed monthly average TN concentration Comments tau p tau p Aurajoki Southwest Finland 74 11 0.002 0.985 0.002 0.974 High mineral fertilizer, manure usage Eurajoki Satakunta 68 16 -0.23** 0.001** -0.088 0.106 High manure usage, lake effect Lapuanjoki South Ostrobothnia 58 29 -0.355** 0** 0.035 0.566 High manure usage Porvoonjoki Uusimaa 53 5 -0.406** 0** -0.248** 0** Point load stops in 2004, only mineral fertilizer used Kalajoki North Ostrobothnia 52 35 -0.093 0.143 0.082 0.302 High manure usage Vantaanjoki Uusimaa 43 8 -0.31** 0 -0.233** 0** High point load share in 90-ties, only mineral fertilizer used Perhoonjoki Ostrobothnia 40 43 -0.095 0.212 0.07 0.369 Low % of agriculture Summanjoki Southeast Finland 35 18 -0.279** 0** -0.089 0.247 Low manure usage, low % of agriculture Siikajoki North Ostrobothnia 34 56 -0.111 0.131 0.209*** 0.013*** Low % of agriculture * = % of peat soils includes % of peat soils (with OM content >40%) from non-agricultural areas and % of organic soils (with OM content >20% in plough layer) from agricultural soils; ** = decreasing trend; *** = increasing trend Agricultural and Food Science (2023) 32: 112–127 120 Verification of simulated TN concentrations against river concentrations Agriculture is only one of the sources of TN loading and there is great variation between catchments. It is not possible to verify agricultural loading alone against river observations. However, in this project, we have verified total loading simulated against river observations, since agriculture contributes an important share of total load- ing. The hypothesis is that if total concentrations and loadings are correctly simulated, then agricultural loading should also be simulated well; this hypothesis may be valid in river catchments with a high ATNL share (above 50%). Table S2 shows simulated and observed mean TN concentrations, biases, Nash and Sutcliffe efficiency (NSE) criteria and source apportionment (as % of total loading) for the period from 2013–2020 for 9 monitoring sites located in different ELY-centers with different intensities of agriculture. The difference between mean simulated and ob- served TN concentration ranges from 0 to 16%, and the NSE criteria for daily simulated TN concentrations varies from 0.46 to 0.75, which is a reasonably good simulation result. The verification is done by comparing model bias and the NSE (Table S2) of mean TN concentrations for the 9 chosen river points. Figure 6 shows the simulated and observed daily TN concentrations for the Aurajoki catchment (NSE criteria was 0.59). Observed daily maximum concentrations or daily observed loads in the Aurajoki river are not decreasing, but annual loads present a slight decrease. The results for other rivers are shown in Figures S3, a-i. Table 1 summarizes the trends of observed and simulated monthly average TN concentrations. There is a significantly decreasing trend in observed TN concentra- tions in Porvoonjoki and Vantaanjoki, and an increasing trend in Siikajoki, while the simulated TN concentrations showed decreasing trends in five catchments, Porvoonjoki, Vantaanjoki, Summanjoki, Eurajoki and Lapuanjoki (Table 1). In Aurajoki, Kalajoki and Perhoonjoki, there is no trend in either simulated or observed TN concentrations. Fig. 5. Catchment scale simulated N sub-processes in agricultural soils: mineral N fertilizer, N in manure, N in yield, N loading, N denitrification, N volatilization a) for Aurajoki catchment and b) for Vantaanjoki catchment Table 2. Simulated N surplus kg ha-1 a-1 (Nfert + Nman – Nyield) in Aurajoki and Vantaanjoki catchments Period Aurajoki Vantaanjoki 1990–1999 51.9 23.8 2000–2009 41.0 5.4 2010–2019 43.8 9.3 I. Huttunen et al. 121 Table 2 summarizes also pecentage of peat soils in the catchments, however peat soils for agricultural areas inclu- de also organic soils in the peat soil percentage calculation. It could be concluded that in mineral soil dominated catchments with mostly mineral fertilizer use (Porvoonjoki, Vantaanjoki, Summanjoki), TN loading is decreasing, probably due to a reduction in the use of mineral fertilizer. However, there might be other factors decreasing the N leaching, such as an increase of grass areas replacing cereal areas. The Increase of no-till practices and the in- creased area of catch crops have also led to decreased N leaching. In peat soil and animal husbandry dominated catchments, simulated TN loading decrease is not so pronounced, probably due to multiple factors, such as mineral fertilizer having a smaller share of total N input; manure N input slightly increasing; and forest and forestry loading also playing an important role as agriculture has a lower load share in these catchments. In addition, the growth day degree temperature sum of growing season increased by about 180 °C on average in Finnish watersheds for the period from 1990–2019, which might cause higher organic N mineralization in peat soil dominated catchments. Spatial variation of agricultural TN loading The VEMALA-ICECREAM modelling system provides the possibility to simulate the ATNL for each field depending on its own characteristics and N input-output fluxes. The spatial variation of the ATNL ranges from 5 to over 30 kg ha-1 a-1 (Fig. 7). The highest ATNL is in Southwest Finland area due to this area having the highest mineral fertilizer use, high manure input and being an area of high cultivation of cereals and special crops (including cabbage, oth- er vegetables, sugar beets, etc.). The Satakunta and Ostrobothnia regions, and the northern part of North Savo have high ATNL due to their high combined mineral and manure N inputs. The western part of Southeast Finland has high ATNL due to cereals being the main crop in the area. Fig. 6. Simulated and observed daily concentrations for the Aurajoki catchment for 1990–2019 (NSE criteria is 0.59) Fig. 7. Spatial variation of simulated TN load from fields (kg ha-1 a-1, mean for 2010–2019) Agricultural and Food Science (2023) 32: 112–127 122 Agricultural TN loading results for 1990–2019 The ATNL is considerably higher in the AS catchment compared to other Baltic Sea sub-basin catchments, while the average lowest ATNL was in VUO catchment (Fig. 8). The ATNL for the whole of Finland has decreased since the 1990s, from 17.4 kg ha-1 a-1 to 14.4 kg ha-1 a-1 in the 2010s.. The highest decrease of ATNL, 32% and 25%, is sim- ulated for VUO and GF catchments, respectively. In the BB and BS catchments, the decrease of ATNL was smaller, 18% and 11%, respectively. In the AS and QUA catchments, the ATNL did not significantly change between 1990– 2019 according to the partial Mann-Kendal trend test (Table 1). Absolute values of the annual ATNL (Table S5) varies to a large extent depending on hydrological conditions; the maximum ATNL of 56.8 103 t a-1 was simulated during 1991 (a particularly wet year), and the minimum ATNL 17.7 103 t a-1 was during 2009. ATNL has decreased by 17%, from a mean value of 40.6 103 t a-1 in 1990s to 33.7 103 t a-1 between 2010–2019. Fig. 8. Agricultural TN loading (kg ha-1 a-1) and trendlines for the whole of Finland and catchments of the Baltic Sea sub-basins for the period from 1990–2019. The moving average is calculated according to equation 8. I. Huttunen et al. 123 The differences in decreasing trends in different regions can be explained through the heterogeneity of the catch- ment characteristics and agricultural practices. In mineral soil dominated VUO (the Vuoksi river has 19% of the peat soil of the whole catchment) and GF (the Kymijoki river has 17% of the peat soils of the whole catchment) catchments with relatively low manure N input, the ATNL decrease was highest. There is also a variation in the ATNL within the GF catchment depending on the main cultivation crop; this is because ATNL is higher in cereal cultivation areas. In the AS catchment, the decrease in ATNL is not so pronounced because of high manure N in- puts and cereal crop cultivation. The mild winters over the last few decades are producing high ATNL, especially in soils without winter vegetation cover. The highest simulated TN concentration in the Aurajoki river was during the mild winter of 2019/2020 (Fig. 6). In catchments with high manure N input and high shares of peat soils (BS, BB, ME), the decrease in ATNL is also not so pronounced. The leached N fraction of total N input, used in the Finnish inventory of GHG emissions for the indirect estimation of N2O emissions, varied from 8% to 21% between 1990–2019 (Fig. 9) and was on average 14.4%. The leached fraction did not change much over time (the insignificant trend in Table 3) as it was strongly correlated with the annually varying runoff (r=0.74, p<0.001, df = 28, Fig. S5b) but not with N input (r=0.22, p=0.24, df=28, Fig. S5c), which had a significant decreasing trend over the period from 1990–2019 (Table 2). The ATNL for the whole of Finland correlated with both N input and runoff. Table 3. Mann-Kendall trend tests for simulated agricultural TN loading for different Baltic Sea sub-basins, mean annual leached N fraction and the runoff for mineral fertilizers, and the total N input to soil. Annual mean runoff was used as a covarying environmental factor in the partial Mann-Kendall tests. Response variable Covariate Test statistics* p-value Whole Finland ATNL Runoff -2.50 0.013 Archipelago Sea Runoff -1.70 0.089 Gulf of Finland Runoff -2.66 0.008 Vuoksi Runoff -2.55 0.011 Bothnian Sea Runoff -1.98 0.048 The Quark Runoff -1.25 0.213 Bothnian Bay Runoff -2.58 0.010 Leached N Fraction Runoff -1.55 0.121 Mineral fertilizer - -0.73 < 0.0001 N input to soil - -0.66 < 0.0001 Runoff - 0.04 0.775 *Z for the Partial Mann-Kendall Trend tests (variables with covariates), tau for the Mann-Kendall Trend test (variables without covariates). Fig. 9. Mean leached fraction from the total N input to the agricultural fields for the whole of Finland during the period from 1990–2019 0.00 0.05 0.10 0.15 0.20 0.25 19 90 19 92 19 94 19 96 19 98 20 00 20 02 20 04 20 06 20 08 20 10 20 12 20 14 20 16 20 18 Le ac he d fr ac tio n Leached fraction Moving average (Leached fraction) Lin. (Leached fraction) Agricultural and Food Science (2023) 32: 112–127 124 Discussion This study shows that for the whole of Finland, simulated ATNL has decreased by 17%, from 17.4 kg ha-1 a-1 to 14.4 kg ha-1 a-1 since the 1990s. An earlier study of ATNL for the period from 1985–2006 has shown an increased ATNL trend, despite the agri-environmental measures implemented since the start of the Finnish Agri-Environmental Program in 1995 (Ekholm et al. 2015). Rankinen et al. (2016) analyzed P and N loading trends for the period from 1985–2012, showing that the TN load from agriculture increased until 2000–2006, but then started to decrease. Clearing of new fields explained 50% of the increase in the TN load to the Baltic Sea between 1995–1999 and 2000–2006 (Rankinen et al. 2016). In these previous studies, only 20 river basins were included, which covered 30% of agricultural fields, whereas in this study we simulated all Finnish watersheds. In this study the simulated ATNL was most affected by mineral fertilizer input, which has reduced remarkably since 1990 but stabilized in the past decade. There is a question about whether N leaching in the real world also depends on mineral fertilizer rates. Higher N fertilizer rates are known to lead to higher N surpluses, and these are associated with an increasing risk of N leaching in Finland (Salo and Turtola 2006). N leaching initially remains stable before following a linear relationship with N surplus, suggesting that 57% of N surplus would be leached. The higher the expected baseline of N leaching, the higher the N surplus required before it starts to affect N leaching (Salo et al. 2013). It is apparent that balancing N input with plant uptake is needed to prevent high NO3 --N leaching (Bergström and Brink 1986). The reasons behind the decrease of N surplus since 1990s and the decrease of simu- lated ATNL is demonstrated in this study. Further decreases of N fertilization in catchments with low N surpluses might not lead to further N leaching decreases. However, the mechanism of N leaching is complicated because only a minor fraction of leached N comes from fertilizer before it has passed through the immobilization-mineralization cycle. A high amount of the NO3 --N susceptible to leaching comes from the mineralization of organic N (Yläranta et al. 1993). Denitrification (anaerobic loss of N2O, NO and N2 from soils) is difficult to measure. McNeill et al. (2005) suggested losses of 20–40% of N applied. Our simulated average denitrification for two selected river catchments (Aurajoki and Vantaanjoki) ranged from 26% to 35%, which corresponds with these estimates. Surey et al. (2020) reports that the denitrification potential, defined as the combined amount of N released as N2O + N2, was larger for plots receiving farmyard manure than pure mineral fertilization due to the organic matter. The changing climate, particularly changes in runoff, may have a much greater effect on the TN transport out of the soils in the future. Longer periods without snow cover are increasing TN transport from the soils. Elevated temperatures enhance the mineralization of organic N and increases NO3 --N leaching (Huttunen et al. 2021). Our simulations demonstrated that runoff had high annual variability but no clear trend over the past 30-years in Finland. Runoff is, however, expected to increase alongside the warming climate, which means that the leached N fraction of the total N input might also increase in forthcoming decades as the leached N fraction depends strongly on the runoff. Nevertheless, our study suggests that the 2006 IPCC default for leached N fraction (FrachLeach is 30% of N input) in greenhouse gas inventory is too high for Finnish conditions, and that a 30-year mean (i.e., 14.4% of N input is leached) would be a better FrachLeach estimate for Finland’s national greenhouse gas inventory for now as no clear increasing/decreasing trend could be detected. Special attention should be paid to agricultural practices on peat soils, since they are an important source of elevated N loading. Peat soils have naturally high organic N content, which is mineralized into mineral compounds which are then leached and taken up by plants. High N mineralization and the addition of mineral or organic fertilizers to peat soils causes elevated N leaching. Very high N leaching has been measured in peat soils, even at negative N balances (Lemola et al. 2000). TN loads seem to be particularly higher in thick peat field plots, as shown by Yli-Halla et al. (2022), where the mean TN loads during a hydrological year were 15.4 and 9.2 kg ha−1 from the thicker and thinner peat plots, respectively. TN leaching from peat soils very much depends on the hydrological conditions and dryness of the soil. Organic N mineralization is especially intensive during dry periods with a low water table, but it seems that under re-wetting conditions in peat soils, groundwater NO3 --N concentrations can decrease very quickly (Tiemeyer and Kahle 2014). N leaching from peat soils can also be mitigated by the right crop choices as leaching from grass fields is lower than from cereal crops (Huhta and Jaakkola 1993). Huhta and Jaakkola report 18 kg ha-1a-1 and 37 kg ha-1a-1 TN leaching measured in a Tohmajärvi experimental fields with grass and barley cultivation, respectively. The ICECREAM model simulates hydrology and N processes for mineral and organic soils separately. Mean simulated N leaching from organic fields for 2010–2019 in Finland is I. Huttunen et al. 125 24.4 kg ha-1 a-1, but from mineral fields, it simulates 12.3 kg ha-1 a-1. Further validation of the ICECREAM model hydrology and N leaching from organic soils against observations would be necessary in future work. The BSAP country allocated reduction target for Finland for TN is set at 3030 t a-1 (HELCOM 2013), the reference loading is the mean for the period from 1997–2003. Simulated ATNL has not decreased much compared to that reference period, therefore further reduction of the ATNL is required to achieve the BSAP targets. Mineral fertilizer input has decreased over recent decades in Finland but is still adequate according to crop yield statistics. Further mitigation of N leaching would require more advanced mineral fertilizer usage, like split applications to react yield potential of growing seasons, precision application techniques, site- and time specific N recommendations. The positive developments in avoiding excess mineral fertilizer input in recent decades is partly counteracted in animal husbandry cluster areas, where a considerable amount of organic N is applied as manure (39% of total N fertilizer on national scale) and it is partially mineralized when there is no crop uptake. The ATNL trends described in this article, however, contain several types of model uncertainty: input uncertainty, structure uncertainty (process descriptions), parameter uncertainty and technical uncertainty (Refsgaard et al. 2007). The uncertainty of the results is also increased by the chaining of the models. The hydrological sub-model has its own model structure and parameter uncertainties related to its ability to simulate heterogeneity of the catchments. The main uncertainties in the ICECREAM model are caused by insufficiencies in the process description, and input data uncertainty. Mineral fertiliser inputs at the farm level are not available for the models as input data, even if this information is collected by regional authorities for the implementation of the agri-environmental subsidy programmes. The VEMALA model itself has model structure and parameter uncertainties simulating P and N transport and retention in water bodies. However, the use of deterministic models, which are based on physical process descriptions, enables us to produce consequent and consistent changes in output variables runoff and TN loading. We believe that the consistency of the results and validation examples are the basis for further use of the results by decision makers in water management work. Conclusions In this novel study, we conducted a detailed spatial simulation of ATNL on a national scale. Our findings revealed an average decrease of 17% in simulated ATNL since the 1990s. This decline could potentially be attributed to a decrease in mineral fertilizer input to the soils. However, validating these statements poses a challenge due to the lack of national-scale ATNL observations and the complexity of the nitrogen cycle in soils. Interestingly, we observed variations in the magnitude of the decreasing trends between the Finnish catchments, influenced by differences in hydrology, agricultural practices, and catchment characteristics. These results indicate the need for region-specific assessments and highlight the potential for using simulated national-scale ATNL estimates in compiling data for the HELCOM, as well as implementing the WFD through regional Environmental centres. To further enhance the accuracy of ATNL simulations, it is crucial to undertake additional research aimed at improving the models that simulate the intricate N cycle in agricultural soils. Such studies will help refine our understanding of the factors influencing ATNL and contribute to more effective environmental management strategies. Further joint efforts are needed to create a national fertilizer database collecting field scale fertilizer application data. Acknowledgements This work was funded by Ministry of Agriculture and Forestry. We are also thankful to the BlueAdapt-project (Con- tract No. 312650 BlueAdapt) funded by the Strategic Research Council of the Academy of Finland for supporting the VEMALA-ICECREAM model development. We thank Marie Korppoo and Noora Veijalainen for their valuable comments on the manuscript, and Nasim Fazel for help with the data visualisation. References Abrahamsen, P. & Hansen, S. 2000. Daisy: an open soil-crop-atmosphere system model. Environmental Modelling & Software 15: 313–330. https://doi.org/10.1016/S1364-8152(00)00003-7 Bechmann, M., Greipsland, I., Riley, H. & Eggestad, H. 2012. Nitrogen losses from agricultural areas - a fraction of applied ferti- lizer and manure (FracLEACH). Bioforsk Report nr 7. Agricultural and Food Science (2023) 32: 112–127 126 Bergström, L., & Brink, N. 1986. Effects of differentiated applications of fertilizer N on leaching losses and distribution of inorganic N in the soil. Plant and Soil: 333–345. https://doi.org/10.1007/BF02374284 Chatskikh, D., Olesen, J.E., Berntsen, J., Regina, K. & Yamulki, S. 2005. Simulation of effects of soils, climate and management on N2O emission from grasslands. Biogeochemistry 76: 395–419. https://doi.org/10.1007/s10533-005-6996-8 COUP manual. Coupled heat and mass transfer model for soil-plant-atmosphere systems (Jansson, P.E. & Karlberg, L. eds.) https://drive.google.com/file/d/0B0-WJKp0fmYCZ0JVeVgzRWFIbUk/view?resourcekey=0-6vzr2awC6k1x4jVfF2ReCw. Accessed: 8 April 2022. Ekholm, P., Rankinen, K., Rita, H., Räike, A., Sjöblom, H., Raateland, A., Vesikko, L., Cano Bernal, J.E. & Taskinen, A. 2015. Phos- phorus and nitrogen fluxes carried by 21 Finnish agricultural rivers in 1985-2006. Environmental Monitoring Assessment 187: 17 p. https://doi.org/10.1007/s10661-015-4417-6 Finér, L., Lepistö, A., Karlsson, K., Räike, A., Tattari, S., Huttunen, M., Härkönen, L., Joensuu, S., Kortelainen, P., Mattsson, T., Piirain- en, S., Sarkkola, S., Sallantaus, T. & Ukonmaanaho, L. 2020. Metsistä ja soilta tuleva vesistökuormitus 2020 - MetsäVesi-hankkeen loppuraportti. Valtioneuvoston selvitys- ja tutkimustoiminnan julkaisusarja 2020: 6. 77 p. (in Finnish). Government Decree 2014. Government Decree on Limiting Certain Emissions from Agriculture and Horticulture. 14 p. Grönroos, J. 2014. Maatalouden ammoniakkipäästöjen vähentämismahdollisuudet ja –kustannukset (Reduction possibilities and costs of agricultural ammonia emissions). Ympäristöministeriön raportteja 26. 92 p. (in Finnish). https://helda.helsinki.fi/han- dle/10138/152766 Grönroos, J., Mattila, P., Regina, K., Nousiainen, J., Perälä, P., Saarinen, K. & Mikkola-Pusa, J. 2009. Development of the ammonia emission inventory in Finland. Finnish Environment Institute. The Finnish Environment 8/2009. 60 p. Grönroos, J., Munther, J. & Luostarinen, S. 2017. Calculation of atmospheric nitrogen and NMVOC emissions from Finnish agri- culture - Description of the revised model. Finnish Environment Institute. 60 p. HELCOM 2013. Summary report on the development of revised Maximum Allowable Inputs (MAI) and updated Country Allocat- ed Reduction Targets (CART) of the Baltic Sea Action Plan. https://helcom.fi/wp-content/uploads/2019/08/Summary-report-on- MAI-CART-1.pdf. Accessed: 25 April 2023. Huhta, H. & Jaakkola, A. 1993. Viljelykasvin ja lannoituksen vaikutus ravinteiden huuhtoutumiseen turvemaasta Tohmajärven huuhtoutumiskentällä v. 1983-87. (The effect of crop and fertilization on nutrient leaching from peat soil in Tohmajärvi leaching field 1983-1987. Maatalouden tutkimuskeskus. Tiedote 20/91. 66 p. (in Finnish). Huttunen, I., Huttunen, M., Piirainen, V., Korppoo, M., Lepistö, A., Räike, A., Tattari, S. & Vehviläinen, B. 2016. A national scale nutrient loading model for Finnish watersheds - VEMALA. Environmental Modeling & Assessment 21: 83–109. https://doi. org/10.1007/s10666-015-9470-6 Huttunen, I., Hyytiäinen, K., Huttunen, M., Sihvonen, M., Veijalainen, N., Korppoo, M. & Heiskanen, A.-S. 2021. Agricultural nutri- ent loading under alternative climate, societal and manure recycling scenarios. The Science of the Total Environment 783. 15 p. https://doi.org/10.1016/j.scitotenv.2021.146871 Knisel, W. 1993. GLEAMS: Groundwater loading effects of agricultural management systems. Version 2.10. Publication No. 5. Athens, Georgia: University of Georgia, Department of Biological and Agricultural Engineering, Coastal Plain Experiment Station. Korppoo, M., Huttunen, M., Huttunen, I., Piirainen, V. & Vehviläinen, B. 2017. Simulation of bioavailable phosphorus and nitrogen loading in an agricultural river basin in Finland using VEMALA v.3. Journal of Hydrology 549: 363–373. https://doi.org/10.1016/j. jhydrol.2017.03.050 Kämäri, M., Huttunen, I., Valkama, P., Huttunen, M., Korppoo, M., Tattari, S. & Lotsari, E. 2019. Modelling inter- and intra-annual variation of riverine nitrogen/nitrate losses from snowmelt-affected basins under agricultural and mixed land use captured with high-frequency monitoring. CATENA 176: 227–244. https://doi.org/10.1016/j.catena.2019.01.019 Lemola, R., Turtola, E. & Eriksson, C. 2000. Undersowing Italian ryegrass diminishes nitrogen leaching from spring barley. Agricul- tural and Food Science in Finland 9: 2101–215. https://doi.org/10.23986/afsci.5661 Libiseller, C. & Grimvall, A. 2002. Performance of partial Mann-Kendall tests for trend detection in the presence of covariates. Environmetrics 13: 71–84. https://doi.org/10.1002/env.507 Lilja, H., Uusitalo, R., Yli-Halla, M.J., Nevalainen, R., Väänänen, T., Tamminen, P. & Tuhtar, J. 2017. Suomenn maannostietokanta: Käyttöopas versio 1.1. . Translated title of the contribution: Geographical soil database of Finland: Users’ manual version 1.1. Luonnonvara- ja biotalouden tutkimus 6. 70 p. Luke 2021. Official Statistics of Finland (OSF): Number of Livestock Helsinki: Natural Resources Institute Finland. http://stat.luke. fi/en/number-of-livestock (Accessed: 16 May 2021). Current URL: https://www.luke.fi/en/statistics/number-of-livestock Luke 2023a. CAP indicators / Nitrogen and phosphorus balance. Helsinki: Natural Resources Institute Finland. https://www.luke. fi/en/statistics/indicators/cap-indicators/nitrogen-and-phosphorus-balance. Accessed: 27 April 2023. Luke 2023b. Sales of fertilizers to farms. Data table in: CAP indicators / Nitrogen and phosphorus balance. Helsinki: Natural Re- sources Institute Finland. https://www.luke.fi/en/statistics/indicators/cap-indicators/nitrogen-and-phosphorus-balance. Accessed 27 April 2023. Luostarinen, S., Grönroos, J., Hellstedt, M., Nousiainen, J. & Munther, J. 2017. SUOMEN NORMILANTA - laskentajärjestelmän kuvaus ja ensimmäiset tulokset. Luonnonvara- ja biotalouden tutkimus 47/2017. Luonnonvarakeskus. Helsinki. 54 p. Mavi, 2009. Opas ympäristötuen ehtojen mukaiseen lannoitukseen 2007-2013. Maaseutuvirasto, Helsinki. Maaseutuviraston julkaisusarja: Hakuoppaita ja ohjeita. 27 p. McNeill, A.M., Eriksen, J., Bergström, L., Smith, K.A., Marstorp, H., Kirchmann, H. & Nilsson, I. 2005. Nitrogen and sulphur man- agement: challenges for organic sources in temperate agricultural systems. Soil Use and Management 21: 82–93. https://doi. org/10.1079/SUM2005303 I. Huttunen et al. 127 Rankinen, K., Keinänen, H. & Bernal, J.E. 2016. Influence of climate and land use changes on nutrient fluxes from Finnish rivers to the Baltic Sea. Agriculture, Ecosystems & Environment 216: 100–115. https://doi.org/10.1016/j.agee.2015.09.010 Refsgaard, J.C., van der Sluijs J.P., Højberg, A.L. & Vanrolleghem, P.A. 2007. Uncertainty in the environmental modelling process - A framework and guidance. Environmental Modelling & Software 22: 1543–1556. https://doi.org/10.1016/j.envsoft.2007.02.004 Rekolainen, S., Kämäri, J., Hiltunen, M. & Saloranta, T. 2003. A conceptual framework for identifying the need and role of models in the implementation of the Water Framework Directive. International Journal of River Basin Management 1: 347–352. https:// doi.org/10.1080/15715124.2003.9635217 Rekolainen, S. & Posch, M. 1993. Adapting the CREAMS model for Finnish conditions. Nordic Hydrology 24: 309–322. https://doi. org/10.2166/nh.1993.10 Salo, T., Turtola, E., Virkajärvi, P., Saarijärvi, K., Kuisma, P., Tuomisto, J., Muurinen, S. & Turakainen, M. 2013. Nitrogen fertiliser rates, N balances, and related risk of N leaching in Finnish agriculture. MTT Report 102. 37p. Statistics Finland 2021. Greenhouse gas emissions in Finland 1990 to 2019. National Inventory Report under the UNFCCC and the Kyoto Protocol. 15 March 2021. Sundholm, S. 2021. Simulating runoff, erosion, and phosphorus transport from a field with the ICECREAM model. Master’s The- sis 2021. Aalto University. 79 p. Surey, R., Lippold, E., Heilek, S., Sauheitl, L., Henjes, S., Horn, M.A., Mueller, C.W., Merbach, I., Kaiser, K., Böttcher, J. & Mikutta, R. 2020. Differences in labile soil organic matter explain potential denitrification and denitrifying communities in a long-term ferti- lization experiment. Applied Soil Ecology 153: 103630. https://doi.org/10.1016/j.apsoil.2020.103630 Swedish Environmental Protection Agency 2021. National Inventory Report Sweden 2021. Greenhouse Gas Emission Inventories 1990 to 2019. Submitted under the United Nations Framework Convention on Climate Change and the Kyoto Protocol. Tattari, S., Bärlund, I., Rekolainen, S., Posch, M., Siimes, K., Tuhkanen, H.-R. & Yli-Halla, M. 2001. Modeling sediment yield and phosphorus transport in Finnish clayey soils. Transactions of the ASAE 44: 297–307. https://doi.org/10.13031/2013.4691 Tiemeyer, B. & Kahle, P. 2014. Nitrogen and dissolved organic carbon (DOC) losses from anartificially drained grassland on organic soils. Biogeosciences 11: 4123–4137 Trolle, D., Nielsen, A., Andersen, H.E., Thodsen, H., Olesen, J. E., Børgesen, C. D., Refsgaard, J. Chr., Sonnenborg, T.O., Karls- son, I. B., Christensen, J.P., Markager, S. & Jeppesen, E. 2019. Effects of changes in land use and climate on aquatic ecosystems: Coupling of models and decomposition of uncertainties. Science of The Total Environment 657: 627–633. https://doi.org/10.1016/j. scitotenv.2018.12.055 Yli-Halla, M., Lötjönen, T., Kekkonen, J., Virtanen, S., Marttila, H., Liimatainen, M., Saari, M., Mikkola, J., Suomela, R. & Joki- Tokola, E. 2022. Thickness of peat influences the leaching of substances and greenhouse gas emissions from a cultivated organic soil. Science of The Total Environment 806: 150499. https://doi.org/10.1016/j.scitotenv.2021.150499 Yläranta, T., Uusi-Kämppä, J. & Jaakkola, A. 1993. Leaching of nitrogen in barley, grass ley and fallow lysimeters. Agricultural Science in Finland 2: 281–291. https://doi.org/10.23986/afsci.72651 National-scale nitrogen loading from the Finnish agriculturalfields has decreased since the 1990s Introduction Materials and methods Model description VEMALA-ICECREAM modelling system description Agricultural N loading simulation with ICECREAM National scale input data for the model Mineral fertilizer data Manure data Database of agricultural field characteristics Results Calibration of the denitrification parameter N sub-processes in agricultural soils Verification of simulated TN concentrations against river concentrations Spatial variation of agricultural TN loading Agricultural TN loading results for 1990–2019 Discussion Conclusions Acknowledgements References