Long‐Term Nitrate Trajectories Vary by Season in Western European Catchments

Human alteration of nutrient cycles has caused persistent and widespread degradation of water quality around the globe. In many regions, including Western Europe, elevated nitrate (NO3−) concentration in surface waters contributes to eutrophication and noncompliance with environmental legislation. Discharge, NO3− concentrations and the vulnerability of the aquatic ecosystems to eutrophication often exhibit a distinct seasonality. Understanding spatial patterns and long‐term trends in this seasonality is crucial to improve water quality management. Here, we hypothesized that NO3− concentrations during high‐flow periods would respond faster to changes in nutrient inputs than low‐flow concentrations because of greater connectivity of shallow diffuse NO3− sources with the river network. To test this hypothesis, we compiled long‐term NO3− and discharge time series from 290 Western European catchments. To characterize the long‐term trajectories of seasonal NO3− concentration, we propose a novel hysteresis approach comparing low‐ and high‐flow NO3− concentration in the context of multi‐decadal N input changes. We found synchronous winter maxima of NO3− and discharge in 84% of the study catchments. However, contrary to our hypothesis, there were surprisingly diverse long‐term trajectories of seasonal NO3− concentration. Both clockwise (faster high‐flow NO3− response) and counterclockwise hysteresis (faster low‐flow NO3− response) occurred in similar proportions, potentially due to a high complexity in the underlying processes. Spatial variability of seasonality in NO3− concentration across the catchments was more pronounced and better predictable than its long‐term variability. This work demonstrates the value of seasonal and inter‐annual hydrochemical analysis and provides new tools for water quality monitoring and management.

Patterns of seasonal N concentration variability are controlled by various interacting physical, climatic, biological, and anthropogenic factors, which vary among catchments. First, the relative contributions of point and diffuse sources can affect concentration seasonality (Van Meter et al., 2020). Point source dominance often leads to highest concentrations during low flow, while diffuse N sources are often exported with higher discharges Van Meter et al., 2020). Second, hydrological connectivity can variably activate heterogeneously distributed solute sources within the catchment as discharge fluctuates at seasonal or event scales (Aubert et al., 2013;Dupas et al., 2016;Seibert et al., 2009). For example, the vertical distribution of solutes in the subsurface can control the concentration-discharge (C-Q) relationship (Botter et al., 2020;Zhi & Li, 2020) at event, seasonal, and inter-annual scales Moatar et al., 2017;Musolff et al., 2015;Zarnetske et al., 2018). Third, seasonal variations in hydroclimatic drivers and biogeochemical controls such as riparian or in-stream retention processes can create seasonal variations of riverine N concentrations (Casquin et al., 2020;Kumar et al., 2020;Lutz et al., 2020). For example, during summer low-flow conditions, biological N uptake and removal processes in the benthic and hyporheic zones can more strongly influence water column chemistry because residence times and temperature are higher (Moatar et al., 2017;Wollheim et al., 2017). Conversely, during high-flow periods, the larger flux of nutrients and the disturbance of aquatic nutrient uptake via scouring and sedimentation can result in minimal biological influence on stream chemistry (Blaszczak et al., 2019;Raymond et al., 2016). Among other things, this makes high-flow conditions the most important to nutrient flux, but summer low-flow conditions the most sensitive to eutrophication (e.g., Minaudo et al., 2020;Withers & Jarvie, 2008).
The seasonal variability of N concentrations can be characterized in terms of annual amplitude and timing of the minimum and maximum concentrations. Applied to multi-annual concentration time series this allows identification of long-term controls of seasonality, including strength and type of nutrient loading, which can change through time (Howden et al., 2010;Westphal et al., 2020). For example, decreases in point and diffuse source strength are expected to have greater influence on low-flow and high-flow concentrations, respectively . More generally, decreases in point sources may have an immediate effect while diffuse sources are likely to show time lags as described above Westphal et al., 2020). This delayed response can differ between seasons (Ehrhardt et al., 2019;Van Meter & Basu, 2017) according to changes in hydrological connectivity and transit time, with a larger contribution from younger water during high flow periods (Benettin et al., 2017;Yang et al., 2018). Accordingly, Dupas et al. (2016), Ehrhardt et al. (2019) and Winter et al. (2020) demonstrated how changes in N input loading to a catchment can alter long-term N export dynamics, seasonal amplitudes, and even reverse the seasonal timing, in some cases. In other cases, the timing of N seasonality was largely constant even under various concentration trends, potentially due to invariant vertical distributions of the nutrient source zones Dupas et al., 2018). Finally, changes in N seasonality can also be caused by changes in in-stream processes related to other factors, including P concentration (Bowes et al., 2011;Minaudo et al., 2015). As the inter-annual trajectories of low-and high-flow concentrations integrate the input history, retention and release processes at catchment scale, they can be a valuable characteristic and indicator of how catchments respond to changes in nutrient inputs.
Despite a growing body of work on N seasonality, the spatial variability of N concentration seasonality and its long-term trend are still poorly understood across catchments. In this context, we investigated spatio-temporal patterns of nitrate (NO 3 − ) seasonality in 290 catchments from two unique datasets of long time series across France and Germany. These catchments cover a wide range of ecohydrological and land use conditions, allowing robust quantification of seasonality and trend parameters on decadal timescales. We hypothesized that in these Western European landscapes where NO 3 − sources are primarily diffuse, most catchments would show (a) a NO 3 − maximum during the winter high-flow season and (b) a faster response of high-flow NO 3 − concentrations to changes in loading compared to low flow because of shallower N sources being affected first and being activated during high flow. To test the hypotheses, we propose a hysteresis approach for long-term trends of NO 3 − seasonality using available low frequency data. We classified the catchments based on their long-term average NO 3 − and discharge seasonality as well as long-term trajectories of NO 3 − seasonality and assessed links between these metrics and potential controls such as mean precipitation and fraction of agricultural land use. This large sample approach allowed us to determine archetypes of spatiotemporal NO 3 − export patterns at continental scale and evaluate underlying processes.

