the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Bending of the concentration discharge relationship can inform about instream nitrate removal
Fanny Sarrazin
Rohini Kumar
Jan H. Fleckenstein
Andreas Musolff
Nitrate (${{\mathrm{NO}}_{\mathrm{3}}}^{}$) excess in rivers harms aquatic ecosystems and can induce detrimental algae growths in coastal areas. Riverine ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ uptake is a crucial element of the catchmentscale nitrogen balance and can be measured at small spatiotemporal scales, while at the scale of entire river networks, uptake measurements are rarely available. Concurrent, lowfrequency ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ concentration and streamflow (Q) observations at a basin outlet, however, are commonly monitored and can be analyzed in terms of concentration discharge (C–Q) relationships. Previous studies suggest that steeper positive log (C)–log (Q) slopes under low flow conditions (than under high flows) are linked to biological ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ uptake, creating a bent rather than linear log (C)–log (Q) relationship. Here we explore if networkscale ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ uptake creates bent log (C)–log (Q) relationships and when in turn uptake can be quantified from observed lowfrequency C–Q data. To this end we apply a parsimonious massbalancebased river network uptake model in 13 mesoscale German catchments (21–1450 km^{2}) and explore the linkages between log (C)–log (Q) bending and different model parameter combinations. The modeling results show that uptake and transport in the river network can create bent log (C)–log (Q) relationships at the basin outlet from log–log linear C–Q relationships describing the ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ landtostream transfer. We find that within the chosen parameter range the bending is mainly shaped by geomorphological parameters that control the channel reactive surface area rather than by the biological uptake velocity itself. Further we show that in this exploratory modeling environment, bending is positively correlated to percentage of ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ load removed in the network (L_{r.perc}) but that networkwide flow velocities should be taken into account when interpreting log (C)–log (Q) bending. Classification trees, finally, can successfully predict classes of low (∼4 %), intermediate (∼32 %) and high (∼68 %) L_{r.perc} using information on water velocity and log (C)–log (Q) bending. These results can help to identify stream networks that efficiently attenuate ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ loads based on lowfrequency ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ and Q observations and generally show the importance of the channel geomorphology on the emerging log (C)–log (Q) bending at network scales.
 Article
(6044 KB)  Fulltext XML

Supplement
(642 KB)  BibTeX
 EndNote
