A century of nitrogen dynamics in agricultural watersheds of Denmark

Intensive agriculture has been linked to increased nitrogen loads and adverse effects on downstream aquatic ecosystems. Sustained large net nitrogen surpluses have been shown in several contexts to form legacies in soil or waters, which delay the effects of reduction measures. In this study, detailed land use and agricultural statistics were used to reconstruct the annual nitrogen surpluses in three agriculture-dominated watersheds of Denmark (600–2700 km2) with well-drained loamy soils. These surpluses and long-term hydrological records were used as inputs to the process model ELEMeNT to quantify the nitrogen stores and fluxes for 1920–2020. A multi-objective calibration using timeseries of river nitrate loads, as well as other non-conventional data sources, allowed to explore the potential of these different data to constrain the nitrogen cycling model. We found the flux-weighted nitrate concentrations in the root zone percolate below croplands, a dataset not commonly used in calibrating watershed models, to be critical in reducing parameter uncertainty. Groundwater nitrate legacies built up in all three studied watersheds during 1950–1990 corresponding to ∼2% of the surplus (or ∼1 kg N ha yr−1) before they went down at a similar rate during 1990–2015. Over the same periods active soil nitrogen legacies first accumulated by approximately 10% of the surplus (∼5 kg N ha yr−1), before undergoing a commensurate reduction. Both legacies appear to have been the drivers of hysteresis in the diffuse load at the catchments’ outlet and hindrances to reaching water quality goals. Results indicate that the low cropland surpluses enforced during 2008–2015 had a larger impact on the diffuse river loads than the European Union’s untargeted grass set-aside policy of 1993–2008. Collectively, the measures of 1990–2015 are estimated to have reset the diffuse load regimes of the watersheds back to the situation prevailing in the 1960s.