Databases
We used long-term riverine NO 3 − -N (NO 3 -N, hereafter) concentration data and daily discharge (Q) data from stations in Germany and France to analyze spatiotemporal patterns of NO 3 -N seasonality linking to export dynamics. From the original databases for Germany (Musolff, 2020;Musolff et al., 2020) and France (http://naiades.eaufrance.fr/; http://hydro.eaufrance.fr/; Dupas et al., 2019;Ehrhardt et al. (2021);Minaudo et al., 2019;Moatar et al., 2017) we selected stations based on the following criteria: (a) time series of NO 3 -N covering at least 20 years, (b) with at least 150 NO 3 -N samples, (c) a maximum sampling gap of 20% of the time series, (d) seasonal coverage, that is at least 10% of the data in each season (all possible quarterly divisions by months), and (e) available corresponding gauge with daily Q data. Discontinuous Q time series at 36 German stations were filled through the support of simulations from the grid-based distributed mesoscale hydrological model mHM (Kumar et al., 2013;Samaniego et al., 2010). Model results with a regression coefficient (R 2 ) greater than 0.6 with the observed Q were accepted. Next, we used a piece-wise linear regression to reduce bias in modeled Q data, used for gap-filling, following Ehrhardt et al. (2020). The criteria resulted in a set of 290 catchments (165 in France, 125 in Germany) with a median of 336 NO 3 -N samples per station and 114,479 in total, which cover various temperate ecoregions, topographic, climatic and land cover settings (Figures 1 and S2). The median time series length was 32 years and the maximum length 46 years starting in 1969. The majority of time series started around 1980 and ended in 2014 or 2015.
We reconstructed daily concentrations (C) and flow-normalized concentrations (C FN ) using Weighted Regression on Time, Discharge and Season (WRTDS, Hirsch et al., 2010) as implemented in the R package EGRET (version 3.0.2, Hirsch & De Cicco, 2015). Flow-normalization estimates are calculated by averaging C estimates from all observed Q values of the specific day of the year throughout the time series. To focus on concentration trends independent of inter-annual discharge variations we use C FN for our analysis. For data gaps of observed NO 3 -N larger than two years, corresponding daily C FN estimates were excluded due to high uncertainties, following Hirsch and De Cicco (2015). Observed and interpolated time series of the Wupper catchment are shown in Figures 2a-2c as an example.

Low-and High-Flow Seasons and Average Nitrate Seasonality
We determined the high-flow (HF) and low-flow (LF) seasons for each catchment by calculating the longterm average discharge for three consecutive months and defined the HF and LF seasons as the three wettest and driest months, respectively. Based on the timing of HF we classified the catchments into winter maximum (winMax, center of three wettest months within November-March) and spring or early summer  (Table S1, Cornes et al., 2018), (c) mean NO 3 -N concentrations and land cover classes (CLC, 2006), (d) three selected catchments with mean NO 3 -N, land cover, and N surplus time series. N surplus is the diffuse N input from agriculture, atmospheric deposition, and biological fixation in excess of N uptake by plants (for details refer to Ehrhardt et al. (2020) and   Table S1). In panel (d) boxplots represent annual N surplus averages of all catchments, gray lines represent mean values over all French and German study catchments, red lines represent single catchment time series.

Figure 2.
Conceptual framework of hysteresis between annual low-flow (C LF ) and high-flow (C HF ) concentrations. Upper box: Time series of the Wupper catchment for (a) observed daily discharge and weekly to monthly NO 3 -N concentrations, (b) observed and daily WRTDS-interpolated NO 3 -N concentrations, and (c) daily flow-normalized NO 3 -N concentrations with annual averages during low-flow (C LF ) and high-flow (C HF ) periods. Lower box: Synthetic (d) time series of annual C LF and C HF , and the seasonal logarithmic ratio log(C HF /C LF ), and (e) corresponding long-term hysteresis loops. The three distinct examples of long-term C-Q trajectories represent from left to right synchronous responses of C LF and C HF , asynchronous responses with C HF preceding and asynchronous responses with C LF preceding. The trajectories with a seasonally varying response times follow a hysteresis loop, while seasonally synchronous responses show no hysteretic behavior. Note that Figure 3 shows trajectories from three study catchments, including the Wupper catchment (as panels a-c), and Figure S1 presents more examples of possible trajectories. maximum (sprMax, April-July; metric seasQmax, see Table 1). The LF season was calculated but not used for an additional classification of average Q seasonality. Similarly, we determined the average seasonal concentration pattern at each station. We calculated the maximum monthly averages of C FN to classify the catchments into the ones with their intra-annual maximum during winter (winMax, November-March) and summer (sumMax, April-October; metric seasCmax, see Table 1). The combination of both classifications yields the long-term average C-Q seasonality.

Hysteresis in Multi-Annual Time Series of Low-and High-Flow Concentrations
We calculated mean C FN during LF (C LF ) and HF (C HF ) periods (see Section 3.2) for each hydrological year (November-October) and catchment from the WRTDS-interpolated daily C FN time series (Section 2.1, see Figure 2c as example). We only calculated these values with at least 80% daily data coverage over the corresponding time period. As a metric of the strength of C-Q seasonality within each year, we calculated the annual seasonal ratio in logarithmic form log(C HF /C LF ). Example trajectories of C LF , C HF and log(C HF /C LF ) are shown in Figure 2d.
To identify and characterize the seasonally asynchronous (i.e., phase/time shift between peaks) responses to changes in inputs, we plotted the time series of annual mean C LF and C HF against each other and applied a hysteresis approach (Figure 2e). Accordingly, the hysteresis captures decadal trajectories of the seasonality represented by time series of annual average low-and high-flow concentrations. Systematic changes in the trajectories are considered to be linked to changes in point and diffuse source inputs with mainly increasing inputs before 1990s and decreases mainly after 1990s (Figure 1d) in response to European regulations as drivers (EEC, 1991a(EEC, , 1991b. If C HF (y-axis) responds faster than C LF (x-axis) hysteresis becomes clockwise (CW), while if C LF responds faster than C HF hysteresis will be counter-clockwise (CC). If C HF and C LF respond simultaneously, no significant hysteresis loop is observed (n.s.). We used a hysteresis method based on the derivative of the X variable (dX), i.e. here C LF (Krueger et al., 2009;Minaudo et al., 2019). Its advantage for long-term water quality trajectories is that it also works for incomplete hysteresis loops. Incomplete loops are common for example, if (a) sampled time series are too short to cover the complete loop or (b) catchments do not return to their original state due to continued differences in the pressures or changes in their functioning. To this end, we fit the following nonlinear model using the R package minpack.lm (version 1.2-1).
Timing of NO 3 -N seasonal maximum of long-term months average Period of maximum monthly NO 3 -N C FN Classes: winMax, sumMax C HF , C LF Mean flow-normalized concentrations during HF and LF periods, calculated annually mg NO 3 -N l −1 log(C HF /C LF ) Logarithmic seasonal ratio, calculated annually log(C HF /C LF ) -rCseas Mean relative NO 3 -N seasonal amplitude between HF and LF seasons mean(abs(C HF -C LF ))/mean(C FN ) c Hysteresis between C LF and C HF See Equation 1, Classes:

Table 1 Metrics Used to Characterize Patterns of Q and NO 3 -N Seasonality
where t represents the time in years, d is the intercept, m is the slope between C HF ' and C LF ', and c is the hysteresis coefficient describing the direction and width of the hysteresis (c < 0 corresponds to a CC, c > 0 to CW, and c ≈ 0 to a non-significant hysteresis based on a confidence level of 95%). Prior to fitting, we linearly scaled C LF and C HF to vary between 0 and 1 (C LF ' and C HF ') to make the hysteresis coefficient c of different catchments comparable independent of differences in concentration (Lloyd et al., 2016).
Variations in seasonal response times cause trends in concentration seasonality and export dynamics, that is here log(C HF /C LF ) ( Figure 2a). Therefore, we use the Mann-Kendall test (R package Kendall, version 2.2) and Sen's slope (Sen, 1968) (R package trend, version 1.1.2) to detect and quantify monotonic trends in the logarithmic seasonal ratio. Missing annual values at 40 stations were filled by linear interpolation (on average 4.25 interpolated values, maximum eight), for example caused by large gaps in NO 3 − N samples (see Section 2.1).The ratio C HF /C LF is linked to the slope b of the logC-logQ regression and is thus another metric of C-Q relationships. The slope b is a widely used metric to characterize solute and particulate export dynamics in terms of C-Q relationships (Godsey et al., 2009) with a being the intercept and b the slope of the logC-logQ regression. If we assume a mean Q for HF (Q HF ) and another for LF (Q LF ) periods, the quotient becomes constant and b dependent on the logarithmic seasonal C FN ratio.
If the ratio of Q HF /Q LF equals Euler's number (∼2.718), b becomes equal to log(C HF /C LF ). Generally speaking, b behaves proportional to the logarithmic seasonal C FN ratio normalized by the logarithmic seasonal ratio of discharge log(Q HF /Q LF ). Though of course, slope b would be more scattered when considering the actual inter-annual Q variability instead of the long-term mean.
Additionally, we assessed and compared both the spatial variability and the long-term temporal variability of NO 3 -N seasonality. The spatial variability was determined for each year as the standard deviation of annual log(C HF /C LF ) across all stations. The temporal variability was calculated for each station as the standard deviation of annual log(C HF /C LF ) over the complete time series including the long-term trends in the seasonality trajectories.

Linking Metrics and Catchment Characteristics
We tested for relationships between the two classifications of long-term average C-Q seasonality and its trajectory (hysteresis) and catchment descriptors to examine our hypothesis about dominant controls. For the descriptors we used the data sets published and described in , Ebeling (2021) and extended to include French catchments in Ebeling and Dupas (2021). The descriptors characterize the hydroclimatic settings, topography, land cover, soil properties, geology, and point and diffuse N sources of the study catchments ( Figure S2, Table S1).
The analyzed metrics are partly categorical and partly continuous (Table 1). To test for significant differences among the different classes of long-term average C-Q seasonality and hysteresis trajectories, we used the non-parametric Wilcoxon test for comparing two classes and otherwise the Kruskal-Wallis test. A random forest (RF) classification model was trained for the two significant hysteresis classes (CC and CW) to investigate the joint predictive power of the descriptors while accounting for collinearities (see Figure S3). We used three times repeated tenfold cross-validation to estimate the mean model performance. Subsequently, permutation allows quantifying variable importance and thus identifying the dominant descriptors from the trained RF model. We used the importance permutation by Altmann et al. (2010) to identify significant descriptors from a first RF model using all descriptors (significance level 5%; Table S1), which then served to train a second RF model. We used the R package mlr3 (version 0.9.0, Lang et al., 2019) to train the RF model and ranger (version 0.12.1, Wright & Ziegler, 2017) for permutation variable importance. Spearman rank correlations served to identify relevant descriptors for the continuous metrics of average Q seasonality (rQseas) and NO 3 -N seasonality (mean log(C HF /C LF )).

Spatial Patterns of Long-Term Average Discharge and Nitrate Seasonality
The highest specific Q was observed in mountainous catchments (e.g., the Alps) with values up to 1,670 mm y −1 (median 340 mm y −1 among the study catchments) while lowest specific Q was 63 mm y −1 . The relative seasonal Q amplitude (rQseas) was highest in northwestern France (max 2.50) and lowest mostly in southeastern Germany (min 0.11) with a median of 1.18 (see Figure S4). A summary of calculated metrics across the study catchments is presented in Table 2.
Across the catchments, we found a dominance of winter maxima in Q (90.3%) and in flow-normalized NO 3 -N concentrations (91.0%) ( Figure 3) with occurrence of both in 83.8% of the catchments. The 9.7% of catchments with spring and early summer maxima of Q (sprMax) were apparent in mountainous areas ( Figure 3a). The 9.0% of catchments with summer maxima of NO 3 -N (sumMax) were observed mainly in northwestern and southern France (Figure 3a). In southern France, these summer NO 3 -N maxima coincided with Q minima, while NO 3 -N minima coincided with the Q maxima in spring. Catchments with spring maxima of Q in southeastern Germany have lowest Q during autumn or winter, which coincided with maxima in NO 3 -N. Figure S5 shows these synchronous and asynchronous long-term average seasonal C and Q variations of the four combined C-Q seasonality classes.
Catchments with the largest positive average seasonal ratios (log(C HF /C LF )) were found in western France, northern Germany and northeastern France (Figure 3a). On the other hand, catchments with average log(C HF /C LF ) < 0 were observed mainly in southern Germany, southern and western France. The spatial variability of the NO 3 -N seasonality among the catchments (standard deviations of log(C HF /C LF ) over all stations for each year) ranged from 0.42 to 0.64 with a median of 0.47 and an interquartile range of 0.04.

Long-Term Trajectories of Nitrate Seasonality
We examplarily visualized the seasonal NO 3 -N trajectories of three representative catchments namely the Wupper as an urban and the Leine and the Rance as agricultural rural catchments (Figures 1d and 3c). The Wupper catchment (605.9 km 2 ) has the highest fraction of artificial surfaces (31%), population density, and load by point sources and ranks under the top three for the fraction of point source loads from the total N input loads (86%) among the studied catchments. The HF season was determined for winter (December-February), LF for summer (May-July). The C LF were higher than the C HF in the beginning, but with C LF decreased from 1994 until C LF < C HF from 1998, indicating a switch in the seasonal timing. Note. Q-discharge, FN-flow normalized, LF-low flow, HF-high flow. The Q values used for these statistics was the same time period used for the WRTDS model, that is the overlapping period of NO 3 -N and Q time series.

Table 2 Summary Statistics of Calculated Metrics of All Catchments (n = 290) From Observations and the WRTDS Models: Median and Interquartile Range (IQR) and the Number of Catchments (n) With Positive (+), Negative (−) or Non-Significant (n.s.) Monotonic Mann-Kendall (MK) Trend
in C LF is two years before the peak in C HF , which creates a counterclockwise hysteresis. Accordingly, the seasonal ratio log(C HF /C LF ) follows a positive trend from a seasonal dilution (log(C HF /C LF ) < 0) to enrichment (log(C HF /C LF ) > 0) pattern. On the other hand, the two catchments with clockwise hysteresis, Leine (317.4 km 2 ) and Rance (142.6 km 2 ), are selected from catchments with high agricultural impact (65%-90%), strong change in diffuse agricultural N input from before to after 1990 ( Figure 1d, and small urban areas (2.6%-6.8%) as well as low fraction of point sources from total N input (below 10%)). The HF season in both catchments is in winter (Leine: January-March; Rance: December-February) and LF season in summer to fall (Leine: August-October; Rance: July-September). In the Leine catchment, C HF peaks in 1997 four years before C LF . In the Rance, C HF peaks in 2004 also four years before C LF with C LF being slightly higher than C HF . Both catchments have a negative trend in the seasonal ratio log(C HF /C LF ), which for the Leine means a trend from enrichment to a more neutral export pattern and a decrease in seasonal amplitude and for the Rance from neutral to dilution with an increase in the absolute seasonal amplitude between C LF and C HF .
Although the input trajectories of N surplus are relatively synchronous among the catchments (including the three example catchments) following the EU and national legislation (Figure 1d), the seasonal NO 3 -N trajectories observed in the streams differ more prominently among the catchments ( Figure S6). Exceptions to this generally synchronous variation in N surplus can be found in East German catchments with a drastic drop in 1990 after the German reunification, including the upper Leine catchment.
For all study catchments, counterclockwise and clockwise hysteresis occurred both in similar proportions (n = 87 CC, n = 85 CW, about 30% each), while the nonsignificant hysteresis class occurred slightly more often (n = 116 n.s., 40%). The fitted hysteresis models had a median performance of R 2 = 0.84 and a R 2 > 0.6 in 73.8% of the catchments. A spatial organization of the hysteresis patterns was apparent with clusters being more pronounced in France (Figure 3b): counterclockwise hysteresis dominated in central, southern and southwestern France and northern Germany, whereas clockwise hysteresis dominated in Brittany and northeastern France including the Seine catchment and southeastern Germany. In the majority of the catchments (80.3%, n = 233) the timing of seasonal C FN maxima and minima remained constant throughout the time series, whereas in 19.7% (n = 57) it switched between HF and LF periods (e.g., in the Wupper catchment, Figure 3c). The majority of the switching catchments (54%, n = 32) changed once from a seasonality with C LF > C HF to C LF < C HF with the median switch in the year 1991. Most of these catchments had counterclockwise hysteresis patterns (63%, n = 20). Some of those catchments (32%, n = 18) changed the seasonal timing of concentration maxima more than once.
The C LF -C HF -hysteresis patterns were connected to long-term trends in the annual seasonal ratios log(C HF / C LF ) (significance level 5%, Figure S6): 60.9% of the counterclockwise catchments had a positive trend in log(C HF /C LF ) and 49.4% of clockwise catchments had a negative trend in log(C HF /C LF ). The nonsignificant hysteresis was mostly linked to a significant increase in the seasonal ratio (43.0%). The different trajectories imply a long-term temporal variability in the NO 3 -N seasonality. The temporal variability of log(C HF /C LF ) at each station varied between 0.01 and 0.54 with a median of 0.09 and an interquartile range of 0.09.

Linking Seasonality Metrics and Catchment Characteristics
First, we linked the classifications of long-term average C-Q seasonality and the long-term trajectories to each other: The dominant class with NO 3 -N and Q winter maxima had similar proportions of all three hysteresis classes (25.9% for CW, 32.1% for CC and 42.0% for n.s.). The two asynchronous NO 3 -N and Q classes with Q winter and NO 3 -N summer maxima or Q spring/summer and NO 3 -N winter maxima both had higher proportions of clockwise hysteresis (63.2% and 47.6% respectively) compared to counterclockwise (10.5% and 14.3%) hysteresis. On the other hand, the asynchronous C-Q seasonality class in southern France with Q spring/summer and NO 3 -N summer maxima (LF period) had mostly counterclockwise (57.1%) and no clockwise (0%) hysteresis.
The timing of Q seasonality linked to elevation ( Figure 4a) and related hydroclimatic descriptors, that is catchments with Q spring/summer maxima had significantly higher mean elevations, precipitation and specific discharge and lower mean temperatures. The relative Q amplitude (rQseas) also correlated with the hydroclimatic drivers, especially the summer to winter precipitation P_SIsw (r = −0.63, Spearman rank, Figure S3) and mean temperature (r = 0.63).
The combined timing of NO 3 -N and Q seasonal maxima (Figure 4a) linked to more characteristics in addition to the controls of only Q seasonality. Compared to catchments with maximum NO 3 -N concentrations in summer, catchments with winter maxima of NO 3 -N tended to be more abundant in catchments with higher ratios of summer to winter precipitation (P_SIsw) and discharge (Q_seasR, not shown) and lower Q seasonality (rQseas), although the class of concurrent NO 3 -N and Q winter maxima contained many catchments with equilibrated summer to winter precipitation (P_SIsw∼1) and smaller precipitation seasonality (P_SI, not shown). Summer maxima in NO 3 -N emerged in catchments with consistently lower summer to winter precipitation (P_SIsw < 1, with only one exception) and lower depths to bedrock (dtb). We also found a link of average NO 3 -N concentration to long-term average C-Q seasonality: The mean C FN was lowest for the classes with Q spring/summer maxima and highest for Q winter maxima with NO 3 -N summer maxima (especially in the Armorican Massif). The major class of concurrent NO 3 -N and Q winter maxima has medium mean C FN values but also a high range. . As expected from the methods, the classes were clearly linked to the strength of the NO 3 -N-Q seasonality, quantified by the mean seasonal ratio log(C HF /C LF ), with the class of concurrent NO 3 -N and Q winter maxima showing mostly positive and the other classes negative values. The mean ratio log(C HF /C LF ) was correlated to topographical and hydroclimatic descriptors (slo.mean r = −0.45, twi.mean r = 0.43, AI r = 0.41, specific Q r = −0.39) and land use (f_agric r = 0.39). Note that the descriptors also correlate with each other ( Figure S3).
The RF model for the two significant hysteresis classes reached only a mean accuracy of 0.69 ± 0.12 from cross-validation (compared to 0.5 for random guess). The descriptor relative Q seasonality (rQseas) ranked highest in the feature importance, followed by evapotranspiration (PET_mm), mean C FN , precipitation seasonality (P_SI), ratio of summer to winter discharge (Q_seasR) and fraction of sand (f_sand), although they showed overall relatively low variable importances ( Figure S7). Expected links to descriptors of N loading were not found dominant: although N surplus descriptors were significant in the RF model, they only ranked 8th or lower. As shown in Figure 4b, the clockwise hysteresis was mostly noticed in catchments with higher Q and precipitation seasonality (rQseas and P_SI), and mean C FN and lower evapotranspiration (PET_mm), fraction of sand (f_sand), and temporal variability in NO 3 -N seasonality (standard deviation of log(C HF /C LF )). The clockwise catchments covered a large variability in the fraction of agriculture, N surplus, the difference in N surplus before and after 1990 (dNsurp71_91), the fraction of point source loads (N_WW_ frac) and artificial surfaces (f_artif). The nonsignificant hysteresis class contained more catchments with higher fractions of point source loads, albeit with a large variability of values.

Spatial Patterns and Controls of Long-Term Average Nitrate and Discharge Seasonality
The long-term average NO 3 -N seasonality was linked to average seasonality of Q in combination with climate, topography and geology. Accordingly, we established distinct archetypal patterns across the large range of catchments and provide possible links to the underlying catchment characteristics. The different timing of Q seasonality (seasQmax) across the catchments was controlled by topography and related hydroclimatic variables. This agrees with Kuentz et al. (2017), who found climate was the main driver for several flow characteristics at the European scale, and Gnann et al. (2020), who found catchment aridity to control the timing and magnitude of Q seasonality relative to climate seasonality. The dominance of Q winter maxima (90.3%) and summer low flow ( Figure S5) suggests that the seasonal cycle of evapotranspiration strongly affects the Q cycle in the study region. Precipitation seasonality (P_SIsw) with relatively high winter precipitation (as prevalent in France) can add to this pattern, especially in northwestern France with the lowest ratios of summer to winter precipitation (P_SIsw) and relatively low annual evapotranspiration and aridity. In that region, the opposing seasonal precipitation and evapotranspiration cycles in a moist climate lead to strong seasonal Q patterns (rQseas) and short hydrological transport times (Gnann et al., 2020;Kumar et al., 2020), which might also be enhanced by shallow aquifers with shallow depths to bedrock (dtb). The Q maximum is shifted to spring or summer in fewer catchments (9.7%), where winter precipitation is likely retained as snow-pack and discharged with higher temperatures during thawing period. Here, this effect seems to exceed the evapotranspiration control. The spatial patterns of high-flow seasonality (Figures 3a  and S4) largely agree with the flood timing reported across the study region, especially the geographical differences between the mountainous, snow-impacted catchments (in the Alps and Pyrenees) and others (Berghuijs et al., 2019;Blöschl et al., 2017).
In combination with the NO 3 -N seasonality, we found three archetypal patterns of long-term average C-Q seasonality with different controls and underlying processes: 1. Both NO 3 -N and Q have a winter maximum (83.8% of the catchments) and vary mostly synchronously (Figures 3a and S5), which indicates prevalent enrichment patterns (transport limitation) and NO 3 -N mobilization processes at a seasonal scale . This supports our first hypothesis that overall synchronous NO 3 -N and Q seasonality dominate. The dominance of NO 3 -N and Q winter maxima and C-Q enrichment patterns agrees with previous research on temperate catchments Minaudo et al., 2015;Moatar et al., 2017;Musolff et al., 2015;Zhi & Li, 2020). Moreover, the spatial patterns of the strength of nitrate-discharge seasonality, represented by the average seasonal ratio log(C HF /C LF ), are also in line with a previous study of Germany-wide C-Q relationships showing strong enrichment patterns in northern Germany and less pronounced enrichment or neutral patterns in central-eastern Germany . This supports that the seasonal ratio is comparable to the widely used slope b for inter-annual C-Q relationships (see Equation 4). This is plausible as seasonal NO 3 -N variations exceed storm-event induced variations, when considering rather large catchments are usually monitored at low (typically monthly) frequency . Furthermore, our analysis showed the dominant control of topography, climate and land use in shaping the spatial variability of the seasonal ratio, with a tendency of higher ratios in agricultural lowland catchments. With dominant diffuse N sources, a mechanism leading to this synchronous archetype can be variable discharge generating zones, where shallow N sources are activated by younger water dominating during HF and lack of connectivity of sources along with longer reaction times (potentially higher removal) for longer flow paths dominates during LF (Benettin et al., 2020;Musolff et al., 2016;Yang et al., 2018;Zhi & Li, 2020). This interpretation is supported by , who additionally found a dominant control of the vertical concentration gradient on the NO 3 -N export dynamics in German catchments. Reduced riparian or instream NO 3 -N removal during winter following the seasonal biogeochemical cycle  could also enhance the pattern of synchronous C and Q seasonality. 2. The NO 3 -N maximum coinciding with the Q minimum and vice versa (9.7% of the catchments, Figures 3a and S5) occurred in mountainous catchments, where asynchronous dilution patterns of C-Q dominate. Spring and summer discharge from mountains stemming from snow melt and summer precipitation could dilute N sources. In these archetypal catchments, mean NO 3 -N concentrations are low (Figures 1c and 4a) due to relatively low N input in combination with strong dilution by high specific discharge (Figure 1b). The dominant dilution pattern agrees with C-Q relationships observed previously in mountainous, snow-dominated German meso-scale catchments . The dilution of N sources is lowest during LF season, that is late summer in southern France and autumn or winter in southern Germany. The dilution pattern of this archetype likely results from the spatial separation of discharge generating zones upstream and dominant agricultural sources downstream, potentially masking other export dynamics in the downstream areas. 3. NO 3 -N maxima are reached during summer LF and minima mostly during autumn before the winter HF period (6.6% of the catchments, Figures 3a and S5). This type of NO 3 -N seasonality with relatively weak dilution patterns, high mean concentrations and small relative seasonal differences between LF and HF concentration seems to be specific to the Armorican Massif and has been described for several catchments in this region previously Guillemot et al., 2021;Martin et al., 2004). The strong Q seasonality in that region and fast hydrological responses (e.g., high flashiness and CVQ) are not reflected in the NO 3 -N seasonality, indicating a rather chemostatic export regime (low concentration relative to high Q variability) for this archetype. Reasons for this pattern include large legacy stores  and bottom-loaded profiles with NO 3 − -rich groundwater which is diluted during HF (Martin et al., 2006) in combination with potential bypassing of the riparian zones where removal by denitrification can be high during LF . However, riparian zones can also act as NO 3 − sources during LF (Duncan et al., 2015), offering a nonexclusive alternative mechanism for this pattern. Van Meter et al. (2020) found that NO 3 − -rich groundwater can lead to an out-of-phase or aseasonal NO 3 − pattern. In case of bottom-loaded profiles, we would expect the catchments to have developed the reverse seasonality over time, as described in for example Ehrhardt et al. (2019), which was indeed observed for some catchments (see Figure S6, crossing the 1:1 line). Alternatively, multiple contributing aquifers with variable response times and sources could be horizontally distributed differently with a longitudinal gradient, with a dominance of nonagricultural wetland soils upstream and agricultural areas downstream. With the upstream discharge diluting downstream sources during HF, this process resembles archetype 2. Another contributing mechanism to this pattern could be the seasonal variation of N inputs as fertilizer is primarily applied in spring or of retention in the soils.

Long-Term Trajectories of Nitrate Seasonality Integrate Complex Controls
We observed both long-term trajectory types equally often, with annual C LF (counterclockwise hysteresis) or C HF (clockwise hysteresis) changing first. Our hypothesis of dominant clockwise hysteresis for NO 3 − is therefore not supported by the studied catchments. Surprisingly, 20% of the catchments even switched their seasonal timing at least once during the observation period, which imposes substantial changes to aquatic ecosystems over time. As described in the introduction, this could create shifts in nutrient availability and limitation during sensitive time periods for aquatic organisms and ecosystems (Minaudo et al., 2015).
The C HF and C LF are influenced by various controls that may lead to the observed long-term trajectories and hysteresis patterns. For example, in the densely populated Wupper catchment, changes in point sources could have affected C LF more and earlier than C HF , as discussed by . Reductions in inputs from point sources usually have a more immediate effect on riverine concentrations and have typically been achieved more easily and earlier than reductions in diffuse-source inputs in Western Europe (Le Moal et al., 2019;Westphal et al., 2020). The example catchments Leine and Rance (Figure 3c) are agriculturally managed with few point sources only and show the expected stronger delay in low-flow concentrations (clockwise hysteresis). This suggests that our hypothesis of changes in diffuse sources being reflected first in changes in C HF may be supported for these two catchments but cannot be verified across the entire set of catchments.
Despite the large sample size and rich ancillary data on catchment characteristics, we did not detect fundamental or general controls on hysteresis patterns. Over the whole range of study catchments, the source descriptors (i.e., land cover, N surplus, point sources, population density) did not show consistent relationships with hysteresis. Small differences between the hysteresis classes were apparent in catchment properties describing the hydroclimate and soil data, though they had a relatively low explanatory power. This suggests overall complex controls on the long-term hysteresis patterns over the large extent of the study, with several interacting controls potentially determining the observed trajectories that will be addressed in the following.
The complexity in controls of responses across a wider range of catchments is in line with recent research on the diversity of nutrient retention capacity in various surface and subsurface catchment components (Frei et al., 2020;Kolbe et al., 2019). In particular, the high hydrochemical complexity of the subsurface in combination with relatively sparse subsurface data availability makes regional to continental predictions exceedingly uncertain Jawitz et al., 2020;Marçais et al., 2018). Moreover, the subsurface reactivity can experience long-term changes if availability of electron donors and thus denitrification potential decreases (Bouwman et al., 2013). Because soil, vadose, and aquifer hydrology and biogeochemistry directly influence both hydrological time lags and active nutrient retention, better characterization of subsurface parameters should be a major research priority (Condon et al., 2020;Li et al., 2021).
The diversity of human activities can also obscure the linkage between catchment characteristics and nutrient dynamics. Catchments with anthropogenic impact often have both point and diffuse N sources, creating complex and time-variant interactions that are not captured by the catchment attributes taken into account here. For example, in diffuse-source-dominated catchments, smaller point sources could still have changed more than diffuse sources and vice versa, hampering predictability of hysteresis by catchment descriptors representing more the source strengths. Point sources could mask changes caused by diffuse source changes especially in larger, more populated catchments (Dupas et al., 2017).
Another unaccounted factor could be long-term changes in in-stream nutrient uptake processes. For example, an increasing C LF while C HF is already decreasing could result from reduced uptake rates during summer. Reduced algal biomass production could be caused by reductions in phosphorus and increasing limitations (Bowes et al., 2011;Minaudo et al., 2015). For example, in the Frome catchment, changes in instream processes were even considered to dominantly control changes in NO 3 − seasonal amplitude because concurrent seasonal changes of NO 3 − in groundwater boreholes were not observed (Bowes et al., 2011). Additionally, NO 3 − could be released to the stream by decomposition of biomass in winter and further affect seasonality (Bowes et al., 2011).
The tendency of clockwise catchments to have higher precipitation seasonality (P_SI, Figure 4b) and low evapotranspiration (PET_mm) could relate to their occurance in the Alps and the Armorican Massif (Figure 3b). Conversely, the high relative Q seasonality (rQseas) was common for catchments in northeastern France and the Armorican Massif. This might indicate that pronounced hydroclimatic variability in relatively wet catchments enhances HF concentrations to react faster, whereas high PET_mm links to dominant counterclockwise hysteresis in central and southern France, although this does not seem to be an overall dominant control as predictive power remained low.
Interestingly, we found distinct hysteresis patterns in the regions with asynchronous archetypal patterns of long-term average C-Q seasonality (16.2% of the catchments). For example, in southern France, the strong seasonal dilution pattern was mostly combined with counterclockwise hysteresis, whereas in the Alps and in the Armorican Massif, clockwise hysteresis dominated (Figure 3b). Dilution of downstream sources from upstream discharge generation in mountainous catchments could generally cause a stronger trend in less diluted low-flow concentration when sources change. However, these diluting effects interact with other mechanisms at downstream locations, for example mobilization processes. Possibly, in southern France, the upstream signal dominates, while in the Alps with smaller relative seasonal Q variations, downstream signals could dominate the overall trajectories, and leading to faster high-flow concentration responses. For the Armorican Massif, high-flow concentration responding first to changes in dominant diffuse N-input could cause the clockwise hysteresis patterns. This is supported by generally higher deep and lower shallow groundwater concentrations causing dilution patterns which could result from long-term N-inputs and large legacy stores . On the other hand, within the dominant archetype of synchronous average C-Q seasonality, the heterogeneity in trajectories and controls was high and no dominant general pattern and control was detectable in this study (as discussed above).
In essence, the seasonality trajectory and its corresponding hysteresis integrate all the above mentioned processes. The missing dominance emphasizes that over the study domain multiple controls are relevant and hierarchies vary. Many settings can lead to the same response pattern and seasonality trajectory.