Transport and transformation of nitrate (${{\mathrm{NO}}_{\mathrm{3}}}^{}$) in river networks are major controls of downstream exports to receiving lakes, reservoirs and coastal systems (Alexander et al., 2000; Billen et al., 1991; Peterson et al., 2001; Seitzinger et al., 2002; Seybold and McGlynn, 2018). Increased ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ concentrations in surface waters can induce detrimental algae growths (Beusen et al., 2016; Canfield et al., 2010; Galloway et al., 1995), compromise river ecosystem health and jeopardize drinking water supplies. Since the beginning of the 20th century, human activities such as agricultural expansion and fossil fuel burning have mobilized additional reactive nitrogen (N), initiating and later exacerbating this problem (Seitzinger et al., 2002). In arable landscapes, which include large parts of Europe, the efficient management of aquatic ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ at network scales is complicated by the spatiotemporal variability in loading patterns and hydrologic regimes as well as the lack of understanding of nutrient pathways, connected transit times and removal processes from input to export. Nevertheless, nitrate concentration and load variability can be predicted at catchment scales when relying on detailed process understanding regarding transport and biogeochemical processing (Alexander et al., 2009; Schlesinger et al., 2006; Wollheim et al., 2008). Moving beyond smallscale variability and characterizing nitrate processing at the catchment scale however remains a challenge (McDonnell et al., 2007; Li et al., 2020).
Within river reaches and streams, reactive solutes like ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ are affected by complex interactions of physical, biological and chemical processes. Physical transport is driven by local discharge and channel geomorphology and dictates the ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ residence time in a reach, thus influencing the timescales at which biogeochemical processing can take place (Kirchner et al., 2000; Runkel and Bencala, 1995; Zarnetske et al., 2011). ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ is removed and transformed by denitrifying bacteria in the anoxic river sediment (Birgand et al., 2007; Peterson et al., 2001); ammonified; or retained through assimilation processes in the oxic or anoxic river compartments by bacteria, fungi and primary producers such as algae and macrophytes, potentially entering higher trophic levels. In the latter case, N in the form of DON (dissolved organic nitrogen) and more commonly DIN (dissolved inorganic nitrogen), together with phosphorus (P), may be released to the water column later on (Durand et al., 2011; Vanni, 2002; Vanni and McIntyre, 2016). The nutrient spiraling model (Newbold et al., 1981; Stream Solute Workshop, 1990) that formally describes these processes has been used to quantify and compare ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ transport and uptake (the net result of all removal and release processes) in river reaches (Peterson et al., 2001; Mulholland et al., 2008; Hall et al., 2009) and stream networks (Ensign et al., 2006; Doyle, 2005; Marce and Armengol, 2009). Quantifying in situ ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ uptake is laborintensive and may involve local nutrient additions, potentially altering the ambient uptake rate (Hensley et al., 2014; Mulholland and Tank, 2002). Other methods require highfrequency measurements (Jarvie et al., 2018; Kunz et al., 2017) that are mostly limited to small spatial scales (i.e., reach scale) and can vary considerably between measuring points (Boyer et al., 2006). At the scale of entire river networks contrarily, uptake measurements are rarely available (but see Wollheim et al., 2017), and models are applied instead to predict spatiotemporal uptake patterns (Boyer et al., 2006; Yang et al., 2018). These models account for the spatial configuration of the stream network, an important aspect for stream biogeochemistry that is often ignored in smallscale experimental approaches (Fisher et al., 2004). Spatially distributed models, however, require calibration of uncertain spatiotemporal parameters and may not reflect the essential features of the system despite fitting observed data well (Klemes, 1986).
River networks link terrestrial source zones to coastal areas and integrate biogeochemical and hydrological catchment functions across scales (Bouwman et al., 2013; Helton et al., 2018). Small streams (usually headwaters) are known to influence the export signal disproportionally because of their overall (high) contribution to total stream length and effective ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ removal capacity (Alexander et al., 2000; Horton, 1945), explained by high sedimentsurfacetowatervolume ratios. Generally, high removal efficiencies have been reported for river network areas with lower specific discharges (Hall et al., 2009; Hensley et al., 2014) under favorable circumstances such as high light availability and heavy instream vegetation (Hensley et al., 2014; Rode et al., 2016) as well as for streams with a high capacity for lateral and hyporheic exchange (GomezVelez et al., 2015; Kiel and Cardenas, 2014). The scaling of instream uptake processes beyond the river reach has been approached by combined experimentalmodeling studies with defined explicit scaling relationships (e.g., Basu et al., 2011; Aguilera et al., 2013; Bertuzzo et al., 2017; Lindgren and Destouni, 2004) and theoretical frameworks explaining how the river network capacity regulates solute export (e.g., Wollheim et al., 2018). Abbott et al. (2018) shows how spatially heterogeneous patterns of water chemistry stabilize while the temporal variability in nutrient concentrations persists when moving downstream, facilitating the temporal scaling of headwater measurements. Nevertheless, insights into linking the interplay of nitrate removal processes at the network scale to downstream export patterns in space and time are largely missing.
Concentrations (C) for instream solutes such as carbon, major ions, particulates and nutrients are commonly monitored concurrently with discharge (Q) at the basin outlet. C–Q relationships integrate the effect of biogeochemical and hydrological processes within the catchment and have mainly been discussed in terms of land–stream transfer and source configuration in catchments as well as subsurface retention processes (Godsey et al., 2009; Musolff et al., 2017; Bieroza et al., 2018). The shape of longterm (multiple years) C–Q relationships in the log–log space is typically described by the slope of a linear regression model (Godsey et al., 2009). Here, three archetypes have been distinguished: (i) a positive log (C)–log (Q) slope, indicating enrichment, occurs when an increasing discharge additionally mobilizes solutes; (ii) a negative C–Q slope or dilution pattern is commonly linked with source limitations; and (iii) a neutral or chemostatic slope implies low variability in instream concentrations across a high range of discharges, a pattern observed for example for solutes derived from weathering bedrock (Ameli et al., 2017; Godsey et al., 2009). The potential information loss associated with linear and monotonic ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ log (C)–log (Q) analysis was addressed by Moatar et al. (2017) and Minaudo et al. (2019) for more than 200 French catchments and by Diamond and Cohen (2018) for 44 rivers in Florida, USA. These studies identified distinct linear lowflow and highflow ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ log (C)–log (Q) regression slopes for a majority of the cases using lowfrequency monitoring data. Moatar et al. (2017) found that stronger positive slopes under lowflow conditions correlate positively with chlorophyll a concentrations (associated with biological processes) and attributed this condition to biological ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ concentration mediation in the stream. This is consistent with the findings of Hall et al. (2009) and Hensley et al. (2014) among others that instream uptake is more efficient under lowflow than under highflow conditions. Furthermore, Wollheim et al. (2017) illustrate nonlinear ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ C–Q relationships conceptually for storm flow dynamics in a river network, showing high retention capacities in the headwater catchments that decrease under increasing flows, changing the slope of C–Q relationships from dilution to enrichment. Based on these studies we hypothesize that the magnitude (or efficiency) of instream ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ uptake is encoded within observed C–Q relationships, and their analysis therefore can improve our understanding of instream uptake processes by providing an alternative to elaborate fieldwork and modeling work aimed at quantifying ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ removal in stream networks. Lowfrequency ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ observations are widely available (e.g., biweekly to monthly grab sampling; Ebeling et al., 2020; Minaudo et al., 2019; Moatar et al., 2017), but if and how these data can be utilized to characterize catchmentscale instream processing has yet to be investigated.
In this paper, it is postulated that networkscale uptake effects can be inferred from the degree of nonlinearity or the amount of bending of lowfrequency, multiannual concentration (C) and discharge (Q) observations. To test this hypothesis, we apply a parsimonious river network model (similar to Bertuzzo et al., 2017; Helton et al., 2018, 2010; and Mulholland et al., 2008) in 13 German catchments to explore the catchmentscale transport and uptake processes that influence downstream log (C)–log (Q) patterns. The specific objectives are to (i) introduce the maximum curvature (Curv_{max}) as a robust metric to quantify bending of lowfrequency C–Q time series in the log–log space, (ii) explore the sensitivity of Curv_{max} to hydrological and instream biogeochemical parameters (e.g., channel shape, water velocity and biological ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ uptake velocity), (iii) explore how C–Q bending is linked to networkscale instream uptake, and (iv) provide guidelines if and under what circumstances the C–Q bending can offer conclusive information on effective instream uptake. In this proofofconcept exploratory study, we demonstrate how (existing) lowfrequency monitoring data can be effectively utilized to quantify nitrate uptake in river networks and show how smallscale uptake processes shape emerging patterns at catchment scales.
2.1 Maximum curvature – Curv_{max}
The shape of log (C)–log (Q) relationships is often described as linear (Bieroza et al., 2018; Godsey et al., 2009; Musolff et al., 2017) or segmented linear (Meybeck and Moatar, 2012; Moatar et al., 2017; Marinos et al., 2020), implying a limit on the possible log (C)–log (Q) shapes (Tunqui Neira et al., 2020) and setting assumptions such as “fixed breaking points”. Here, we introduce the concept of maximum curvature (Curv_{max}) to quantify rather than describe the shape of “brokenstick” log (C)–log (Q) relationships without the assumption of a fixed form. In a strict geometrical sense the curvature (−∞; +∞) is the instantaneous rate of change in direction of a point that moves on a curve. A straight line for example has a curvature of zero, and a large circle has a lower absolute curvature than a small circle (Pressley, 2001). Here, Curv_{max} identifies the magnitude and direction of the log (C)–log (Q) section with the largest instantaneous change. To calculate Curv_{max} for an observed (noisy) log (C)–log (Q) relationship, a smoothed spline, f, is iteratively fitted with increasing degrees of freedom (df) to capture the general log (C)–log (Q) shape accurately but avoid overfitting (Fig. B1). Initially, df=3, and the log (Q) region of the largest instantaneous change is identified as Q_{m}±0.05, with ${Q}_{\mathrm{m}}=\mathrm{arg}{max}_{\mathrm{log}Q}\left{f}^{\prime \prime}\right$. Then, df is increased until, at df=i, the log (Q) corresponding to the largest instantaneous change is not within the initial Q_{m} region anymore. Consequently, Curv_{max} is calculated for a smoothed spline fit, f, with
The Curv_{max} metric, computed for a log (C)–log (Q) relationship, could be considered to be a complementary metric to the slope of the linear regression model (Godsey et al., 2009) and could serve as an alternative for segmented linear regression fits (Meybeck and Moatar, 2012; Moatar et al., 2017; Marinos et al., 2020) (Fig. S1 in the Supplement) as it quantifies the degree of nonlinearity as the amount of bending.
We assume here that a multiannual (6 to 15 years) lowfrequency (biweekly to monthly) C–Q relationship without temporal (significant) trends in a given station has one Curv_{max}. To verify this assumption in a realistic setting, Curv_{max} was computed for observed nitrate (C) and Q data (1995–2010) of French water quality stations with biweekly to monthly sampling frequencies (Dupas et al., 2019). Following the removal of C outliers (falling outside of μ±3.5σ in the log space, with μ and σ representing the sample mean and standard deviation, respectively), 444 stations were selected that satisfy the following four criteria: (i) the station should have at least 70 coupled C and Q observations (the maximum number of samples within one station was 402), (ii) a minimum of 6 years of data are represented, (iii) there is no bias in the intraannual distribution of the data (i.e., never less than 10 % of the C–Q observations in one season, Fig. S2 in the Supplement), and (iv) the station C observations had no significant temporal trends (Mann–Kendall test, p value>0.05) (Ebeling et al., 2021). We then assessed the robustness of Curv_{max} to the lowfrequency C–Q observations in a time series by selecting different subsamples of C–Q data from the entire available time series for a given station. More specifically, 100 random time series subsamples (each with a minimum length of 70) without replacement but with overlap were taken for each station, with the subsamples passing the four criteria above, and Curv_{max} was calculated for each subsampled time series. On average, the subsamples represented nearly 80 % of the complete time series for a station.
2.2 Network model
In this work an explorative gridbased (100 m×100 m) mass balance network model (comparable to Bertuzzo et al., 2017; Helton et al., 2018, 2010; and Mulholland et al., 2008; conceptually shown in Fig. B2) was used to simulate instream nitrate transport and biological removal on a daily basis. The model was developed in R (R Core Team, 2013).
2.2.1 Stream network and hydrological properties
Following Bertuzzo et al. (2017) and Helton et al. (2018), each river network node (i.e., grid cell) i ($\mathrm{1}\le i\le N$) has a drainage area A_{i} [m^{2}] that is calculated as the sum of the total upstream drainage area ∑_{j}W_{ji}A_{j} [m^{2}] and the direct drainage area a_{i} [m^{2}] (e.g., laterally contributing drainage area) to grid cell i (Eq. 1):
where W_{ji} [–] is an element in the connectivity matrix W (N×N) such that W_{ji}=1 if j is directly neighboring and flowing into i, and W_{ji}=0 otherwise. A_{j} [m^{2}] is the total drainage area to node j.
It is the spatially varying contribution of the total upstream and direct drainage area to each river network node that creates the spatial variability in the river network. The temporal variability in the river network is driven by the timevarying but spatially homogeneous specific discharge Q_{t.sp} [m^{3} s^{−1}], calculated as the ratio of the daily discharge at the catchment outlet and the total number of catchment grid cells. The total local discharge Q_{i} [m^{3} s^{−1}] at a given grid cell i is thus proportional to the total drainage area at that grid cell, A_{i}, (following Bergstrom et al., 2016, and Bertuzzo et al., 2017) (Eq. 2). It is Q_{i}, which in turn dictates the downstream and atastation hydraulic geometry relationships of river geomorphic parameters channel width, (w_{i}; m) and average channel depth (d_{i}; m) (Leopold and Maddock, 1953) (Eqs. 2a and b). The local velocity in a grid cell v_{i} [m s^{−1}] is calculated according to Eq. (2c), and the corresponding travel time, T_{i} [d], is computed in Eq. (2d):
where the flow length through a grid cell i, l_{i} [m], equals 100 or $\mathrm{100}\sqrt{\mathrm{2}}$ m for horizontal and vertical or diagonal flow directions, respectively. Parameters a_{w} [–] and K_{w} [–] are the respective exponent and coefficient parameters in the river width–discharge relationship (Eq. 2a), while a_{d} [–] and K_{d} [–] compose the exponent and coefficient parameters of the depth–discharge relationship (Eq. 2b), respectively. The ratio of a_{d} to a_{w} corresponds to a parameter $r\phantom{\rule{0.25em}{0ex}}\text{[\u2013]}\in {\mathrm{R}}^{+}$, which prescribes the crosssection geometry relation such that a triangular channel crosssection is represented by r=1, a parabolic channel crosssection by r=2, and channel crosssections with progressively flatter bottoms and steeper banks by increasing values of r (Dingman, 2007). The width–discharge relation in Eq. (2a) is conceptually illustrated in Fig. B3 for two sets of a_{w} and K_{w}, where a low a_{w} corresponds to the width of a channel that does not change much with varying discharge, while a high a_{w} can result in highly varying channel widths.
2.2.2 Nitrate uptake
Similar to Eq. (1) the incoming load, L_{in,i} [mg s^{−1}], to a river network grid cell i is composed as the sum of upstream load contributions L_{in.up,i} [mg s^{−1}] and direct landtostream loading L_{in.ls,i} [mg s^{−1}], given that L=CQ (Eq. 3). The contribution of direct landtostream loading concentration can be expressed as a power law (Musolff et al., 2017) with the exponent b [–], the slope in the log (C)–log (Q) relationship that is an indicator of the C–Q archetype (Godsey et al., 2009) and coefficient c [–]. Here, b is assumed to be constant over the seasons, which considers that ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ loading is transportlimited rather than sourcelimited as explicitly shown for agricultural catchments (Basu et al., 2011; Winter et al., 2021). Following Jawitz and Mitchell (2011), the coefficient c is calculated to yield the longterm mean instream input concentration C_{mean} [mg L^{−1}] (Eq. A1). Additional ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ sources such as the load resulting from ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ release within the stream network and point sources are not considered here (similar to Bertuzzo et al., 2017, and Wollheim et al., 2006). This assumption is supported by largescale assessments in France (Moatar et al., 2017) and Germany (Ebeling et al., 2021), where it was found that, with the exception of pure urban catchments, diffuse sources rather than point sources control the C–Q shape in typical European mixedlanduse settings. Also, we do not consider other loading processes that may create bending at the catchment outlet (e.g., shifts in transport pathways and solute sources; Marinos et al., 2020).
The modeled instream ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ uptake follows firstorder removal kinetics (Alexander et al., 2000; Boyer et al., 2006; Ensign and Doyle, 2006), such that the outgoing load from grid cell i L_{i} [mg s^{−1}] is a fraction of the incoming load L_{in,i} (Eq. 4), and the absolute removed load L_{r,i} [mg s^{−1}] can be described as (Eq. 5). Here, L_{r,i} is influenced by separate hydrological (${P}_{i}\cdot {l}_{i}/{Q}_{i}$) and biological (v_{f}) components (similar to Bertuzzo et al., 2017).
where P_{i} is the wetted perimeter of the crosssection calculated from the Manning equation (using the bed slope S_{i} and assuming a fixed roughness coefficient = 0.03 [m m^{−1}]) in open channels (Eq. A2). The uptake velocity parameter v_{f} [m d^{−1}] refers to the vertical movement of ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ molecules from the water column towards the biofilm at the pelagic–benthic interfaces and the sediments, where the instream processing chiefly occurs with v_{f}=k_{i}d_{i}, and k_{i} is the firstorder removal constant (Ensign and Doyle, 2006; Wollheim et al., 2006; Marcé et al., 2018). The parameter v_{f} accounts for the processes altering the rate and form of downstream ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ delivery (Doyle, 2005) (therefore it is not limited to denitrification only). We assume that v_{f} is independent of the instream ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ concentration C_{mean} (Pennino et al., 2014; O'Brien et al., 2007) such that the areal uptake rate $U={v}_{\mathrm{f}}\cdot {C}_{\text{mean}}$ [$\mathrm{mg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$] is tightly linked with C_{mean} in a firstorder relationship. Others (e.g., Hensley et al., 2014; Mulholland et al., 2008; O'Brien et al., 2007) contrarily found explicit scaling relationships, where v_{f} decreases nonlinearly for increasing C_{mean} (10^{−4}–10^{1} mg L^{−1}) when considering distinct catchments. However, in Germany, the ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ concentration range across a range of catchments is small (10^{−1}–10^{1} mg L^{−1} according to Ebeling et al., 2021), and rivers generally have minor longitudinal concentration variability (Hensley et al., 2014; Ensign and Doyle, 2006), which suggests independent definitions of v_{f} and C_{mean}.
The Damköhler number Da [–] is calculated as the ratio between transport (τ_{T}) and reaction (τ_{R}) (Eq. 6) timescales and is often used to characterize the relative importance of hydrological and biogeochemical processes in hydrologically connected systems (Oldham et al., 2013; Kumar et al., 2020):
where τ_{T} represents the effective travel time, TT [d], or the exposure timescale under advective conditions. We estimated the catchmentwide TT as the spatiotemporal median of the sum of all downstream T_{i} (Eq. 2d) for a grid cell in the network $\left({\sum}_{i}^{\text{Out}}{T}_{i}\right)$ (similar to Bergstrom et al., 2016), whereas τ_{R} represents the reactive timescale of biological processes. It is approximated as k^{−1} [d^{−1}], with the effective catchmentwide k estimated as the spatiotemporal median of the gridscale firstorder reaction constant ${k}_{i}={d}_{i}/{v}_{\mathrm{f}}$.
2.3 Exploring Curv_{max} with Monte Carlo simulations
Monte Carlo simulations are performed to explore how Curv_{max} evolves from a range of model input parameter combinations in a variety of catchments (Sect. 2.3.1 below). These simulations utilize the same ensemble of 11 107 unique parameter combinations in each of the individual study catchments with distinct observed discharge time series to account for the differences one parameter combination may have in each of the catchments. The unique parameter combinations are generated by Latin hypercube sampling from uniform parameter ranges that are set according to literature values (Table 1). Some physical constraints were also imposed such that the channel geometry parameters a_{w} and a_{d} must obey continuity principles (${a}_{\mathrm{w}}+{a}_{\mathrm{d}}<\mathrm{1}$ and a_{w}>a_{d}, following Leopold and Maddock, 1953). Similar to Bertuzzo et al. (2017) the chosen parameter combination is kept constant in time and uniform in space for simplicity during one model run. The simulated variables are (i) Curv_{max} [–], deduced from simulated log (C)–log (Q) relationships when a minimum of 80 % of the C data are above the “detection limit” of 0.002 mg L^{−1} ${{\mathrm{NO}}_{\mathrm{3}}}^{}$; (ii) the networkwide percentage of load removed L_{r.perc} [%], which is calculated as the median of the ratio between the daily absolute removed load and the daily absolute incoming load in the river network; (iii) the median network travel time, TT [d]; (iv) the Damköhler number Da [–]; (v) the slope of the linear regression fit of the log (C)–log (Q) relationship at the catchment outlet, b_{out} [–]; and (vi) the median concentration at the catchment outlet, C_{out} [mg L^{−1}], and the median water velocity, v [m s^{−1}]. While all outputs can be spatially and temporally explicit on a daily time step, Curv_{max} is examined at the catchment outlet, integrating both spatial and temporal aspects. The Monte Carlo results are subsequently subjected to a global sensitivity analysis with the PAWN method (Pianosi and Wagener; 2015) to elucidate influential model parameters. Furthermore a correlation analysis is conducted to explore how these influential parameters impact simulated Curv_{max}. Finally, a classification and regression tree algorithm (CART; Breiman et al., 1984) allowed us to visualize parameter interactions as detailed in Sect. 2.3.2 below.
2.3.1 Catchment selection
The river networks of 13 mesoscale catchments across Germany (areas between 21 and 1450 km^{2}; Ebeling, 2020a; Table 2) are extracted together with the corresponding daily discharge data at the outlet (uninterrupted for ∼10 years between 1995 and 2010; Musolff, 2020) as inputs for the exploratory river network model. The selected catchments include three nested subcatchments for the Selke as well as the Holtemme river system, both part of the Bode, a wellstudied river system near the Harz mountains in central Germany (Fig. 1; Ehrhardt et al., 2019; Rode et al., 2016; Winter et al., 2021; Mueller et al., 2018). All catchments have distinct geophysical settings as stream order, median discharge and catchment shape (quantified with the Horton form factor; Horton, 1945; Table 2), which is needed to obtain a realistic range of simulated Curv_{max} using the explorative network model in the Monte Carlo mode. The selected catchments were delineated in ArcMap (ESRI, 2011) from a 100 m×100 m digital elevation model (DEM) (EEA, 2013; Ebeling et al., 2021). A flow direction, flow accumulation and valley slope grid in the same resolution were established. The channel threshold drainage area for the network delineation was set to 150 grid cells (1.5 km^{2}), which agreed well with the observed river network, resulting in a treeshaped river network with N grid cells or nodes.
2.3.2 Model evaluation
To verify the network model's ability to reproduce realistic concentration time series and Curv_{max}, the observed monthly nitrate concentrations and the simulated time series were compared in 1 of the 13 selected catchments, the Selke catchment (at Meisdorf gauging station, 282 km^{2}; Table 2; Winter et al., 2021), where extensive field campaigns and modeling studies have been conducted related to instream processes (Rode et al., 2016; Dupas et al., 2017; Yang et al., 2019, 2018). Note that such a model validation was performed only in the Selke Meisdorf catchment as the model aim is to explore how certain input parameter combinations may result in log (C)–log (Q) bending at the catchment outlet rather than to reproduce catchmentspecific concentrations by performing parameter calibration. This explorative approach is comparable to Bertuzzo et al. (2017) and Helton et al. (2018), who applied a similar model in synthetic river networks but did not validate the model structure in a realistic setting. The Meisdorf catchment is a relatively homogeneous upstream part of the Selke, consisting of forest and cropland, and is characterized by constant export regimes (Winter et al., 2021). For an input parameter combination set to reasonable values for this catchment (Table C1; Rode et al., 2016), the landtostream ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ inputs averaged 1.2 $\mathrm{kg}\phantom{\rule{0.125em}{0ex}}\mathrm{N}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{km}}^{\mathrm{2}}$, which is similar to the 1.9 $\mathrm{kg}\phantom{\rule{0.125em}{0ex}}\mathrm{N}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{km}}^{\mathrm{2}}$ reported by Winter et al. (2021) for the Selke river (Meisdorf), and it is well within the general 0.001 to 100 $\mathrm{kg}\phantom{\rule{0.125em}{0ex}}\mathrm{N}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{km}}^{\mathrm{2}}$ range established by Mulholland et al. (2008). The simulated flow velocity had a spatiotemporal median value of 0.47 m s^{−1}, which is also comparable with measured flow velocities (RisseBuhl et al., 2017). Additionally, the Selke catchment was used to gain an insight into how the interplay of transport and uptake processes at every network grid cell can result in a curved C–Q pattern at the catchment outlet for one parameter combination, while in the other catchments the network model outputs were only considered at the catchment outlet for the entire range of input parameter combinations in the Monte Carlo approach. Finally, the simulated range of Curv_{max}, obtained from applying the network model in the 13 different catchments with >10 000 input parameters, was compared to a range of observed Curv_{max} computed from the C–Q relationships of 444 French catchments (Dupas et al., 2019; Sect. 2.1). Other Monte Carlo simulation outputs, such as ranges of L_{r.perc} and Da, were compared to literature values for validation purposes.
2.3.3 PAWN sensitivity analysis and correlation analysis
We performed a global sensitivity analysis (GSA) using the momentindependent PAWN method (Pianosi and Wagener, 2015). The method allowed for estimating the effect of the parameter inputs on the entire model output distribution and can be applied to rank the inputs and identify the uninfluential ones. The resulting PAWN sensitivity indices were estimated from generic input–output samples created with the numerical approximation strategy proposed by Pianosi and Wagener (2018). With this strategy, the range of variation in each input x_{i} is partitioned into a number n_{i} of equally probable “conditioning” intervals (I_{i,k}, $k=\mathrm{1},\mathrm{\dots},{n}_{i}$); i.e., each interval contains the same number of data points. Given a scalar model output y (here Curv_{max}), the PAWN method compares the output cumulative distribution function (CDF) (F_{y}(y)), computed by concurrently varying all the inputs, and the n_{i} conditional CDFs for that input (${F}_{y\mathrm{}{x}_{i}}(y\mathrm{}{x}_{i}\in {I}_{i,k})$). Each conditional CDF is obtained by varying all inputs within their entire range except for x_{i}, whose values are contained within one of the n_{i} conditioning intervals. The Kolmogorov–Smirnov statistic (KS) is then calculated as the maximum vertical distance between the conditional and unconditional CDFs, while the PAWN sensitivity index (S_{i}) for input x_{i} aggregates the results over all conditional CDFs through a summary statistic as presented in Eq. (7):
where $\text{KS}\left({I}_{i,k}\right)={max}_{y}\left{F}_{y}\right(y){F}_{y\mathrm{}{x}_{i}}(y\mathrm{}{x}_{i}\in {I}_{k,i}\left)\right$.
In this study, we applied Eq. (7) using n_{i}=10 conditioning intervals for each input parameter and used the maximum KS value, KS_{max} (ranging from 0 to 1), as a summary statistic, which is appropriate for screening noninfluential input parameters. For a given parameter, the highest value of KS_{max} of 1 would indicate a direct dependence of the model output (in this case Curv_{max}) on that parameter, while a value of KS_{max} of 0 would mean that the parameter is completely noninfluential. We estimated confidence intervals of the sensitivity indices using 15 000 bootstrap resamples and checked the robustness of the results. The PAWN analysis was carried out using the Python version of the SAFE toolbox for global sensitivity analysis (Pianosi et al., 2015).
To explore the direction of change in the C–Q bending at the catchment outlet resulting from variations in the model parameters and the catchment instream uptake, a Spearman rank correlation analysis was performed including all the simulated catchment responses and parameter combinations. These correlations were visualized in a correlation matrix using the “corrplot” package in R (Wei and Simko, 2017).
2.3.4 Identify parameter and model output interactions with classification tree
Finally, we aim to determine if, within this modeling framework, C–Q bending at the catchment outlet (specifically Curv_{max}) informs about the networkwide instream uptake. Thereto, a recursive modeling approach is proposed, using the classification and regression trees algorithm (CART; Breiman et al., 1984), which allows for the identification of nonlinear synergistic interactions among model parameters and output variables. This nonparametric method segregates classes for a response variable by progressively splitting selected predictor variables in a binary way. The resulting decision tree is intuitive to interpret and can facilitate the fast characterization of river networks. The response variables include the effective catchmentwide removal efficiency L_{r}, the Damköhler number Da and the uptake velocity v_{f}, while the predictors are Curv_{max}, the median network velocity v and all of the model input parameters except for v_{f} (Table 1). For each response variable, three classes are defined representing low, intermediate and high ranges found in the literature (Table 3) that each contain 5 % of the simulation outputs (obtained by distributing the nonmissing model simulations over 20 percentiles). The overall CART accuracy for each response variable is assessed by attributing 80 % of the simulation outputs in the low, intermediate and high classes to a training sample and assigning the remaining 20 % to a test sample. The training sample is then used to construct the classification tree, while the test sample is needed to assess the prediction accuracy and calculate the performance statistics for each class. The CART analysis was performed using the “caret” package in R with the Gini impurity measure as a splitting criterion (Kuhn, 2020).
3.1 Empirical Curv_{max}
The estimated Curv_{max} for the French observed ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ log (C)–log (Q) data (Dupas et al., 2019) ranges between −5.25 and 3.88 (median is −0.23; Fig. B4), and 77 % of the stations are characterized by Curv_{max}≤0 or a linear or concave shape (similar to Moatar et al., 2017). The time series subsamples for each station generally had a small Curv_{max} variability (interquartile range, IQR, for a given station below 1) for 93 % of the stations, with some exceptions demonstrating a larger IQR up to 8. This indicates that Curv_{max} quantification for most lowfrequency C–Q time series is robust. The Spearman rank correlation (ρ=0.53, $p\text{value}\mathrm{2.2}\times {\mathrm{10}}^{\mathrm{16}}$) between the absolute observed Curv_{max} and IQR for each station is significant and positive, implying that C–Q relationships with a higher absolute Curv_{max} have a higher uncertainty when quantifying the C–Q bending. However, Curv_{max} variability (IQR) in the subsamples for each station has no significant correlation with the number of data points available for one station. This implies that Curv_{max} tends to be temporally robust when the C–Q data obey the four criteria in Sect. 2.1 so that the length of the lowfrequency time series does not impact the estimated Curv_{max}. Overall, the proposed Curv_{max} metric is suitable to quantify bending in multiannual, temporally stable log (C)–log (Q) relationships.
3.2 Model evaluation in the Selke river (Meisdorf)
To evaluate the network model performance in a realistic setting, we implemented the model with a fixed parameter combination (Table C1) in the Selke catchment and aimed to capture C–Q dynamics at the basin outlet. The simulated ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ concentration time series (simulated C with uptake) for the Meisdorf station in Fig. 2a shows a seasonal pattern that follows the observation concentration data reasonably well (Nash–Sutcliffe efficiency, NSE=0.50; percent bias, $\text{pbias}=\mathrm{0.4}\phantom{\rule{0.125em}{0ex}}\mathit{\%}).$ This seasonality is also reflected in simulated daily percentage of load removed (the ratio between the daily total removed load and the daily total incoming load in the river network) and ranges from almost 0 % to 3.4 % in this case, with the median L_{r.perc} value equal to 0.41 %. The highest removal efficiencies are simulated in fall and summer and coincide with low simulated ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ concentrations at the catchment outlet. The conservative ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ concentration (simulated C no uptake) is stable around 3 mg L^{−1}. The observed nitrate concentrations generally show an enrichment export pattern in the log (C)–log (Q) space (b_{out}=0.40, R^{2}=0.56) and a Curv_{max} of −0.35, which agrees well with the simulated Curv_{max} of −0.28 (Fig. 2b) and deviates significantly from the conservative scenario of simulated C without uptake (b=0.014; Table C1). The observed low nitrate concentrations coincide with low discharges in fall and summer, while high concentrations occur mainly in winter, when discharges are higher.
Within the Selke Meisdorf river network the simulated Curv_{max} is largely contained within −1.12 to −0.29 (10th and 90th quantiles, respectively) for the given parameter combination (Table C1, Fig. 3). Low Curv_{max} ($<\mathrm{1.12}$) is found exclusively at grid cells with a low total drainage area (A_{i}<9 km^{2}), and Curv_{max} stabilizes at higher values with increasing drainage areas (inset Fig. 3, Fig. S3 in the Supplement). The incoming (L_{in.ls} and L_{in.up}; Eq. 3), removed (L_{r}; Eq. 5) and outgoing absolute load (L_{i}; Eq. 4 with L=CQ) as a function of Q in the log–log space are shown in Fig. 3 for three selected grid cells on the main river stem with low (C), intermediate (B) and high (A) drainage areas. The corresponding log (C)–log (Q) relationship for the outgoing load (L_{i}) at the outlet (A) is presented in Fig. 2b for simulated and observed values. Note that Curv_{max} is calculated from log (C)–log (Q) relationships rather than log (L)–log (Q). The loads in grid cell A, B and C generally increase with discharge, while the load removal efficiency decreases with discharge. The highest removal efficiencies are found in the headwater grid cell C (39 % for low discharge), followed by midstream grid cell B (3 % for low discharge) and the outlet A (0.5 % for low discharge). The total absolute load removed (L_{r}, sum per year per grid cell) is largest for firstorder grid cells (average 24.1 kg N yr^{−1}) that represent 55 % of the river network, followed by second and thirdorder grid cells (averaging both around 20 kg N yr^{−1}) that represent 20 % and 25 % of the network (inset Fig. 3). Finally, the total yearly incoming load (L_{in.ls}+L_{in.up}, sum per year per grid cell) increases with stream order from 1329 kg N yr^{−1} on average in a firstorder grid cell to 5128 and 42 124 kg N yr^{−1} in secondorder and thirdorder river cells.
With uniform, constant parameters the network model does not account for a spatiotemporal parameter variability. Nevertheless, it successfully (see NSE and pbias) reproduces the seasonality of the observed concentrations over the 2000–2010 period for the Selke Meisdorf catchment (Fig. 2a). For comparison, Yang et al. (2018) found a similar performance (NSE=0.47, $\text{pbias}=\mathrm{3.35}$ %) when applying a fully distributed model with 16 calibrated parameters in this catchment between 1997 and 2009. The uptake velocity v_{f} for our simulation was set to 0.098 m d^{−1} to closely match the observed (assimilatory) uptake range of 0.009 to 0.103 m d^{−1} for the Selke Meisdorf river network (Rode et al., 2016); the median annual percentage of load removed equals 4.7 %, which is within a comparable range reported in prior studies (Rode et al., 2016, and Yang et al., 2018, found annual means of 4.8 % and 7.6 %, respectively). Note that the annual percentage of load removed accounts for load taken up throughout the entire river network, which may be higher in the headwaters (15 t) than in downstream locations (7 and 5 t for second and thirdorder stream sections; inset Fig. 3). Yang et al. (2018) reported very high uptake efficiencies (up to 75 %) for summer seasons that were caused by low ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ concentrations (0.21 mg N L^{−1}) and load (L=CQ), which are not represented in our model simulation (the lowest simulated ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ concentration equaled 0.4 mg N L^{−1}). Additionally, due to the parsimonious structure of the proposed model, we did not account for the temporally changing effects of environmental factors like temperature and light availability that might (seasonally) influence uptake efficiencies in the river network (Kadlec and Reddy, 2001). Nevertheless, these high uptake efficiencies under low flows in summer are not a main contributor to the annual percentage of load removed that is dominated by high flows, generally recorded during winter. Thus for the Monte Carlo simulations we calculated L_{r.perc} as the median of the daily percentage of load removed rather than the total removal efficiency for the entire simulated time period to better represent an effective longterm networkwide removal capacity.
The interplay of incoming, removed and outgoing load at each network grid cell shapes the log (L) and log (C)–log (Q) relationships and consequently Curv_{max} at the catchment outlet (Fig. 3). Landtostream loading (L_{in.ls}) that varies linearly with direct incoming discharge at a given grid cell in the log space (Eq. 3 with L=CQ; Curv_{max}=0) can lead to a bent outgoing log (C)–log (Q) relationship where concentration or load (L_{i}) varies nonlinearly with discharge (Curv_{max}≠0). The onset of a bent log (C)–log (Q) pattern (${\text{Curv}}_{\text{max}}=\mathrm{0.37}$) is illustrated in the headwater grid cell C in Fig. 3, where L_{in.ls} is the only incoming load (upstream incoming load, L_{in.up}, equals 0 in this case). The absolute removed load in a grid cell is higher under increasing Q, while the percentage of load removed is lower, which explains observed C–Q patterns with higher log (C)–log (Q) slopes for low flows than for high flows (Moatar et al., 2017; Wollheim et al., 2008, 2017; Doyle, 2005; Basu et al., 2011). This decreased ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ load removal efficiency in the downstream direction (spatial scale) or during events (temporal scale) can arise because stream morphology characteristics such as depth and water velocity that correlate with varying discharge constitute higher surfacetovolume ratios at low flows (generally in the headwaters) than at higher flows (at the outlet) (Peterson et al., 2001; Hensley et al., 2014). The uptake and landtostream loading at the downstream grid cells (B and A in Fig. 3) have a decreasing local impact on the outgoing load due to the large upstream load contributions that increase in the downstream direction (see explicit scaling relationship for input flux in Bertuzzo et al., 2017). This is also explained by Wollheim et al. (2018), who suggest that the river network saturates as supply exceeds biological “demand”, causing the log (C)–log (Q) relationship to approach the slope of the loading function (presented as the simulated C without uptake in Fig. 2b). Dupas et al. (2017) on the other hand show how ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ uptake effects are decreasingly visible in C–Q observations downstream, and concentrations largely matched those estimated by a conservative mixing model. The saturation effect with the accumulation of large load is reflected in the Curv_{max} converging to a constant value when moving from upstream to downstream or from a lowerorder to a higherorder river reach (Figs. 3 and S3). This also corroborates the recent findings of Abbott et al. (2018), who found that the temporal variability (here reflected in the C–Q relationship) in nutrients is preserved moving downstream in a river network. Overall the Selke example shows that the network model can realistically reproduce the bending of observed ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ C–Q relationships that evolve from the decreasing removal efficiency at higher discharges.
3.3 Monte Carlo simulation results
The overview of the model outputs for each of the 13 study catchments in Table 4 shows that catchments 1, 5 and 11 display the lowest 10th quantile Curv_{max} values of −1.61, −1.40 and −1.24 (more bending), while the catchments 4 and 6 registered higher (less bending) and less variable Curv_{max} (10th quantiles at −0.31 and −0.35) (Fig. B5). Catchments 3, 4 and 8 are characterized by high runoff and Q (Table 2) at the catchment outlet and demonstrate low percentages of load removed, L_{r.perc} (90th quantile at 29.8 %, 32.1 % and 19.3 %, respectively). The highest L_{r.perc} values are found in catchments 1 and 10 (98.4 % and 95.1 % for the respective 90th quantiles). The regression slope of the log (C)–log (Q) relationship at the basin outlet, b_{out}, is positively skewed for all the catchments (most positive slopes found in catchment 5), while the slope b of the landto stream loading function had no positive or negative preference (Table 1, Eq. 3). The distribution of the concentrations at the catchment outlet, C_{out}, are generally similar across all catchments (10th and 90th percentiles within 0 to 6.2 mg L^{−1}) and are significantly less variable than the landtostream incoming concentration (parameter C_{mean}) that varied from 10^{−4} to 20 mg L^{−1} across all the simulations (Table 1). The highest C_{out} values are found in catchment 8, the largest catchment. The median water velocity v (Eq. 2c) is between 0.01 and 0.5 m s^{−1} for the 10th and 90th quantiles of all the study catchments, and the largest v is simulated for catchments 3 and 4, which also have the highest discharge. The median river network travel time, TT, for all simulations and catchments ranges from 0.1 to 4 d between their respective 10th and 90th quantiles and remarkably has no clear relationship with catchment properties such as the total river network length (Table 2). Finally, the Damköhler number, Da (Eq. 6), is variable around 1 with the highest values, indicating reactiondriven conditions, found for catchments 2 and 12 (respective ranges from 0.6 to 10.3 and 0.7 to 10.8 for the 10th and 90th quantiles). The lowest Da values are found for catchments 4 and 10 (90th quantile <2), implying more transportdriven conditions.
The Monte Carlo output in Table 4 shows reasonable values for the different variables, taking into account that the goal of this modeling exercise was not to reproduce catchmentspecific conditions but rather to explore how ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ uptake influences C–Q bending for a range of parameter combinations that represent a spectrum of possible catchment conditions. The simulated Curv_{max} values for all 13 German study catchments and parameter combinations (80 % of the values between −0.70 and −0.012; Table 4 and Fig. B5) are comparable with the range of Curv_{max} from ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ log (C)–log (Q) relationships in the French catchments (80 % of the values between −0.41 and −0.067; Fig. B4) (Dupas et al., 2019). Simulated Curv_{max} is always smaller than or equal to zero as explained in Sect. 2.2.2. For the model output L_{r.perc}, a wide range of uptake efficiencies were captured from almost 0 % to near to 100 % for some simulations and a median value of 14.4 % across simulations. This simulated range exceeds the proposed range by Birgand et al. (2007) of 10 % to 70 % of N removal for agricultural drainage networks at annual timescales. High removal percentages (median over the simulated time period of daily percentage of load removed in the network exceeding 95 %) are registered for 3.4 % of all simulations, while very limited load removal (L_{r.perc}<5 %) occurred for 32.1 % of all the simulations. Other simulation outputs such as the effective velocity v surprisingly rendered similar distributions across the catchments (Table 4) given that the median Q varied for almost 3 orders of magnitude at the basin outlet (Table 2). Their specific discharges (Sect. 2.2.1) however are similar, and by taking the spatiotemporal median v as an effective catchment value for each simulation the (more numerous) headwater grid cells were better represented than the grid cells close to the basin outlet. A similar effect is found for the range of the effective travel time TT. Generally these similar v and TT distributions from model simulations between catchments align with the notion of Langbein and Leopold (1964) that drainage networks evolve naturally to transport water (and sediment) most efficiently such that an equilibrium between channel form and water and sediment load is imposed (Leopold and Maddock, 1953). Also Damköhler numbers Da exhibited realistic ranges, mostly distributed around 1 (Fig. S3; Oldham et al., 2013; Ocampo et al., 2006), with 36.5 % of the simulations <0.8 % and 50.8 % > 1.2, indicating that more simulations are reactiondriven than transport driven. Finally, other variables such as the range of the modeled river widths are found to be reasonable to a large degree (Fig. S4 in the Supplement).
As for the simulations, the same >10 000 parameter input sets are applied in each catchment, differences between the catchments result from the different network structures and hydrological regimes that control transport and uptake processes for each parameter input set in each catchment. From Tables 4 and 2 it is clear that differences in model outputs (i.e., Curv_{max}, L_{r.perc} and Da) between the catchments cannot be attributed to a single catchment property such as total network length or basin area. For example Curv_{max} has the highest variability between simulations in the catchment 1, the smallest catchment, compared to the other catchments, which could be attributed to the variability in local loading and uptake patterns in the network (driven by Q) that are still visible at the catchment outlet. Following the simulated Selke Meisdorf example in Sect. 3.1 (Figs. 3 and S3), we show that Curv_{max} tends to converge to a constant value with increasing drainage areas (similar to Abbott et al., 2018, for nutrient concentrations; Dupas et al., 2017, for nutrient uptake; and Bertuzzo et al., 2017, for DOC removal). Drainage area is however not the only catchment property influencing Curv_{max} at the outlet. For example, catchment 6 is the secondlargest catchment (Table 2) and has the least bent (and least variable) log (C)–log (Q) relationships. The network structure could possibly play a role here as the largest catchment, catchment 8, has some large tributaries near the basin outlet (Fig. 1), which could bypass removal and transport high load during events, introducing a more variable Curv_{max} (Mineau et al., 2015; Helton et al., 2018). The percentage of load removed, L_{r.perc}, is notably lower for catchments with high runoff or Q – like 3, 4 and 8 (Table 4) – which corresponds with findings in Sect. 3.1 that uptake efficiency decreases with increasing Q because of increasing loads to the system (Wollheim et al., 2018; Mulholland et al., 2008) and that increasing Q results in less efficient uptake within the reactive surface area (Peterson et al., 2001; Hensley et al., 2014). The high L_{r.perc} in small catchments 1 and 10 could then be attributed to their low Q; however why the small catchment 5 does not have similar uptake performance is not evident. Generally the model output variability between the catchments (as a result of different network and discharge properties) is minor compared to the output variability within the catchments (due to the effect of the chosen input parameter set).
3.4 Curv_{max} sensitivity analysis and model parameter correlation
The PAWN sensitivity index KS_{max} in Fig. 4 and Table C2 shows that across all catchments Curv_{max} is most sensitive to the exponents in the width–Q relation a_{w} (KS_{max}=0.62) and depth–Q relation a_{d} (KS_{max}=0.51). Here, there is little variability between the catchments (KS_{max} has a low coefficient of variation, CV, of 0.06 and 0.22, respectively). Overall, the slope of the linear loading function, b, is least important in shaping Curv_{max} (KS_{max}=0.14). Nevertheless, a high variability in KS_{max} is observed (CV=0.76) that is caused by larger sensitivities for catchments 1 and 12 (KS_{max} near 0.45). Curv_{max} is equally sensitive to v_{f} and C_{mean} (KS_{max} 0.18 and 0.19), but v_{f} exhibits higher variability in KS_{max} than C_{mean} (CV of 0.59 and 0.47). Furthermore, over all the catchments Curv_{max} is sensitive to the median velocity v and the Damköhler number Da (KS_{max} equals 0.64 and 0.31, respectively; CV of 0.26 and 0.38). When considering the catchments individually, basin 1, with the smallest discharge, has the highest median KS_{max} (0.59) across all input parameters, while catchment 4, which has the highest discharge, exhibits the lowest median KS_{max} (0.13). Additionally, Curv_{max} is very sensitive to the velocity v in catchment 1 (KS_{max}=0.95), while it is least sensitive to v in catchment 4. Nevertheless, the other catchments show no clear order in Curv_{max} sensitivity according to catchment properties such as Q. For example in nested catchments 10, 11 and 12 (Fig. 1), the largest catchment, catchment 12, has the highest KS_{max} (0.50) and lowest CV (0.26) over all the input parameters, indicating that here Curv_{max} is more sensitive to the input parameters than in the smaller subcatchments 10 and 11.
In a next step the estimated Curv_{max} across all simulations is correlated to the model input parameters as well as to output variables like the percentage of load removed (L_{r.perc}) the log (C)–log (Q) slope at the catchment outlet (b_{out}), the median concentration at the basin outlet (C_{out}) and the uptake constant (k) to identify the strength and direction of their relationship. The resulting Spearman correlation matrix (Fig. 5) reflects the PAWN sensitivity findings, with the highest Curv_{max} correlation found with parameters a_{w} (ρ=0.68) and a_{d} (ρ=0.56) and input variable v (ρ=0.57). Curv_{max} is independent of v_{f} ($\mathit{\rho}=\mathrm{0.04}$) but shows a negative correlation with L_{r.perc} ($\mathit{\rho}=\mathrm{0.36}$), suggesting that lower Curv_{max} (more bending) can be related to a higher L_{r.perc}. Furthermore, Curv_{max} is negatively correlated to the log (C)–log (Q) regression slope at the catchment outlet b_{out} ($\mathit{\rho}=\mathrm{0.28}$) such that higher bending coincides with more positive b_{out}. The variable v is additionally strongly negatively correlated with L_{r.perc} ($\mathit{\rho}=\mathrm{0.87}$), so a high percentage of load removed occurs at low velocities. Da on the other hand is positively correlated to L_{r.perc} (ρ=0.58), which indicates that higher Da values are occurring together with higher load removed. Da thereby seems to be controlled more tightly by variation in k^{−1} ($\mathit{\rho}=\mathrm{0.71}$) than by TT (ρ=0.48). Finally C_{out} is negatively correlated with L_{r.perc} ($\mathit{\rho}=\mathrm{0.82}$) and Da ($\mathit{\rho}=\mathrm{0.61}$).
The PAWN and correlation analysis results suggest that the input parameters dictating the channel morphology, a_{w} and a_{d} (Sect. 2.3), are controlling factors for the magnitude of the bending in log (C)–log (Q) relationships at the catchment outlet. More specifically parameters a_{w} and a_{d} influence the response of the wetted perimeter (P_{i}; Eq. A2) in a given reach in the network and drive the reactive surface area (P_{i}⋅l_{i}) with changes in discharge (Eq. 2a and b, Figs. 5 and B3). The absolute load removed L_{r,i} (Eq. 5) can be written with the width and depth exponents a_{w} and a_{d} explicitly (Eq. A3) so that ${L}_{\mathrm{r},i}\sim \mathrm{1}/\left({Q}^{\mathrm{1}{a}_{\mathrm{w}}{a}_{\mathrm{d}}}\right)$. When the denominator is large (small a_{w} and a_{d}) the effect of low and high Q's on the local absolute removed load increases and can lead to a lower Curv_{max} (more bending; Sect. 3.1, Fig. B6). Networkbased modeling studies often set the width exponent a_{w} to a value of 0.5 that was found to be representative of rivers globally (Bertuzzo et al., 2017; Rode et al., 2016; Wollheim et al., 2018). This a priori fixed a_{w} may, however, strongly affect the simulated C–Q dynamics at the basin outlet as is demonstrated here. Curv_{max} finally shows the lowest sensitivity to the loading parameters b and C_{mean} that influence the incoming load to a grid cell (Eq. 3) and thus impact the local absolute load removed L_{r,i} (Eq. 5) rather than the percentage of removed load L_{r.perc}. This indicates that the contribution of local incoming load in the downstream direction has a limited impact on the log (C)–log (Q) bending at the catchment outlet. For example in the Selke Meisdorf catchment, the locally contributing Q's are generally smaller (or equal for the headwaters) than the total Q in a given reach so that the influence of the loading parameters b and C_{mean} on the total load decreases in downstream reaches (Sect. 3.1, Fig. 3).
Although Curv_{max} only has an intermediate sensitivity to the uptake velocity v_{f}, and they do not correlate well, v_{f} is an important “boundary condition” for log (C)–log (Q) bending at the catchment outlet. No biological demand (low v_{f}) would mean that none of the incoming load would be removed from the river network. The outlet signal would in this case be solely driven by the dischargecontrolled transport processes, and no bending would be observed (Curv_{max}=0). Although increasing v_{f} does correlate with decreasing concentrations ($\mathit{\rho}=\mathrm{0.34}$) and increasing load removed (ρ=0.34), it does not always lead to more bending, as illustrated in Fig. B6 for the Selke example. Because v_{f} is defined as a constant within one simulation that is independent of the local nutrient concentration (Sect. 2.2.2), the percentage of load removed in the network is mainly controlled by the varying hydrological conditions here, represented by the effective networkwide velocity v (L_{r.perc} and v, $\mathit{\rho}=\mathrm{0.82}$). This confirms that discharge and channel morphology are among the most important predictors of removal (Alexander et al., 2000; Seitzinger et al., 2002; Wollheim et al., 2006). The role of v was further examined in the context of restored and channelized streams (Kunz et al., 2017) and agrees with our findings that decreased v influences N cycling (Peterson et al., 2001).
The PAWN and correlation analysis results show that Curv_{max} is sensitive to the Damköhler number Da (KS_{max}=0.31; Fig. 4, Table C2), which has a high positive correlation with the percentage of load removed L_{r.perc} (ρ=0.58; Fig. 5). This indicates that high Da occurs concurrently with more efficient removal and is in line with others (Ocampo et al., 2006) who found sometimes almost 100 % ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ removal in the riparian zones of an agricultural catchment with Da exceeding 2. The transport timescale TT that makes up Da (ρ=0.48; Fig. 5) together with the inverse of the firstorder uptake constant k^{−1} ($\mathit{\rho}=\mathrm{0.71}$; Eq. 6) is examined for classes of low, median and high Da (defined in Table 3) in Fig. 6a to disentangle which values of k^{−1} and TT occur together and can constitute a certain Da range (each class contains 5 % of all simulations). It is shown here that low Da values are driven by both low TT and high, variable k^{−1}, implying a transportdriven system with limited ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ removal (median L_{r.perc} equals 2.4 % in Fig. 6a for low Da). High Da values, contrarily, have high TT and low k^{−1}, fostering intermediate uptake percentages (median L_{r.perc}=27.1 %). Although also v_{f} clearly differentiates for classes of low, medium and high Da in Fig. 6a, the corresponding Curv_{max} values are similar in their range and mean. Nevertheless, this does not mean that Da is not influencing Curv_{max} at the basin outlet as there could be interactions with other inputs that are not captured here (which is supported by the PAWN findings, where Da appears to be influential).
From the Curv_{max} perspective (Fig. 6b) we identify model output ranges of L_{r.perc}, Da and input variable v_{f} that constitute low, median and high Curv_{max} classes (Table 3). High Curv_{max} (less bending, $\sim \mathrm{0.02}$) is thereby linked to low L_{r.perc} (median 4.8 %), while low Curv_{max} (more bending, $\sim \mathrm{0.60}$) is connected to higher and more variable L_{r.perc} (median 33.6 %), generally indicating that more bent systems are more efficient in terms of removal and vice versa. To explore some cases when this latter statement might not be true, we examine the input parameter ranges where more bent simulations (${\text{Curv}}_{\text{max}}<\mathrm{0.51}$, 13.1 % of all simulations) occur concurrently with low percentage of removal (L_{r.perc}<5.2 %, 0.9 % of all simulations) on the one hand and high percentage of removal (L_{r.perc}>63.0 %, 4.9 % of all simulations) on the other hand in Fig. B7a. Here, it is seen that highbending, lowuptake cases mainly occur for simulations with a high effective velocity v (driven by lower values for the channel shape parameters K_{w}, K_{d}, a_{w} and a_{d}). Low a_{w} and a_{d} are correlated with more bending (low Curv_{max}), and Curv_{max} is most sensitive to these parameters. However, low a_{w} and a_{d} do not lead to a more efficient ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ uptake if the other channel shape parameters K_{w} and K_{d} cause relatively high velocities (median v>0.1 m s^{−1}) throughout the network. The latter case is shown to be true for a minor percentage of all simulations (0.9 %) and explains why low Curv_{max} (more bending) can be connected to a wider range of L_{r.perc}. Figure B7b shows that concurrent simulations of less bending (${\text{Curv}}_{\text{max}}>\mathrm{0.03}$) and high removal (L_{r.perc}>63.0 %) are even rarer (0.1 % of all simulations) compared to concurrent less bending (${\text{Curv}}_{\text{max}}>\mathrm{0.03}$) and low removal (L_{r.perc}>5.2 %; 7.4 % of all simulations). Deviations from the expected high Curv_{max}–low L_{r.perc} pattern are also here driven by (very low) v. In the latter case however, a_{w} and a_{d} are generally high (leading to high Curv_{max}), and the different v values stem from coefficients K_{w} and K_{d}, which are higher in highremoval simulations. Finally, Fig. 6d illustrates that low, medium and high uptake velocities v_{f} lead to distinct Da and L_{r.perc} but do not show up in the bent signal at the catchment outlet.
3.5 Predicting instream processing with Curv_{max}
To determine if observed C–Q bending at the catchment outlet (here Curv_{max}) can be utilized to quantify instream uptake in the upstream river network and to visualize model parameter interactions, a classification tree was established for low, medium and high values (Table 3, Fig. 6) of the response variables L_{r.perc}, Da and v_{f} (Fig. 7). The prediction accuracy metrics in Table C3 and the probability histograms in Fig. 7 show that L_{r.perc} can be predicted relatively well (overall accuracy of 0.66) compared to the other response variables Da (accuracy 0.51) and v_{f} (accuracy 0.40). The fitted CART models all perform significantly better than a random allocation of simulation results to each class for each response variable (accuracy > no information rate, $p\text{value}\mathrm{2.2}\times {\mathrm{10}}^{\mathrm{16}}$). While the classes for L_{r.perc} and v_{f} are partitioned using only the network effective velocity v and Curv_{max}, predicting Da in our case requires information on the channel geomorphology parameters width coefficient K_{w} and depth exponent a_{d}. The histograms for each of the response variables in Fig. 7 indicate the probability of a test sample being of a certain class when following the partition rules in the respective decision tree. For example, for L_{r.perc} the probability that the daily percentage of load removed is small (around 8 %) exceeds 0.95 when the effective velocity v in the catchment is larger than 0.22 m s^{−1}, while the probability that L_{r.perc} is high (around 70 %) in this case is close to 0 (Fig. 7a). For v_{f} the lowermost (1) and highest (20) classes are predicted most accurately (0.58 and 0.56, respectively; Table C3) and indicate that when the velocity is not very small, and Curv_{max} is smaller than −0.51 (more bent), v_{f} is most likely high (probability 0.59). For Da, the lower and higher classes can be predicted most accurately (0.69 and 0.68, respectively); for example, Da is small with a probability of 0.58 when K_{w} is relatively low (<6.8). When on the other hand K_{w} exceeds 6.8 and a_{d} is larger than 0.4 or when a_{d} is smaller than 0.4 but Curv_{max} is smaller than −0.45 and v is very small (<0.04 m s^{−1}), it is most likely that Da is large.
These findings demonstrate that log (C)–log (Q) bending at the catchment outlet, together with the median velocity and the response of the width and the depth of a channel to discharge (parameters K_{w}, K_{d}, a_{w} and a_{d}) can help to classify the instream daily percentage of load removed L_{r.perc}, the Damköhler number Da and to a certain extent the uptake velocity v_{f}. This conclusion depends of course on the initial assumptions of our model setup (e.g., linear landtostream loading vs. Q with a constant slope b and no dominant influence of waste water sources as well as no seasonality in uptake velocity v_{f}). The velocity may be computed from the channel shape, discharge (Eq. 2c) and the topography with the channel shape parameters that are sometimes available from rating curve information or detectable from highresolution satellite pictures. The CART models could help obtain an initial probability of ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ removal efficiency in a river network, especially in a context where networkwide uptake measurements are scarce (Wollheim et al., 2017; Hensley et al., 2014), and physical, fully distributed models are not always feasible to apply (Boyer et al., 2006; Klemes, 1986). Although the CART models are developed using “only” the 13 German catchments included in the Monte Carlo analysis, in Sect. 3.3 and Table 4 we have shown that the output variability between the catchments (as a result of different catchment properties) is minor compared to the output variability within the catchments (due to the effect of the input parameter set), thus justifying the CART model use. Nevertheless, the prediction performance of these CART models might be influenced in unknown ways when applied to catchments with dissimilar catchment sizes, network structures or hydrological regimes.
In this study, we explore how lowfrequency ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ log (C)–log (Q) relationships, observed at a basin outlet, can display bending as a result of networkscale instream uptake processes. We established a parsimonious gridbased river network model for 13 distinct German catchments and investigated the influence of instream loading, transport and uptake parameters on the bending of log (C)–log (Q) relationships. Based on our exploratory analysis we conclude the following.