Introduction
Reactive Nitrogen (N) is a key element for life, with a biogeochemical cycle spanning the atmosphere, the biosphere on land, soils, ground-and surface water, and the marine environment.The major increase in production of industrial fertilizers since the second half of the 20th century has been a major shift in the N cycle (Sutton et al 2011).Intensive agricultural practices, running with a large net N surplus, are prone to increasing loads to downstream components of the cycle, leading to adverse effects on aquatic environments (Sutton et al 2011).The accumulated N along the whole cascade contributes to time lags between input reductions and water quality improvements in downstream water bodies (e.g., Van Meter et al 2017a, Basu et al 2022).
Several uncertainties on the quantification of N fluxes and stores on land persist owing to the complexity of landscapes and the equifinality of several processes (Sarrazin et al 2022), e.g. it can be difficult to assert if the gap between the N field surplus and the river load under a certain period is due to permanent removal by denitrification or to delays in transit via biogeochemical or hydrological legacies (Dupas et al 2020, Basu et al 2022).As such, estimating the effects of agricultural practices and policies cannot be limited to calculating the associated net N surplus, but will also require to build an understanding of (i) the volatilization and atmospheric deposition to and from the atmosphere and lands neighboring agricultural areas, (ii) how much gets temporarily accumulated into soils, (iii) how much drains to surface water or leaches into groundwater, (iv) how long is the water travel time in the surface and groundwater systems and (v) how prevalent are the conditions leading to denitrification in those systems (e.g.Hansen et al 2019, Kumar et al 2020)).Obtaining the correct, dynamic picture of the stores and fluxes of reactive N at the scale of a landscape or a watershed thus requires an in-depth knowledge of both the farming and water systems, temporally and spatially (Cellier et al 2011).One way to build confidence in this understanding of the regional N cycle is to study its evolution over long time periods and over shifts in land use, agricultural practices and N input levels (Van Meter et al 2017a, 2017b, 2018, 2021, Ilampooranan et al 2019, 2022, Ascott et  Agriculture has historically been and remains a key land use in Denmark with around 59% of the land area being under crop rotations in 2021 (Statistics Denmark 2022a, 2022b).Eutrophication fueled by N pollution became a national problem in the 1980s, and led to the adoption of an extensive monitoring program of soil, drain and stream waters, and a series of national action plans (e.g.Windolf et al 2012, Blicher-Mathiesen et al 2014, Hansen et al 2019).With long agricultural and hydro-climatic records, and more than three decades of specific data gathering and experiences with remediating N pollution, Denmark offers great opportunities to study the N cycle in landscapes dominated by intensive agriculture, to identify the sources and sinks of reactive N, as well as potential legacy effects, and to understand the effect of different measures on the pollution loads.
In this study, with river water chemistry records dating back to the late 1960s or mid-1970s and a century of detailed agricultural statistics, using parsimonious models of N dynamics, we propose to summarize and test our long-term understanding of the N cycle and potential legacies at the watershed scale (∼1000 km 2 ) in three regions of Denmark dominated by agriculture.Taking advantage of the vast soil water chemistry monitoring effort under arable land since the early 1990s (Blicher-Mathiesen et al 2014, 2023) we investigate how these data on N-leakage at the root zone complement the data on N loads collected at the river outlets, to inform our understanding of N fluxes and stores, and we present historical reconstructions in three Danish watersheds for the past century.We also propose estimates of the effects of two groups of policies to reduce N leakage implemented in Denmark, namely: the European Union's mandatory set-aside schemes from 1993 to 2008, and the combined strict measures on fertilization limits, increased nitrogen utilization in manure and catch crop requirements in place nation-wide between 2008 and 2015.

Study areas
With a low elevation and flat topography, the climate of Denmark is mainly influenced by dominant westerly winds (Cappelen and Jørgensen 1999) and eastward moisture transport from the North Sea, giving rise to an east-west gradient of specific discharge (figures 1 and 2).The seasonal time-series of discharges and nitrate river loads (figure 2) show that the main runoff/leaching season (fall and winter) is separated from the main growing season (spring and summer).One can thus use the N surplus at the end of the growing season to model the current year's total potential for legacies and leaching, and then use the annual runoff to allocate this total as either leaching or legacies.This simplification of the seasonal dynamics, allows to work with the annual time series published in agricultural yearbooks, which go back more than a century in Denmark, and to capture the effect of processes with time-scales in the range of a decade or more.
Under these assumptions, potential study areas had to meet several conditions: • documented agricultural statistics, • streamflow records, • long time-series of N riverine loads, • measurements for N concentrations in root zone percolate under cultivated fields, • measurements of topsoil N content.
These constraints led to the selection of three gauged watersheds: the Suså River, the Odense River and the River Gudenå (figure 1 and table 1) over the period 1920-2020.These watersheds are fairly spatially homogeneous in terms of landscape, farming structures and practices.Agricultural specialization has led to a dominance of arable farms (∼75% of total agricultural area-TAA) in the region of Zealand and the Suså river watershed, combined with low overall livestock density over the latest decade: ∼8% and ∼4% of TAA allocated to pig farms and dairy farms respectively (table 1, Statistics Denmark 2022c).The province of Funen including the Odense River watershed and the province of East Jutland including the River Gudenå watershed are both characterized by approximately 50% of TAA dedicated to arable farms, 20% to pig farms, and 15% to dairy and cattle farms (table 1, Statistics Denmark 2022c).

Landscape model of the N cycle
This study makes use of the ELEMeNT-N model developed and introduced by Van Meter et al (2017a).ELEMeNT-N encompasses the major pools and fluxes of N on land, simplified into two subsystems (figure 3).The soil system goes down to the end of the root zone and to potential tile-drainage systems (∼1 m in these regions); it is discretized in parcels of land (0.1% of the watershed area) that together follow the historical evolution of key land use classes: cropland in rotation (some of it set aside between 1993 and 2008), permanent grasslands and nonagricultural lands (figure 4).A specific N surplus calculation is carried out for each land use and watershed (figure 5), and each parcel of land is characterized by two pools of Soil Organic N (SON), active and protected (aSON and pSON), and one pool of dissolved soil N (chiefly nitrate, figure 3).The surface water and shallow groundwater systems, together referred to as the watershed system, are described by parametric storage/transfer functions (Van Meter et al 2017a, Sarrazin et al 2022).Nitrogen point sources are added at the outlet to the diffuse load from the watershed system to express the total modeled river nitrate load that is to be compared to the river load calculated from the monitoring data (figure 3).

Constraining datasets and model ensemble selection
The parameter values to obtain a realistic representation of the N cycle in each watershed were not known a priori but had to be estimated by inverse modeling,  ] 1920 1940 1960 1980 2000 20201920 1940 1960 1980 2000 20201920 1940 1960 1980  Annual specific discharge, river load (Nitrate + nitrite and Ammonium + Ammonia), cropland average root zone nitrate flux-weighted concentration in selected watersheds and associated intensively monitored arable land mini-catchments for a water year spanning June 1st-May 31st, and distributions of monthly specific discharge and monthly area-normalized nitrate river load.Data sources: river discharge and river load (GRDC 2022, DCE 2022, Thodsen et al 2023), root zone nitrate concentration (DCE 2022, Blicher-Mathiesen et al 2023).The N load in Danish rivers is largely dominated by oxidized inorganic N (nitrate and nitrite), accounting on average for 80% of total N river load across all rivers for the period 1995-2020 (HELCOM 2022).
i.e., combinations of values had to be tested randomly first and then selected or discarded based on the associated model response (figure 3).For each tested input combination, the model results X m were compared to the corresponding monitoring dataset Four objective functions, i.e., composite measures of goodness-of-fit, were built to estimate the added information carried by C c and SON c (on top of J T ) to constrain the N cycling model: Table 1.Characteristics of the studied watersheds, associated monitored agricultural mini-catchments (LOOP) and regions used for N surplus calculations.Land fractions and runoff (mean±standard deviation) are given as averages for the period 1990-2018 and relative to the total land area.The data on LOOP catchments is reproduced from Petersen et al (2021).Tile-drainage coverage is calculated from Møller et al (2018). 1 LU (Livestock unit) is equivalent to 100 kg N yr −1 in manure ex-storage (i.e.net of losses in stables and storage), which corresponds approximately to the excretion of one dairy cow.For the time-coverage of data series of nitrate concentrations in soil water and river loads, see figure 2 The number of samples for grasslands and woodlands SON was quite limited in the regions of interest with only 8 and 13 samples respectively (figure S4(a)), which is why they were not used.For cropland SON with data for 2009, 2015 and 2018, only the mean bias indicator β [SON c ] was used.The use of ρ [C c ] is discussed in the next section.
Uniform distributions of the input parameters spanning a priori ranges (table S1) were randomly sampled following a GLUE-like methodology (Beven and Binley 1992) until 1000 models surpassed a chosen threshold for each watershed and objective function (requiring between 1.5 and 5 million random trials).The parameter search was implemented using the model-independent software OSTRICH (Matott 2017).The thresholds were set to 0.9 for F JT , 0.85 for F JT+SONc , 0.735 for F JT+Cc and 0.7 for F JT+Cc+SONc .