Spatial Variability Larger and Better Predictable Than Long-Term Temporal Variability of Nitrate Seasonality
Spatial variability of seasonality among the catchments is larger than the long-term temporal variability of most catchments, even when including significant trends in the logarithmic seasonal ratio. Although the number of significant long-term trends in the seasonal ratio log(C HF /C LF ) (76%) and hysteresis patterns (59%) reveal distinct changes in NO 3 − seasonality on multi-decadal time scales in the majority of catchments, inter-catchment variability across the study area largely exceeded the temporal variability within a catchment. Dupas et al. (2019) and Abbott, Gruau, et al. (2018) also report examples of higher spatial than temporal variability, resulting in spatial persistence, that is stable spatial patterns among catchments through time. These two studies, based on shorter periods of six  and 12 years (Abbott, Gruau, et al., 2018), include mainly seasonal variations. Our finding thus extends the spatial persistence concept for NO 3 − seasonality to decadal trajectories across a wide range of Western European catchments.
At the same time, we can better explain the spatial variability of long-term average seasonality and C-Q dynamics among the catchments (Section 4.1). In contrast, the differences in trajectories resulting from different response times of low-flow and high-flow concentrations were only poorly predictable from the available data, and only small differences in catchment characteristics were apparent between the hysteresis classes (Section 4.2). This could be due to the fact that trajectories are affected by multiple controls with asynchronous source changes and seasonally variable responses, possibly masking each other. Further investigation may include temporarily variable catchment descriptors, such as long-term data on point source inputs, and time series of ecological in-stream metrics such as chlorophyll-a.