Noisy, multiannual and lowfrequency ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ log (C)–log (Q) relationships at a basin outlet can be described as bent, and the amount of bending can be robustly quantified with the new Curv_{max} metric. Curv_{max} is temporally stable on multiannual timescales and can be computed alongside existing log (C)–log (Q) descriptors, such as the slope of the linear regression model.

A bent log (C)–log (Q) relationship (Curv_{max}<0) at the basin outlet can arise from log–log linear landtostream C–Q relationships if uptake is present within the river network (v_{f}≠0). This supports the hypothesis that more positive slopes under low flow (bent log (C)–log (Q) curves) are linked to biological ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ concentration mediation in the stream (Moatar et al., 2017) and connects Curv_{max} (as a quantitative measure) to observations of increased removal efficiency under low flows (Wollheim et al., 2017). Our findings also stress the need to monitor the entire discharge range and capture low flows as well as high flows in a catchment.

The bending at the catchment outlet is primarily shaped by the channel geomorphological parameters a_{w} and a_{d} (exponents in the respective stream width and depth to discharge relationships, with Curv_{max} sensitivity indices KS_{max} equal to 0.62 and 0.51 and Spearman correlation coefficient ρ equaling 0.68 and 0.56, respectively) and less by the uptake velocity v_{f} (KS_{max}=0.18, $\mathit{\rho}=\mathrm{0.04}$), given that v_{f} differs from zero. In the latter case Curv_{max} would equal zero, and the log (C)–log (Q) relationship would be solely shaped by the accumulation of upstream load. Thus, the change in reactive channel bed area with discharge (mediated by a_{w} and a_{d}) has a greater influence on the bending at the outlet than the biological removal capacity (here v_{f}). Additionally we demonstrate that an a priori fixed a_{w} might strongly affect the simulated C–Q dynamics at the basin outlet. This calls for a better representation of channel geomorphological parameters at stream networks to allow for better modeling assessments.