Data-worth to reduce model equifinality
We were able to model all three watersheds well with KGE ranging between 0. no sizeable effect in comparison, unless in conjunction with C c (figure 6).
It should be noted that the best achievable values of ρ [C c ] varied among the studied watersheds from 0.8 for Gudenå, to 0.7 for Suså, to 0.5 for Odense (figure 6; optimal range [0.9-1]) while all the other performance metrics α and β [SON c ] could be brought into their optimal range of [0.9-1.1].Including ρ [C c ] in the objective function and lowering the selection threshold would have led to significantly poorer behavior on all other indicators, which is why its influence was analyzed outside the objective function, in post-analysis.Several considerations could explain that difficulty to reproduce the timing of the leaching flux: (i) the intrinsic limits of the simple soil N model structure, (ii) the errors caused by approximating the annual infiltration by the specific runoff across all land uses, (iii) the partial representativity of the proxy LOOP sites for the conditions in the wider watersheds (table 1).
Bivariate plots for the model ensembles selected using all calibration datasets can show parameter interaction effects and the equivalent N removal capacity under average runoff conditions in the soil Λ s , watershed Λ w and combined Λ T (figure 7 and supplementary material 3.1).Our results indicate a total removal Λ T of 80%-90% for the River Gudenå, 65%-85% for the Odense River and 70%-85% for the Suså River (figure 7(c)).The watershed removal Λ w seems fairly well understood for the River Gudenå with 70%-80% and for the Suså River with 55%-70%, but a little less clearly for the Odense River with 50%-70% (figure 7 The best models (selected with F JT+Cc+S ŌN and in the top 100 scores of ρ [C c ]) indicate mean lag times (within root zone/from root zone to outlet/in total) of 1.9/0.4/2.3 years for Gudenå, 5.0/1.0/6.0 years for Odense and 3.0/0.4/3.5 years for Suså (medians of best model ensembles).As a reference, Ehrhardt et al (2021) observed a median mode of N travel times of 5.4 years (using a log-normal distribution) for 166 catchments in France and Germany, with 70% of all catchments under 10 years.The relatively low values of mean transport times of nitrates through the studied watersheds fit the local contexts of semiconfined top aquifers with a shallow redox interface, overlaid by loamy soils with a high density of tiledrains and ditches (Højberg et al 2021).Under these circumstances, artificial drains are the main pathway for nitrates to surface waters (Hansen et al 2019) and the contribution of deep groundwater circulation to surface water discharge is both limited and likely to carry only relatively small amounts of nitrates The mean values for our selected model ensembles are −3/+ 12 kg N ha yr −1 (aSON/pSON/SON) for Gudenå, −5/+ 14/+ 9 kg N ha yr −1 for Odense, and −3/+ 9/+ 6 kg N ha yr −1 for Suså over the same time period.The ensembles of selected models seem to be biased toward excessive accumulation of pSON which is compensated by an underestimation of soil denitrification.Models in the upper left region of the lines of indetermination (lower ∆SON/higher soil denitrification) may be more realistic (figures S8-S10).

Assessment of nitrogen legacies
The period 1920-1949 was still characterized by a limited access to mineral N fertilizers, a low N field surplus (∼15-20 kg N ha yr −1 , figure 8(a)), and apparently a slight loss of SON (∼2 kg N ha yr −1 , figure 8(a)).This mining of the soil fertility can be understood as reaping the final benefits of the last land conversions to cropland operated up until the 1890s in the studied regions.
The period 1950-1969 was marked by an explosion of synthetic N fertilizers in Denmark causing a rapid increase in the average N field surplus (from 15 to 60-65 kg N ha yr −1 , figure 5) which opened an era of N legacy building in all three watersheds that lasted until the end of the 1980s (figure 8(a)).This extra provision of N seems to have gone primarily into building aSON legacies (∼5-10 kg N ha yr −1 ) in the 1950s and 1960s and shifted progressively towards more accumulation as pSON in the 1970 s and 1980s (reaching up to ∼20 kg N ha yr −1 , figure 8(a)).Significant accumulation of soil inorganic N is observed only in the soils of the Suså watershed (∼2 kg N ha yr −1 , figure 8(a)).An increase in the watershed N (shallow groundwater) occurred mainly in the 1970s (∼3 kg N ha yr −1 , figure 8(a)).Very little change in aSON was visible in the 1980s compared to the previous decades; the diffuse load kept increasing in absolute terms then though, and started increasing relatively as well (figures 8(a) and (b)).
Measures starting in the 1990s engaged a decrease in mineral fertilizer use and cropland N surplus first from 75-90 to 40-50 kg N ha yr −1 and later further down to ∼30-40 kg N ha yr −1 , which was unseen since ca.1960 (figure 5); this opened a period of legacy reduction and decrease of the diffuse load (figure 8(a)).Common responses across watersheds are a decrease of watershed N (∼1-2 kg N ha yr −1 ), a much-reduced accumulation of pSON, and a release of aSON (∼5-10 kg N ha yr −1 , figure 8(a)).Some loss of mineral soil N is also observed in the Suså watershed (figure 8(a)).
The final period 2016-2020 was inaugurated by the relaxation of strict rules on fertilization and a return to surplus conditions akin to those in place in the late 1990s and early 2000s or previously in the 1960s (∼35-50 kg N ha yr −1 ), with delayed decreases in leaching and diffuse load still manifesting (figure 8(a)).
For the period 1950-2015, the mean N retention, defined as the difference between surplus and diffuse load, averaged at 89% of the surplus for Gudenå, 79% for Odense and 84% for Suså ( figure 8(b) and table 2) which places all three watersheds in the higher end of both the [45-88%] range observed by Dupas et al (2020) for 16 catchments in western France in the period 1976-2015 and the [72 ± 16%] reported by Ehrhardt et al (2021) for 238 catchments in France and Germany.The relative stability of N retention through time hides large variations though in how it was distributed between denitrification and accumulation in soils, across periods of legacy building and reduction (table 2).Long-term, our results indicate that denitrification accounted for 58%-69% of the surplus on average and SON accumulation for 16-21%, with negligible changes (<1%) of mineral soil N and downstream watershed N, which is comparable to the 12%-18% of N accumulation in soils, and 58%-66% of denitrification reported by Sarrazin et al (2022) for the Weser River and tributaries in Germany for 1960-2015, but indicates a larger share of denitrification than reported by Liu et al (2021) for the Grand River in Canada for 1930-2016 for example (table 2).
The average watershed N surplus and the river load exhibit a hysteresis loop in all three river basins between 1960 and 2020 (figure 8(c)).The shift between the 1970s and 1980s follows the period of highest buildup of aSON and watershed N (and mineral soil N for Suså) while the shift around 2015 follows the largest releases of the same legacies (figures 8(a)-( c)).This points to some N legacies at work in shallow groundwaters (±2% of surplus or ∼1 kg N ha yr −1 ), in aSON (±10% of surplus or ∼4-6 kg N ha yr −1 ) and in mineral soil N (Suså; ±3% of surplus or ∼1-2 kg N ha yr −1 ) as the drivers of hysteresis.The observed ratios of approximately 80%-90% of the legacies occurring in soils are in line with the results of Liu et al (2021) for the Grand River (Canada) and Sarrazin et al (2022) for the Weser River (Germany).

Estimating the effects of policies
The selected models offered the opportunity to bracket the effect of certain policies by translating the effect of their implementation on the input data.For example, the adoption of mandatory EU set-aside schemes from 1993 to 2008 was explored by simulating a scenario where all the set-aside areas stayed in rotation and received the corresponding yearly average cropland N field surplus during that period.
Similarly, Denmark imposed very strict measures in the period 1999-2015 (Petersen et al 2021), in particular: • limiting fertilization rates at times to 80% of the economic optimum, • requirements for increased utilization of nitrogen in manure, • mandatory catch crops.The compiled cropland surplus curves revealed a noticeable dip between 2008 and 2015 (included) as a result of the measures, followed by a return to its earlier values after the suboptimal fertilization limits were lifted (figure 5).The mean cropland N surplus was 23 kg N ha −1 cropland yr −1 (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) in the Suså River basin before jumping back to 37 kg N ha −1 cropland yr −1 (2016-2020, figure 5).
Similarly, the mean cropland N surplus went from 33 to 58 kg cropland yr −1 in the Odense River basin and Gudenå 1920-1949 1950-1989 1990-2015 1950-2015 Diffuse from 46 to 68 kg cropland yr −1 in the basin of the River Gudenå (figure 5).Year 2018 was excluded from the calculated average surpluses post-2015 because it saw an exceptionally severe drought in Denmark (and Europe) and unusually high N surplus.The impact of the three measures was estimated by simulating a scenario where the cropland surplus in each region would have stayed constant and equal to the level in effect pre-2008 and post-2015 (figure 5).Permanent grass set-asides showed a load reduction stabilizing close to 6% (figure 9) on par with the percentage of set-aside cropland area (∼7%; table S2).Note that non-food set-asides were still grouped with the cropland category.The set-aside schemes were not only aimed at alleviating environmental issues but also at providing economic control on production and pricing mechanisms.Its implementation was not targeted to high-risk areas for leaching, and instead based on the total area that was set aside, leading largely to the least profitable plots being set aside.As such, targeted set-asides to low retention areas or similar measures like permanently vegetated buffer strips along streams and ditches should be expected to yield larger effects than presented here.
Pathways without the policy-led dips in cropland N surplus (table S2) showed reductions in diffuse river load at the end of the implementation period: approximately 23% (median of the best models' ensemble) for the River Gudenå, 22% for the Odense River, and 21% for the Suså River Values past 2020 are projections assuming a constant runoff equal to the mean value measured for the period 1990-2020.The circles show the estimate of the median effect of the policies at their strictest level had they not been interrupted.Note that the surplus in all three watersheds was strongly driven up by the drought in 2018 (figure 5).Since this was not reflected in the comparative scenarios built to estimate the maximal potential of the measures, one observes large jumps in the projected reductions in 2019.The same effect explains the overshoot projected in the Suså watershed after another marked drought in 2020.
(figure 9).Note that the reduction levels reached by the policy shown in figure 9 are expressed for the whole watershed areas, not only their cropland areas, which explains the lower reduction levels as compared to the cropland relative N surplus reduction.
Scenarios were also built to project the final effect of both groups of measures had they not been interrupted in 2008 or 2016 (figure 9).In all three watersheds, the measures show significant effects about three years after they were implemented, unlock half their final potential after 5-6 years, and would have reached near-full potential after about 10 years of strict implementation (figure 9).Both groups of measures showed a sustained, but diminishing effect after the implementation window, which is understood by a slow re-building of the SON stocks after the reset operated by periods of low N surpluses (figures 8 and 9).It is noteworthy that the Danish regulators after implementing nation-wide limits on fertilization and requirements on utilization of nitrogen in manure and catch crops, are now targeting further N load reduction by experimenting with catchmentspecific (<15 km 2 ) rules based on detailed mappings of zones with a low potential for denitrification and a high risk of leaching (Højberg et al 2021).Further modeling efforts could aim at integrating the catch crops requirements and the higher share of manure in fertilization as an effect on the model parameters of ELEMeNT-N, to separate them from the impact of low N surpluses.Other measures were indirectly at work during the same period, in particular a significant reduction of atmospheric N emissions and deposition (EMEP 2022; figure 5).These efforts were reflected in the data underlying the N surplus calculation and unequivocally contributed to the observed river load reductions, but their effect was not estimated separately.

Conclusions
Long-term hydrological records and detailed landuse and agricultural statistics were used to reconstruct a century of N trajectories in three agriculturedominated watersheds of Denmark.The use of annual measurements of nitrate leaching rates in representative mini-catchments helped build a robust understanding of nitrate reduction from the root zone to the outlet as compared to using only river load data (figures 6 and 7).The reconstructed trajectories of N fluxes show the impact of decades of large field N surpluses and the associated build-up of soil N legacies, and groundwater legacies to a lesser extent, which drove a non-linear response of the diffuse N loads reaching the rivers' outlet (figure 8).By implementing measures on maximum fertilization and catch crops, the Danish agricultural sector has shown that it is possible to reduce the soil N surpluses, resorb the active legacies and restore the regime of diffuse riverine load to the situation previously experienced in the 1960s (figure 8), while producing cropland yields higher by more than 30% (figure 5).Other recent assessments identified biogeochemical (SON) legacies as a major interfering link between implemented measures and missed water quality targets in watersheds dominated by modern agriculture (Dupas et al 2020, Ehrhardt et al 2021, Liu et al 2021, Sarrazin et al 2022), and as a priority area for more research and modeling in order to design adequate nitrogen management plans (Basu et al 2022, Lutz et al 2022, Golden et al 2023).

Figure 1 .
Figure 1.Map of selected watersheds, associated administrative regions and intensively monitored agricultural mini-catchments with soil water monitoring.
Figure2.Annual specific discharge, river load (Nitrate + nitrite and Ammonium + Ammonia), cropland average root zone nitrate flux-weighted concentration in selected watersheds and associated intensively monitored arable land mini-catchments for a water year spanning June 1st-May 31st, and distributions of monthly specific discharge and monthly area-normalized nitrate river load.Data sources: river discharge and river load(GRDC 2022, DCE 2022, Thodsen et al 2023), root zone nitrate concentration(DCE 2022, Blicher-Mathiesen et al 2023).The N load in Danish rivers is largely dominated by oxidized inorganic N (nitrate and nitrite), accounting on average for 80% of total N river load across all rivers for the period 1995-2020 (HELCOM 2022).
X o (where X can stand for the river nitrate load J T , the average cropland flux-weighted nitrate concentration in percolate C c at proxy sites, or the cropland average soil N-SON c , data extrapolated from the LUCAS soil survey by Orgiazzi et al (2018), see figure S4), using the decomposition of the Kling-Gupta model efficiency KGE (Gupta et al 2009) in terms of variability bias α [X] = σ [X m ] /σ [X o ] (with σ the standard deviation), mean bias β [X] = Xm / Xo ( X denoting the mean), and correlation ρ [X] with ρ the Pearson correlation coefficient of X m and X o .

Figure 3 .Figure 4 .
Figure 3. Principle representation of the ELEMeNT-N model (lefthand side) with four land use classes (cropland, set-aside, permanent grassland and non-agricultural land) and two compartments (soil and watershed), and inverse modeling workflow (righthand side).Areas can transition from one land use class to another to reflect the land use evolution and carry the N pools they acquired under past managements.See section 2.2.3 of the supplementary material and tableS1for more details about the model parameters.Gray boxes in the flow chart were handled by the OSTRICH software.

Figure 5 .
Figure 5. Breakdown and evolution of the N surplus by region (columns) and land use (rows).Solid lines show the net N surplus.Note the differences in y-scales.Dashed lines show an alternative scenario without strict policies on maximum fertilization and increased manure use in the period 2008-2015.See supplementary material 2.2.2 for details.

Figure 6 .
Figure 6.Distribution of the Damköhler numbers for nitrates in the soil and watershed compartments: 1000 selected models under each set of calibration constraints and watershed.

Figure 7 .
Figure 7. Bivariate distribution of (a) soil transport parameters, (b) watershed transport parameters, (c) soil and watershed N removal capacities for the models selected under F JT+Cc+SONc .Each model is represented by a marker; marker size, opacity and color are displayed as functions of ρ [Cc].Contour lines show soil, watershed and total removal capacities.

Figure 8 .
Figure 8. Modeled N fluxes and changes in storage by period (ensemble average of F JT+Cc+S ŌN with 100 highest ρ [Cc]), expressed as area-normalized fluxes (a) or percentages of the N surplus (b), and nitrate river load (diffuse and point sources) as a function of soil surplus (c).Fluxes in (a) and (b) are counted as positive while the sign of changes in storage can vary.The algebraic sum by column gives the mean surplus of the period.In (a) the surface area of each polygon is proportional to the magnitude of the accumulated flux or legacy.In (c) square markers show modeled ensemble-average by period while other markers show the individual monitoring data points that were used in the calibration.The upswing in the 1980s corresponds to the peak load from point sources, largely remediated by the implementation of N removal in wastewater treatment at the end of the decade (transition from λp, pre to λp,post).

Figure 9 .
Figure 9.Estimated impact of two policies on the diffuse river load.The shaded area shows the implementation period and the whisker boxes shows the spread of the top 100 models according to F JT+Cc+S ŌN and ρ[Cc].Values past 2020 are projections assuming a constant runoff equal to the mean value measured for the period 1990-2020.The circles show the estimate of the median effect of the policies at their strictest level had they not been interrupted.Note that the surplus in all three watersheds was strongly driven up by the drought in 2018 (figure5).Since this was not reflected in the comparative scenarios built to estimate the maximal potential of the measures, one observes large jumps in the projected reductions in 2019.The same effect explains the overshoot projected in the Suså watershed after another marked drought in 2020. .

Table 2 .
Comparison of diffuse load and retention expressed as percentages of surplus for the period.