Conclusions
We characterized average long-term NO 3 − -Q seasonalities and their trajectories from long-term time series in 290 French and German catchments covering a large variety in hydroclimate, topography, lithology and anthropogenic pressures. We implemented a novel hysteresis approach using low-and high-flow NO 3 -N concentration trajectories. Our main findings are: • We observed a widespread dominance of concurrent maxima of Q and NO 3 -N in winter (84%), supporting our hypothesis. Deviations from this archetype of long-term average C-Q seasonality were linked to topography and hydroclimatic seasonality, and to source heterogeneity or lithology especially in the Armorican Massif • Surprisingly, counterclockwise and clockwise hysteresis patterns of low-flow and high-flow NO 3 -N concentrations occurred equally often, thus we had to reject our hypothesis about dominant occurrence of faster response during high flow. We exemplarily showed that indeed counterclockwise hysteresis can be found in point source dominated and clockwise hysteresis in agricultural catchments. However, the lack of a consistent pattern in the French and German catchments suggests a high level of complexity and potential interactions in the controls over the wide range of catchments. These controls need to be further disentangled in a set of catchments with better known input time series. We point out that scarce data especially on point source loads (including their temporal evolution) is limiting the understanding of catchment functioning across a large sample. Therefore monitoring efforts and data use policies should be improved • Overall, we found that spatial variability of long-term average C-Q seasonality between the catchments was larger and easier to predict than its long-term temporal evolution at a single station Our large sample study has several implications for water quality management. The dominant spatial variability indicates that (a) uniform regulations may not be appropriate but need to take the spatial variability into account and should target the management regionally, and (b) short-term monitoring can already be useful to characterize the overall system functioning. The distinction between low-and high-flow trajectories can guide ecological assessment and water quality management considering that aquatic ecosystems are more prone to eutrophication during low flow, while exported loads are more susceptible to changes in high flows.