Figure 1 . Study catchments (n = 238) based on the quality
criteria with selected catchment characteristics: a – Elevation (EEA,
2013), b – Land cover (CLC, 2000), c – Lithology (BCR & UNESCO,
2014), d – Depth to bedrock (Shangguan et al., 2017).
From this joint database we selected catchments where the following
conditions were given: riverine NO3-N concentration
observations available for at least 20 years of data with data gaps less
than 2 years and the total number of observations being more than 150.
Given these criteria, 238 catchments were selected (Figure 1a). The time
series covered data between 1971 and 2015 with a median length of 30
years and in total 96,443 measurements for NO3-N.
Overall we covered 40% of the total land area of both countries (i.e.,
around 361,000 km², taking nested catchments into account). The selected
catchments encompass contrasting settings in terms of morphology,
climate, geological properties and land use attributes (Supporting
Information Tables S1.1 and S1.2). More than half of the study
catchments have a size of less than 1,000 km² (max. 62,500
km2). The median altitude ranges from 15 m to 1848 m
with a median slope of 3°. Climatic settings of the sites reach from
Atlantic to Continental climate with aridity indices ranging between 0.4
and 1.5. The median annual precipitation across the sites is around
816 mm, and the estimated base flow index (BFI) ranges from 29% to 97%
with a median of 65%.
Most catchments (> 90%) are dominated by sandy soils
(median: 44.6%), with 18 of those located in northwestern Germany. The
bedrock mainly covers fissured and hard rock geology with the latter
being predominant in most of the catchments. The geology is
characterized by crystalline rocks in the Armorican Massif, the Pyrenees
and the Massif Central and in some of the German mountainous catchments;
and younger sedimentary rocks in most parts of France and Germany
(Allain, 1951; BGR & CGMW, 2005). Quaternary sediments are found in the
Northern German Lowlands, the Alpine foothills and north of the Pyrenees
(Allain, 1951; BGR & CGMW, 2005).
Regarding land use, 87% of the catchments had at least one-third of
their area covered by agriculture that mainly incorporates non-irrigated
arable land and pastures (EEA, 2016; Figure 1b). Riverine
NO3-N concentrations in these areas are therefore
predominantly impacted by diffuse agricultural N sources (EEA, 2018).
The median share of forest cover across the study catchments is 37%.
Although the fraction of artificial surfaces was small, the median
population density with 92 inhabitants km-2 in the
study catchments is almost three-times the average European value
(Worldometers.info, 2020).
2.2. Nitrogen input
The N input was selected as diffuse N stemming from agricultural N
surplus, atmospheric deposition and biological fixation in
non-agricultural areas. The N surplus consists of agricultural N input
that is in excess of crop and forage exports (also known as land
nitrogen budget; de Vries et al., 2011). Here, we relied on two national
scale data sets. Agricultural N contribution and atmospheric N
deposition for the French catchments were provided by Poisvert et al.
(2017). The annual agricultural N surplus for German catchments was
provided by Bach and Frede (1998) as well as Häußermann et al. (2019).
It basically consists of two data sets available at a (coarser) state
level (NUTS2) for 1950–1999 and at finer county level (NUTS3) for
1995–2015. Both data sets were harmonized to produce a consistent
long-term data set. The atmospheric N deposition for German catchments
is based on Europe-wide gridded data from a chemical transport model of
the Meteorological Synthesizing Centre-West (MSC-W) of the European
Monitoring and Evaluation Programme (EMEP) (Bartnicky & Fagerli, 2006;
Bartnicky & Benedictow, 2017).
In agricultural areas, biological fixation was already included in the N
budgets. The biologically fixed N fluxes to non-agricultural land use
types for France and Germany were calculated using the European Corine
Land Cover data set from the year 2000 (EEA, 2020), which is most
representative regarding the water quality time series. Terrestrial
biological N mean uptake rates were set for forest (to 16.04 kg N
ha-1 yr-1; Cleveland et al., 1999),
for natural and urban grassland (to 2.7 kg N ha-1yr-1; Cleveland et al., 1999) and other land use
(wetlands, water bodies, open space with little or no vegetation to
0.75 kg N ha-1 yr-1; Van Meter et
al., 2017). A comparison of the two national long-term data sets for
diffuse N with a Europe-wide benchmark estimation for 1997–2003 (West
et al., 2014) indicated an acceptable offset (see Supporting Information
S2 for further information).
Due to the lack of spatially and temporally reliable long-term data on N
input by waste water, we did not consider this point source. For France,
Dupas et al. (2015) estimated the contribution from point sources to
total N flux to be 3% during the period 2005–2009, and we hypothesized
that the negligible contribution of point sources also held for Germany.
2.3. Nitrogen output as riverine NO3-N concentrations
and loads
Gaps in the discharge time series at 30 runoff stations in Germany were
filled through the support of simulations from the grid-based
distributed mesoscale hydrological model mHM (Kumar et al., 2013;
Samaniego et al., 2010). Here, only model simulations resulting in an
R2 greater than 0.6 when compared with the observed
discharge were accepted. A piecewise linear regression was utilized to
correct for potential biases in the modelled data. These bias-corrected
modelled discharge data were finally used to gap-fill the original data
to obtain a continuous daily time series. In France, no such national
hydrological model existed and therefore, we only included catchments
with nearly continuous daily discharge monitoring for which short gaps
in the discharge (max. 7 days) were interpolated by a fixed-interval
smoothing via a state-space model using the R software package
“Baytrends”.
The irregularly sampled, riverine NO3-N concentrations
were used to estimate daily concentrations by using the software packageExploration and Graphics for RivErTrends (EGRET) in the R
environment by Hirsch and DeCicco (2019). The applied Weighted
Regressions on Time, Discharge, and Season (WRTDS) uses a flexible
statistical representation for every day of the discharge record and has
been proven to provide robust estimates (Hirsch et al., 2010; Van Meter
& Basu, 2017). As we focus on changes in concentrations and fluxes
independent of inter-annual discharge variability (Hirsch et al., 2010),
we used flow-normalized concentrations and fluxes for further analyses.
For each catchment median annual flow-normalized NO3-N
concentrations and annual summed NO3-N fluxes were
calculated and scaled to the catchment area.
2.4. Nitrogen transport time
Travel time distributions are commonly derived as the transfer function
between rainfall concentration time series and stream concentrations of
a conservatively transported solute or water isotope (e.g. Kirchner et
al., 2000). We transfer this concept to reactive N transport with the N
input as an incoming time series with annual resolution that is assumed
to yield the median annual riverine NO3-N concentration,
when convolved with a fitted distribution. This transport time
distribution (TTD) can be based on different theoretical probability
distribution functions. To represent the long memory of past inputs,
long-tailed distributions are most suitable at catchment scales
(Kirchner et al., 2000). Therefore, the N input was convolved using a
log-normal distribution (Equation 1; Ehrhardt et al., 2019; Musolff et
al., 2017) to find the optimal fit to riverine NO3-N
concentrations. We alternatively used a gamma distribution (Equation 2;
Godsey et al., 2010; Fiori et al., 2009; Kirchner et al., 2000) as a
transfer function, and we compared the quality of fit
(R2) with both methods.
Equation 1 \(f\left(t\right)\) =\(\frac{1}{\text{tσ\ }\sqrt{2\ \pi}\ }\ exp\ (-\frac{{(ln\ t-\ u)\ }^{2}}{2\sigma^{2}}\ )\)
Equation 2\(f\left(t\right)=t^{-\alpha}\frac{\varepsilon^{-t/\beta}}{\beta^{\text{α\ }}G(\alpha)}\)
The two parameters mu (µ) and sigma (σ) for the log-normal and shape
(\(\alpha\)) and scale (β ) for the gamma distribution,
respectively, were calibrated through optimization based on minimizing
the sum of squared errors between the normalized annual diffuse N input
and normalized annual median riverine NO3-N
concentrations. For this purpose we used the Particle Swarm Optimization
(using the R package “hydroPSO” by Zambrani-Bigiarini & Rojas, 2013)
algorithm in 30 independent runs. We estimated the mode of the selected
best fitted TTD (with max. R2) to represent the peak
TT and at the same time to resemble the peak N export of the mobile,
inorganic N.
2.5. Nitrogen retention and its temporal change
The total cumulative diffuse N input load was compared to the respective
riverine NO3-N load (assumed as N load) to analyze the N
retention in the catchment (Equation 3). The difference between the two
is the load being retained in the catchment as biogeochemical legacy, as
hydrologic legacy or being removed by denitrification. The cumulative
flux differences were calculated based on two approaches: 1) using the
annual frames of the overlapping years in in- and outflux, while
disregarding time shifts; and 2) applying the derived TTs, to compare
the convolved inputs with the corresponding annual exported load.
Equation 3\(Retention=1-\ \frac{\text{Nout}}{\text{Nin}}=1-\frac{\sum_{i=ts}^{\text{te}}{NO3-N\ Flux\ i}}{\sum_{i=ts\ }^{\text{te}}\text{Ninput\ i}}\)
To further characterize the catchment’s reaction to N input changes, we
compared the median diffuse N input in the 1980s (median year of max. N
input: 1986) with the one in the last years of the time series (≥ 2010)
for a subset of stations (n = 120) that sufficiently covered the 1980s
and 2010s. The same was done with the exported riverine
NO3-N loads in the 1980s and the 2010s. To gain robust
estimates for the size of difference, we calculated the bootstrapped (n
= 10,000) median differences between the 1980s and 2010s (for N input
and N output) with their corresponding 95% confidence intervals.
2.6. Statistical analysis for controls in catchment response and
retention
We applied a Partial Least Squares Regression (PLSR) to identify the
main factors controlling N TTs and N retention in a catchment. PLSR is
an established multivariate regression approach to analyze data sets
that are strongly correlated among predictors and noisy (Wold et al.,
2001). The PLSR model finds the variables (catchment characteristics)
that best predict the response variables (retention and TT; Ai et al.,
2015). The importance of each predictor for the dependent variable is
indicated by the measure Variable Importance in the Projection(VIP). Factors with VIPs larger than 1 are considered to be
significantly important for explaining the dependent variable (Ai et
al., 2015; Shi et al., 2013). The corresponding regression coefficient
is used to explain the direction of influence of each independent
variable (Shi et al., 2013). The predictor variables used in this study
characterize the topography, land cover, climate, hydrology, lithology,
soils and population density of the studied catchments (Supporting
Information S1).
3. Results
3.1. Nitrogen transport time scales
Using the gamma distribution yielded comparable results to the results
for log-normal distribution (both with median R2 =
0.8), but less catchments with an acceptable fit (R2 ≥
0.6) between the convolved annual N inputs and riverine concentrations.
Therefore, we only report the results using a log-normal distribution as
a transfer function.
In some catchments (n = 72) no acceptable fit of TTDs could be obtained.
According to a Wilcoxon rank sum test, the variability in
NO3-N concentrations in these catchments (CV: 0.08) is
significantly different (p ≤ 0.01) to the ones in the other catchments
(CV: 0.12 with n = 166). A low temporal variability in the input or
output makes it challenging to derive a reliable transfer function
connecting them.
The median mode (peak) of the TTs for the 166 selected catchments with
an acceptable fit was 5.4 years (Supporting Information Table S3.).
Although the mode ranged from 0.2 to 34.1 years, the majority (70%) had
a mode TT less than 10 years (Figure 2c). Only a few catchments (10%)
showed a mode of at least 20 years, most of them (11/17) located in the
Massif Central (Figure 2a).
Although the TT derivation was not mass conform, on average across the
study catchments, 75% (75%-percentile) of the N input should have been
exported after 18 years (range: 1.4–38.2).
3. 2. Nitrogen retention
The median N retention of
the selected catchments (n = 238) was 72% (sd: 16%; Supporting
Information Table S3.; Figure 2b), meaning that a large part of N was
retained as legacy or denitrified. Despite the wide range (-24–96%,
with one negative outlier Figure 2b), 48% of the catchments had a
retention between 50% and 75%. A convolution of the N inputs according
to the corresponding TT resulted in a slightly lower retention with a
median of 70% (n = 238; 71% with n = 166; Figure 2d).
N retention and TT did not correlate in the study catchments. Almost the
same amount of catchments with retention above the median had TTs below
and above the median (Figure 2e).
The median diffuse N input in the 1980s was 62.6 kg N
ha-1 yr-1 (IQR: 42.0 kg N
ha-1 yr-1), decreasing by around
36%, when assuming the bootstrapped difference in medians of 22.6 kg N
ha-1 yr-1 (95% CI: 20.5–25.6 kg N
ha-1 yr-1) in comparison to the
2010s. Diffuse N input in the 2010s was around 38.4 kg N
ha-1 yr-1 (IQR: 23.1 kg N
ha-1 yr-1). The median N load in the
1980s was 12.4 kg N ha-1 yr-1 (IQR:
6.1 kg N ha-1 yr-1) with a
bootstrapped difference of medians of 1.2 kg N ha-1yr-1 (95% CI: 0.8–1.6 kg N ha-1yr-1) to the 2010s (median N load: 11.2 kg N
ha-1 yr-1; IQR: 5.7 kg N
ha-1 yr-1). The mismatch between N
input and riverine N export decreased from an annual excess of 50.2 kg N
ha-1 in the 1980s to 27.2 kg N ha-1in the 2010s, also reflecting a decrease in apparent retention from 80%
to 71%.