Curv_{max} at the basin outlet can be linked to the networkwide removal efficiency L_{r.perc} ($\mathit{\rho}=\mathrm{0.36})$ under certain conditions, generally showing that systems with more bending in their log (C)–log (Q) relationship are more efficient in terms of removal and vice versa. It is, however, clear that also cases with high bending (${\text{Curv}}_{\text{max}}<\mathrm{0.51}$) and low removal (L_{r.perc}<5.2 %, 0.9 % of all simulations) or low bending (${\text{Curv}}_{\text{max}}>\mathrm{0.03}$) with high removal (L_{r.perc}>63.0 %, 0.1 % of all simulations) exist that are imposed by respective higher and lower networkwide median velocities. This shows how the velocity, v, (calculated from the channel shape parameters a_{w}, a_{d}, K_{w} and K_{d}) may mediate the connection between L_{r.perc} and Curv_{max} and indicates that v should be considered when interpreting log (C)–log (Q) bending. Consequently, anthropogenic impacts in terms of channelization of river networks might lead to lower removal efficiencies.

Classification trees – like CART – can be useful for predicting low, median and high classes of response variables L_{r.perc}, the Damköhler number Da and v_{f}. They provide useful insights on how catchments with lowfrequency concentration and discharge time series (that are generally available) can reveal information on the upstream river network uptake performance.
To evaluate the generality of the results presented here, Curv_{max} should be calculated for ${{\mathrm{NO}}_{\mathrm{3}}}^{}$ concentration observations of a larger range of catchments and linked to the respective catchment properties. Properties such as light and stream ecological state can serve as proxies for uptake performance, and for example topographic gradient can be a proxy for network transport velocity. Finally, including conservative tracers in the analysis can be used to estimate loading scenarios. Such a datadriven exploration would further elucidate the linkages between nutrient uptake efficiency and lowfrequency C and Q observations.
Calculation of c is performed as follows (Jawitz and Mitchell, 2011):
with
Stream channel wetted perimeter P_{i} [L] (where A is the crosssectional area [L^{2}]; R_{H} [L] is the hydraulic radius; and w_{i} [L], d_{i} [L] and v_{i} [L T^{−1}] are the local stream width, average depth and velocity, respectively) is calculated in Eq. (A2). S_{i} [L L^{−1}] is the stream bed slope, and n [–] is the Manning roughness coefficient that is equal to 0.03 for all simulations.
The load removed in a grid cell (Eq. 5) with the width and depth exponents a_{w} and a_{d} stated explicitly is calculated as follows:
The velocity in a grid cell (Eq. 2c) with the width and depth exponents a_{w} and a_{d} stated explicitly is calculated as follows:
The network model code and the CurvMax function can be accessed at https://doi.org/10.4211/hs.da70a09dc6074242ada756c29d12dcb3 (Dehaspe, 2021).
The data we used were published by third parties and referenced in the text where appropriate, with complete access data in the reference list.
The supplement related to this article is available online at: https://doi.org/10.5194/hess2564372021supplement.
JD and AM designed the study together with RK. JD developed the model code, carried out the simulations and interpreted them. FS provided help with the PAWN method. JD and AM prepared the manuscript draft, and all coauthors contributed to reviewing and editing the manuscript.
The authors declare that they have no conflict of interest.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This research was supported by funding from the InStream cohort and discussions with the participating scientists. Support to Fanny Sarrazin was provided by the Reduced Complexity Models project cofunded by the Helmholtz Association, and Fanny Sarrazin and Rohini Kumar acknowledge the Advanced Earth Modeling Capacity (ESM) project funded by the Helmholtz Association. The authors further thank Remi Dupas for providing the French nitrate CQ data, Pia Ebeling for the German catchment characteristics data and Nima Postanchi for his endless support.
The article processing charges for this openaccess publication were covered by the Helmholtz Centre for Environmental Research – UFZ.
This paper was edited by Christian Stamm and reviewed by Wilfred Wollheim and one anonymous referee.
Abbott, B. W., Gruau, G., Zarnetske, J. P., Moatar, F., Barbe, L., Thomas, Z., Fovet, O., Kolbe, T., Gu, S., PiersonWickmann, A. C., Davy, P., and Pinay, G.: Unexpected spatial stability of water chemistry in headwater stream networks, Ecol. Lett., 21, 296–308, https://doi.org/10.1111/ele.12897, 2018.
Aguilera, R., Marcé, R., and Sabater, S.: Modeling nutrient retention at the watershed scale: Does small stream research apply to the whole river network?, J. Geophys. Res.Biogeo., 118, 728–740, https://doi.org/10.1002/jgrg.20062, 2013.
Alexander, R. B., Smith, R. A., and Schwarz, G. E.: Effect of stream channel size on the delivery of nitrogen to the Gulf of Mexico, Nature, 403, 758–761, https://doi.org/10.1038/35001562, 2000.
Alexander, R. B., Böhlke, J. K., Boyer, E. W., David, M. B., Harvey, J. W., Mulholland, P. J., Seitzinger, S. P., Tobias, C. R., Tonitto, C., and Wollheim, W. M.: Dynamic modeling of nitrogen losses in river networks unravels the coupled effects of hydrological and biogeochemical processes, Biogeochemistry, 93, 91–116, https://doi.org/10.1007/s1053300892748, 2009.
Ameli, A. A., Beven, K., Erlandsson, M., Creed, I. F., McDonnell, J. J., and Bishop, K.: Primary weathering rates, water transit times, and concentrationdischarge relations: A theoretical analysis for the critical zone, Water Resour. Res., 53, 942–960, https://doi.org/10.1002/2016wr019448, 2017.
Andreadis, K. M., Schumann, G. J. P., and Pavelsky, T.: A simple global river bankfull width and depth database, Water Resour. Res., 49, 7164–7168, https://doi.org/10.1002/wrcr.20440, 2013.
Basu, N. B., Rao, P. S. C., Thompson, S. E., Loukinova, N. V., Donner, S. D., Ye, S., and Sivapalan, M.: Spatiotemporal averaging of instream solute removal dynamics, Water Resour. Res., 47, W00J06, https://doi.org/10.1029/2010wr010196, 2011.
Bergstrom, A., McGlynn, B., Mallard, J., and Covino, T.: Watershed structural influences on the distributions of stream network water and solute travel times under baseflow conditions, Hydrol. Process., 30, 2671–2685, https://doi.org/10.1002/hyp.10792, 2016.
Bertuzzo, E., Helton, A. M., Hall, R. O., and Battin, T. J.: Scaling of dissolved organic carbon removal in river networks, Adv. Water Resour., 110, 136–146, https://doi.org/10.1016/j.advwatres.2017.10.009, 2017.
Beusen, A. H. W., Bouwman, A. F., Van Beek, L. P. H., Mogollón, J. M., and Middelburg, J. J.: Global riverine N and P transport to ocean increased during the 20th century despite increased retention along the aquatic continuum, Biogeosciences, 13, 2441–2451, https://doi.org/10.5194/bg1324412016, 2016.
Bieroza, M. Z., Heathwaite, A. L., Bechmann, M., Kyllmar, K., and Jordan, P.: The concentrationdischarge slope as a tool for water quality management, Sci. Total Environ., 630, 738–749, https://doi.org/10.1016/j.scitotenv.2018.02.256, 2018.
Billen, G., Lancelot, C., and Meybeck, M.: N, P and Si retention along the aquatic continuum from land to ocean, in: Ocean Margin Processes in Global Change, 1st Edn., John Wiley & Sons, 19–44, 1991.
Birgand, F., Skaggs, R. W., Chescheir, G. M., and Gilliam, J. W.: Nitrogen removal in streams of agricultural catchments – A literature review, Crit. Rev. Env. Sci. Tec., 37, 381–487, https://doi.org/10.1080/10643380600966426, 2007.
Botter, G., Basso, S., RodriguezIturbe, I., and Rinaldo, A.: Resilience of river flow regimes, P. Natl. Acad. Sci. USA, 110, 12925–12930, https://doi.org/10.1073/pnas.1311920110, 2013.
Bouwman, A. F., Bierkens, M. F. P., Griffioen, J., Hefting, M. M., Middelburg, J. J., Middelkoop, H., and Slomp, C. P.: Nutrient dynamics, transfer and retention along the aquatic continuum from land to ocean: towards integration of ecological and biogeochemical models, Biogeosciences, 10, 1–22, https://doi.org/10.5194/bg1012013, 2013.
Boyer, E. W., Alexander, R. B., Parton, W. J., Li, C. S., ButterbachBahl, K., Donner, S. D., Skaggs, R. W., and Del Gross, S. J.: Modeling denitrification in terrestrial and aquatic ecosystems at regional scales, Ecol. Appl., 16, 2123–2142, https://doi.org/10.1890/10510761(2006)016[2123:Mditaa]2.0.Co;2, 2006.
Breiman, L., Friedman, J., Stone, C. J., and Olshen, R. A.: Classification and Regression Trees, Chapman & Hall, London, 1984.
Canfield, D. E., Glazer, A. N., and Falkowski, P. G.: The evolution and future of Earth's nitrogen cycle, Science, 330, 192–196, https://doi.org/10.1126/science.1186120, 2010.
Dehaspe, J.: R codes for a stream network model and a C–Q bending metric, HydroShare [code], https://doi.org/10.4211/hs.da70a09dc6074242ada756c29d12dcb3, 2021.
Diamond, J. S. and Cohen, M. J.: Complex patterns of catchment solutedischarge relationships for coastal plain rivers, Hydrol. Process., 32, 388–401, https://doi.org/10.1002/hyp.11424, 2018.
Dingman, S. L.: Analytical derivation of atastation hydraulic–geometry relations, J. Hydrol., 334, 17–27, https://doi.org/10.1016/j.jhydrol.2006.09.021, 2007.
Doyle, M. W.: Incorporating hydrologic variability into nutrient spiraling, J. Geophys. Res., 110, G01003, https://doi.org/10.1029/2005jg000015, 2005.
Dupas, R., Musolff, A., Jawitz, J. W., Rao, P. S. C., Jäger, C. G., Fleckenstein, J. H., Rode, M., and Borchardt, D.: Carbon and nutrient export regimes from headwater catchments to downstream reaches, Biogeosciences, 14, 4391–4407, https://doi.org/10.5194/bg1443912017, 2017.
Dupas, R., Minaudo, C., and Abbott, B. W.: Stability of spatial patterns in water chemistry across temperate ecoregions, Environ. Res. Lett., 14, 074015, https://doi.org/10.1088/17489326/ab24f4, 2019.
Durand, P., Breuer, L., Johnes, P. J., Billen, G., Butturini, A., Pinay, G., and Wright, R.: Nitrogen processes in aquatic ecosystems, in: The European Nitrogen Assessment: Sources, Effects and Policy Perspectives, edited by: Sutton, M., Howard, C., Erisman, J., Billen, G., Bleeker, A., Grennfelt, P., van Grinsven, H., and Grizzetti, B., Cambridge University Press, Cambridge, 2011.
Ebeling, P.: CCDB – catchment characteristics data base Germany, HydroShare [data set], https://doi.org/10.4211/hs.0fc1b5b1be4a475aacfd9545e72e6839, 2020a.
Ebeling, P.: WQQDB – water quality metrics for catchments across Germany, HydroShare [data set], https://doi.org/10.4211/hs.9b4deeca259b4f7398ce72121b4e2979, 2020b.
Ebeling, P., Kumar, R., Weber, M., Knoll, L., Fleckenstein, J. H., and Musolff, A.: Archetypes and Controls of Riverine Nutrient Export Across German Catchments, Water Resour. Res., 57, e2020WR028134, https://doi.org/10.1029/2020WR028134, 2021.
EEA – European Environmental Agency: DEM over Europe from the GMES RDA project (EUDEM, resolution 25 m) – version 1, October 2013, European Environmental Agency, EEA geospatial data catalogue, EEA [data set], https://sdi.eea.europa.eu/data/66fa7dca87724a5d9d562caba4ecd36a (last access: 14 March 2021), 2013.
Ehrhardt, S., Kumar, R., Fleckenstein, J. H., Attinger, S., and Musolff, A.: Trajectories of nitrate input and output in three nested catchments along a land use gradient, Hydrol. Earth Syst. Sci., 23, 3503–3524, https://doi.org/10.5194/hess2335032019, 2019.
Ensign, S. H. and Doyle, M. W.: Nutrient spiraling in streams and river networks, J. Geophys. Res.Biogeo., 111, G04009, https://doi.org/10.1029/2005jg000114, 2006.
Ensign, S. H., McMillan, S. K., Thompson, S. P., and Piehler, M. F.: Nitrogen and phosphorus attenuation within the stream network of a coastal, agricultural watershed, J. Environ. Qual., 35, 1237–1247, https://doi.org/10.2134/jeq2005.0341, 2006.
ESRI: ArcGIS Desktop: Release 10, Environmental Systems Research Institute, Redlands, CA, 2011.
Fisher, S. G., Sponseller, R. A., and Heffernan, J. B.: Horizons in stream biogeochemistry: Flowpaths to progress, Ecology, 85, 2369–2379, https://doi.org/10.1890/030244, 2004.
Galloway, J. N., Schlesinger, W. H., Levy, H., Michaels, A., and Schnoor, J. L.: Nitrogen fixation: Anthropogenic enhancementenvironmental response, Global. Biogeochem. Cy., 9, 235–252, https://doi.org/10.1029/95gb00158, 1995.
Godsey, S. E., Kirchner, J. W., and Clow, D. W.: Concentrationdischarge relationships reflect chemostatic characteristics of US catchments, Hydrol. Process., 23, 1844–1864, https://doi.org/10.1002/hyp.7315, 2009.
GomezVelez, J. D., Harvey, J., Cardenas, M. B., and Kiel, B.: Denitrification in the Mississippi River network controlled by flow through river bedforms, Nat. Geosci., 8, 941–945, https://doi.org/10.1038/Ngeo2567, 2015.
Hall, R. O., Tank, J. L., Sobota, D. J., and Mulholland, P. J.: Nitrate removal in stream ecosystems measured by ^{15}N addition experiments: Total uptake, Limnol. Oceanogr., 54, 653–665, https://doi.org/10.4319/lo.2009.54.3.0653, 2009.
Helton, A. M., Poole, G. C., Meyer, J. L., Wollheim, W. M., Peterson, B. J., Mulholland, P. J., Bernhardt, E. S., Stanford, J. A., Arango, C., Ashkenas, L. R., Cooper, L. W., Dodds, W. K., Gregory, S. V., Hall, R. O., Hamilton, S. K., Johnson, S. L., McDowell, W. H., Potter, J. D., Tank, J. L., Thomas, S. M., Valett, H. M., Webster, J. R., and Zeglin, L.: Thinking outside the channel: modeling nitrogen cycling in networked river ecosystems, Front. Ecol. Environ., 9, 229–238, https://doi.org/10.1890/080211, 2010.
Helton, A. M., Hall, R. O., and Bertuzzo, E.: How network structure can affect nitrogen removal by streams, Freshwater Biol., 63, 128–140, https://doi.org/10.1111/fwb.12990, 2018.
Hensley, R. T., Cohen, M. J., and Korhnak, L. V.: Inferring nitrogen removal in large rivers from highresolution longitudinal profiling, Limnol. Oceanogr., 59, 1152–1170, https://doi.org/10.4319/lo.2014.59.4.1152, 2014.
Horton, R. E.: Erosional Development of Streams and Their Drainage Basins – Hydrophysical Approach to Quantitative Morphology, Geol. Soc. Am. Bull., 56, 275–370, https://doi.org/10.1130/00167606(1945)56[275:edosat]2.0.co;2, 1945.
Jarvie, H. P., Sharpley, A. N., Kresse, T., Hays, P. D., Williams, R. J., King, S. M., and Berry, L. G.: Coupling HighFrequency Stream Metabolism and Nutrient Monitoring to Explore Biogeochemical Controls on Downstream Nitrate Delivery, Environ. Sci. Technol., 52, 13708–13717, https://doi.org/10.1021/acs.est.8b03074, 2018.
Jawitz, J. W. and Mitchell, J.: Temporal inequality in catchment discharge and solute export, Water Resour. Res., 47, W00J14, https://doi.org/10.1029/2010wr010197, 2011.
Kadlec, R. H. and Reddy, K. R.: Temperature effects in treatment wetlands, Water. Environ. Res., 73, 543–557, https://doi.org/10.2175/106143001x139614, 2001.
Kiel, B. A. and Cardenas, M. B.: Lateral hyporheic exchange throughout the Mississippi River network, Nat. Geosci., 7, 413–417, https://doi.org/10.1038/Ngeo2157, 2014.
Kirchner, J. W., Feng, X., and Neal, C.: Fractal stream chemistry and its implications for contaminant transport in catchments, Nature, 403, 524–527, https://doi.org/10.1038/35000537, 2000.
Klemes, V.: Dilettantism in Hydrology – Transition or Destiny, Water Resour. Res., 22, S177–S188, https://doi.org/10.1029/WR022i09Sp0177S, 1986.
Kuhn, M.: caret: Classification and Regression Training, R package version 6.086, available at: https://CRAN.Rproject.org/package=caret (last access: 16 December 2021), 2020.
Kumar, R., Hesse, F., Rao, P. S. C., Musolff, A., Jawitz, J. W., Sarrazin, F., Samaniego, L., Fleckenstein, J. H., Rakovec, O., Thober, S., and Attinger, S.: Strong hydroclimatic controls on vulnerability to subsurface nitrate contamination across Europe, Nat. Commun., 11, 6302, https://doi.org/10.1038/s41467020199558, 2020.
Kunz, J. V., Hensley, R., Brase, L., Borchardt, D., and Rode, M.: High frequency measurements of reach scale nitrogen uptake in a fourth order river with contrasting hydromorphology and variable water chemistry (Weisse Elster, Germany), Water Resour. Res., 53, 328–343, https://doi.org/10.1002/2016wr019355, 2017.
Langbein, W. B. and Leopold, L. B.: QuasiEquilibrium States in Channel Morphology, Am. J. Sci., 262, 782–794, https://doi.org/10.2475/ajs.262.6.782, 1964.
Leopold, L. B. and Maddock, T.: The hydraulic geometry of stream channels and some physiographic implications, US Government Printing Office, Washington, DC, 1953.
Li, L., Sullivan, P. L., Benettin, P., Cirpka, O. A., Bishop, K., Brantley, S. L., Knapp, J. L. A., Meerveld, I., Rinaldo, A., Seibert, J., Wen, H., and Kirchner, J. W.: Toward catchment hydrobiogeochemical theories, WIREs Water, 8, e1495, https://doi.org/10.1002/wat2.1495, 2020.
Lindgren, G. A. and Destouni, G.: Nitrogen loss rates in streams: Scaledependence and upscaling methodology, Geophys. Res. Lett., 31, 1–4, https://doi.org/10.1029/2004gl019996, 2004.
Marcé, R. and Armengol, J.: Modeling nutrient instream processes at the watershed scale using Nutrient Spiralling metrics, Hydrol. Earth Syst. Sci., 13, 953–967, https://doi.org/10.5194/hess139532009, 2009.
Marcé, R., von Schiller, D., Aguilera, R., Martí, E., and Bernal, S.: Contribution of Hydrologic Opportunity and Biogeochemical Reactivity to the Variability of Nutrient Retention in River Networks, Global. Biogeochem. Cy., 32, 376–388, https://doi.org/10.1002/2017gb005677, 2018.
Marinos, R. E., Van Meter, K. J., and Basu, N. B.: Is the River a Chemostat: Scale Versus Land Use Controls on Nitrate ConcentrationDischarge Dynamics in the Upper Mississippi River Basin, Geophys. Res. Lett., 47, e2020GL087051, https://doi.org/10.1029/2020gl087051, 2020.
McDonnell, J. J., Sivapalan, M., Vaché, K., Dunn, S., Grant, G., Haggerty, R., Hinz, C., Hooper, R., Kirchner, J., Roderick, M. L., Selker, J., and Weiler, M.: Moving beyond heterogeneity and process complexity: A new vision for watershed hydrology, Water Resour. Res., 43, W07301, https://doi.org/10.1029/2006wr005467, 2007.
Meybeck, M. and Moatar, F.: Daily variability of river concentrations and fluxes: indicators based on the segmentation of the rating curve, Hydrol. Process., 26, 1188–1207, https://doi.org/10.1002/hyp.8211, 2012.
Minaudo, C., Dupas, R., GascuelOdoux, C., Roubeix, V., Danis, P.A., and Moatar, F.: Seasonal and eventbased concentrationdischarge relationships to identify catchment controls on nutrient export regimes, Adv. Water Resour., 131, 103379, https://doi.org/10.1016/j.advwatres.2019.103379, 2019.
Mineau, M. M., Wollheim, W. M., and Stewart, R. J.: An index to characterize the spatial distribution of land use within watersheds and implications for river network nutrient removal and export, Geophys. Res. Lett., 42, 6688–6695, https://doi.org/10.1002/2015gl064965, 2015.
Moatar, F., Abbott, B. W., Minaudo, C., Curie, F., and Pinay, G.: Elemental properties, hydrology, and biology interact to shape concentrationdischarge curves for carbon, nutrients, sediment, and major ions, Water Resour. Res., 53, 1270–1287, https://doi.org/10.1002/2016wr019635, 2017.
Mueller, C., Musolff, A., Strachauer, U., Brauns, M., Tarasova, L., Merz, R., and Knöller, K.: Tomography of anthropogenic nitrate contribution along a mesoscale river, Sci. Total Environ., 615, 773–783, https://doi.org/10.1016/j.scitotenv.2017.09.297, 2018.
Mulholland, P. J. and Tank, J. L.: Can uptake length in streams be determined by nutrient addition experiments? Results from an interbiome comparison study, J. N. Am. Benthol. Soc., 21, 544–560, https://doi.org/10.2307/1468429, 2002.
Mulholland, P. J., Helton, A. M., Poole, G. C., Hall, R. O., Hamilton, S. K., Peterson, B. J., Tank, J. L., Ashkenas, L. R., Cooper, L. W., Dahm, C. N., Dodds, W. K., Findlay, S. E., Gregory, S. V., Grimm, N. B., Johnson, S. L., McDowell, W. H., Meyer, J. L., Valett, H. M., Webster, J. R., Arango, C. P., Beaulieu, J. J., Bernot, M. J., Burgin, A. J., Crenshaw, C. L., Johnson, L. T., Niederlehner, B. R., O'Brien, J. M., Potter, J. D., Sheibley, R. W., Sobota, D. J., and Thomas, S. M.: Stream denitrification across biomes and its response to anthropogenic nitrate loading, Nature, 452, 202–205, https://doi.org/10.1038/nature06686, 2008.
Musolff, A.: WQQDB – water quality and quantity data base Germany: metadata, Hydroshare, https://doi.org/10.4211/hs.a42addcbd59a466a9aa56472dfef8721, 2020.
Musolff, A., Fleckenstein, J. H., Rao, P. S. C., and Jawitz, J. W.: Emergent archetype patterns of coupled hydrologic and biogeochemical responses in catchments, Geophys. Res. Lett., 44, 4143–4151, https://doi.org/10.1002/2017gl072630, 2017.
Newbold, J. D., Elwood, J. W., Oneill, R. V., and Vanwinkle, W.: Measuring Nutrient Spiralling in Streams, Can. J. Fish. Aquat. Sci., 38, 860–863, https://doi.org/10.1139/f81114, 1981.
O'Brien, J. M., Dodds, W. K., Wilson, K. C., Murdock, J. N., and Eichmiller, J.: The saturation of N cycling in Central Plains streams: ^{15}N experiments across a broad gradient of nitrate concentrations, Biogeochemistry, 84, 31–49, https://doi.org/10.1007/s1053300790737, 2007.
Ocampo, C. J., Oldham, C. E., and Sivapalan, M.: Nitrate attenuation in agricultural catchments: Shifting balances between transport and reaction, Water Resour. Res., 42, W01408, https://doi.org/10.1029/2004wr003773, 2006.
Oldham, C. E., Farrow, D. E., and Peiffer, S.: A generalized Damköhler number for classifying material processing in hydrological systems, Hydrol. Earth Syst. Sci., 17, 1133–1148, https://doi.org/10.5194/hess1711332013, 2013.
Pennino, M. J., Kaushal, S. S., Beaulieu, J. J., Mayer, P. M., and Arango, C. P.: Effects of urban stream burial on nitrogen uptake and ecosystem metabolism: implications for watershed nitrogen and carbon fluxes, Biogeochemistry, 121, 247–269, https://doi.org/10.1007/s1053301499581, 2014.
Peterson, B. J., Wollheim, W. M., Mulholland, P. J., Webster, J. R., Meyer, J. L., Tank, J. L., Marti, E., Bowden, W. B., Valett, H. M., Hershey, A. E., McDowell, W. H., Dodds, W. K., Hamilton, S. K., Gregory, S., and Morrall, D. D.: Control of nitrogen export from watersheds by headwater streams, Science, 292, 86–90, https://doi.org/10.1126/science.1056874, 2001.
Pianosi, F. and Wagener, T.: A simple and efficient method for global sensitivity analysis based on cumulative distribution functions, Environ. Modell. Softw., 67, 1–11, https://doi.org/10.1016/j.envsoft.2015.01.004, 2015.
Pianosi, F. and Wagener, T.: Distributionbased sensitivity analysis from a generic inputoutput sample, Environ. Modell. Softw., 108, 197–207, https://doi.org/10.1016/j.envsoft.2018.07.019, 2018.
Pianosi, F., Sarrazin, F., and Wagener, T.: A Matlab toolbox for Global Sensitivity Analysis, Environ. Modell. Softw., 70, 80–85, https://doi.org/10.1016/j.envsoft.2015.04.009, 2015.
Pressley, A. N.: Elementary Differential Geometry, 1st Edn., Springer, London, https://doi.org/10.1007/9781848828919, 2001.
R Core Team: R: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, Austria, 2013.
RisseBuhl, U., Anlanger, C., Kalla, K., Neu, T. R., Noss, C., Lorke, A., and Weitere, M.: The role of hydrodynamics in shaping the composition and architecture of epilithic biofilms in fluvial ecosystems, Water Res., 127, 211–222, https://doi.org/10.1016/j.watres.2017.09.054, 2017.
Rode, M., Halbedel Nee Angelstein, S., Anis, M. R., Borchardt, D., and Weitere, M.: Continuous InStream Assimilatory Nitrate Uptake from HighFrequency Sensor Measurements, Environ. Sci. Technol., 50, 5685–5694, https://doi.org/10.1021/acs.est.6b00943, 2016.
Runkel, R. L. and Bencala, K. E.: Transport of reacting solutes in rivers and streams, in: Environmental Hydrology, Kluwer Academic Publishers, Dordrecht, the Netherlands, 1995.
Schlesinger, W. H., Reckhow, K. H., and Bernhardt, E. S.: Global change: The nitrogen cycle and rivers, Water Resour. Res., 42, https://doi.org/10.1029/2005wr004300, 2006.
Seitzinger, S. P., Styles, R. V., Boyer, E. W., Alexander, R. B., Billen, G., Howarth, R. W., Mayer, B., and Van Breemen, N.: Nitrogen retention in rivers: model development and application to watersheds in the northeastern U. S. A., Biogeochemistry, 57–58, 199–237, https://doi.org/10.1023/A:1015745629794, 2002.
Seybold, E. and McGlynn, B.: Hydrologic and biogeochemical drivers of dissolved organic carbon and nitrate uptake in a headwater stream network, Biogeochemistry, 138, 23–48, https://doi.org/10.1007/s1053301804261, 2018.
Stream Solute Workshop: Concepts and Methods for Assessing Solute Dynamics in Stream Ecosystems, J. N. Am. Benthol. Soc., 9, 95–119, https://doi.org/10.2307/1467445, 1990.
Tunqui Neira, J. M., Andréassian, V., Tallec, G., and Mouchel, J.M.: Technical note: A twosided affine power scaling relationship to represent the concentration–discharge relationship, Hydrol. Earth Syst. Sci., 24, 1823–1830, https://doi.org/10.5194/hess2418232020, 2020.
Vanni, M. J.: Nutrient Cycling by Animals in Freshwater Ecosystems, Annu. Rev. Ecol. Syst., 33, 341–370, https://doi.org/10.1146/annurev.ecolsys.33.010802.150519, 2002.
Vanni, M. J. and McIntyre, P. B.: Predicting nutrient excretion of aquatic animals with metabolic ecology and ecological stoichiometry: a global synthesis, Ecology, 97, 3460–3471, https://doi.org/10.1002/ecy.1582, 2016.
Wei, T. and Simko, V.: R package “corrplot”: Visualization of a Correlation Matrix, version 0.84, GitHub, https://github.com/taiyun/corrplot (last access: 16 December 2021), 2017.
Winter, C., Lutz, S. R., Musolff, A., Kumar, R., Weber, M., and Fleckenstein, J. H.: Disentangling the impact of catchment heterogeneity on nitrate export dynamics from event to longterm time scales, Water. Resour. Res., 57, e2020WR027992, https://doi.org/10.1002/essoar.10503228.1, 2021.
Wollheim, W. M., Vörösmarty, C. J., Peterson, B. J., Seitzinger, S. P., and Hopkinson, C. S.: Relationship between river size and nutrient removal, Geophys. Res. Lett., 33, L06410, https://doi.org/10.1029/2006gl025845, 2006.
Wollheim, W. M., Peterson, B. J., Thomas, S. M., Hopkinson, C. H., and Vörösmarty, C. J.: Dynamics of N removal over annual time periods in a suburban river network, J. Geophys. Res., 113, G03038, https://doi.org/10.1029/2007jg000660, 2008.
Wollheim, W. M., Mulukutla, G. K., Cook, C., and Carey, R. O.: Aquatic Nitrate Retention at River Network Scales Across Flow Conditions Determined Using Nested In Situ Sensors, Water Resour. Res., 53, 9740–9756, https://doi.org/10.1002/2017wr020644, 2017.
Wollheim, W. M., Bernal, S., Burns, D. A., Czuba, J. A., Driscoll, C. T., Hansen, A. T., Hensley, R. T., Hosen, J. D., Inamdar, S., Kaushal, S. S., Koenig, L. E., Lu, Y. H., Marzadri, A., Raymond, P. A., Scott, D., Stewart, R. J., Vidon, P. G., and Wohl, E.: River network saturation concept: factors influencing the balance of biogeochemical supply and demand of river networks, Biogeochemistry, 141, 503–521, https://doi.org/10.1007/s1053301804880, 2018.
Yang, X., Jomaa, S., Buttner, O., and Rode, M.: Autotrophic nitrate uptake in river networks: A modeling approach using continuous highfrequency data, Water. Res., 157, 258–268, https://doi.org/10.1016/j.watres.2019.02.059, 2019.
Yang, X. Q., Jomaa, S., Zink, M., Fleckenstein, J. H., Borchardt, D., and Rode, M.: A New Fully Distributed Model of Nitrate Transport and Removal at Catchment Scale, Water Resour. Res., 54, 5856–5877, https://doi.org/10.1029/2017wr022380, 2018.
Zarnetske, J. P., Haggerty, R., Wondzell, S. M., and Baker, M. A.: Dynamics of nitrate production and removal as a function of residence time in the hyporheic zone, J. Geophys. Res., 116, G01025, https://doi.org/10.1029/2010jg001356, 2011.