Bending of the concentration discharge relationship can inform about in-stream nitrate removal

. Nitrate (NO 3 − ) excess in rivers harms aquatic ecosystems and can induce detrimental algae growths in coastal areas. Riverine NO 3 − uptake is a crucial element of the catchment-scale nitrogen balance and can be measured at small spatiotemporal scales, while at the scale of entire river networks, uptake measurements are rarely available. Concurrent, low-frequency NO 3 − concentration and stream-ﬂow ( 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 ﬂow conditions (than under high ﬂows) are linked to biological NO 3 − uptake, creating a bent rather than linear log (C) –log (Q) relationship. Here we explore if network-scale NO 3 − uptake creates bent log (C) –log (Q) relationships and when in turn uptake can be quantiﬁed from observed low-frequency C – Q data. To this end we apply a parsimonious mass-balance-based 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 NO 3 − land-to-stream transfer. We ﬁnd 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 NO 3 − load removed in the network ( L r.perc ) but that network-wide ﬂow velocities should be taken into account when interpreting log (C) – log (Q) bending. Classiﬁcation trees, ﬁnally, 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 efﬁciently attenuate NO 3 − loads based on low-frequency NO 3 − and Q observations and generally show the importance of the channel geomorphology on the emerging log (C) –log (Q) bending at network scales.


Introduction
Transport and transformation of nitrate (NO 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 NO 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 NO 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 small-scale 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 NO 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 NO 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). NO 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 NO 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 Doyle, 2005;Marce and Armengol, 2009). Quantifying in situ NO 3 − uptake is labor-intensive and may involve local nutrient additions, potentially altering the ambient uptake rate (Hensley et al., 2014;Mulholland and Tank, 2002). Other methods require high-frequency 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 small-scale 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 NO 3 − removal capacity (Alexander et al., 2000;Horton, 1945), explained by high sediment-surface-to-watervolume 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 in-stream vegetation (Hensley et al., 2014;Rode et al., 2016) as well as for streams with a high capacity for lateral and hyporheic exchange (Gomez-Velez et al., 2015;Kiel and Cardenas, 2014). The scaling of in-stream uptake processes beyond the river reach has been approached by combined experimental-modeling 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 in-stream 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 long-term (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 in-stream 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 NO 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 low-flow and high-flow NO 3 − log(C)-log(Q) regression slopes for a ma-jority of the cases using low-frequency monitoring data. Moatar et al. (2017) found that stronger positive slopes under low-flow conditions correlate positively with chlorophyll a concentrations (associated with biological processes) and attributed this condition to biological NO 3 − concentration mediation in the stream. This is consistent with the findings of Hall et al. (2009) andHensley et al. (2014) among others that in-stream uptake is more efficient under low-flow than under high-flow conditions. Furthermore, Wollheim et al. (2017) illustrate non-linear NO 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 in-stream NO 3 − uptake is encoded within observed C-Q relationships, and their analysis therefore can improve our understanding of in-stream uptake processes by providing an alternative to elaborate fieldwork and modeling work aimed at quantifying NO 3 − removal in stream networks. Low-frequency NO 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 catchment-scale in-stream processing has yet to be investigated. In this paper, it is postulated that network-scale uptake effects can be inferred from the degree of non-linearity or the amount of bending of low-frequency, multi-annual 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., 2018Helton et al., , 2010and Mulholland et al., 2008) in 13 German catchments to explore the catchment-scale 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 in-stream biogeochemical parameters (e.g., channel shape, water velocity and biological NO 3 − uptake velocity), (iii) explore how C-Q bending is linked to network-scale in-stream uptake, and (iv) provide guidelines if and under what circumstances the C-Q bending can offer conclusive information on effective in-stream uptake. In this proof-of-concept exploratory study, we demonstrate how (existing) low-frequency monitoring data can be effectively utilized to quantify nitrate uptake in river networks and show how small-scale uptake processes shape emerging patterns at catchment scales.

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 "broken-stick" 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 m = arg max log Q |f |. 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 non-linearity as the amount of bending.
We assume here that a multi-annual (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)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) of French water quality stations with biweekly to monthly sampling frequencies . 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 intra-annual 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 low-frequency 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.

Network model
In this work an explorative grid-based (100 m × 100 m) mass balance network model (comparable to Bertuzzo et al., 2017;Helton et al., 2018Helton et al., , 2010and Mulholland et al., 2008; conceptually shown in Fig. B2) was used to simulate in-stream nitrate transport and biological removal on a daily basis. The model was developed in R (R Core Team, 2013).

Stream network and hydrological properties
Following Bertuzzo et al. (2017) and Helton et al. (2018), each river network node (i.e., grid cell) i (1 ≤ i ≤ N) has a drainage area A i [m 2 ] that is calculated as the sum of the total upstream drainage area j W j i 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 j i [-] is an element in the connectivity matrix W (N × N ) such that W j i = 1 if j is directly neighboring and flowing into i, and W j i = 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 time-varying 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, andBertuzzo et al., 2017) (Eq. 2). It is Q i , which in turn dictates the downstream and at-a-station hydraulic geometry relationships of river geomorphic parameters channel width, (w i ; m) and average channel depth (d i ; m) (Leopold and Maddock, 1953) where the flow length through a grid cell i, l i [m], equals 100 or 100 √ 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 [-] ∈ R + , which prescribes the cross-section geometry relation such that a triangular channel cross-section is represented by r = 1, a parabolic channel cross-section by r = 2, and channel crosssections with progressively flatter bottoms and steeper banks by increasing values of r (Dingman, 2007). The widthdischarge 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.

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 land-to-stream loading L in.ls,i [mg s −1 ], given that L = CQ (Eq. 3). The contribution of direct land-to-stream loading concentration can be expressed as a power law  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 NO 3 − loading is transport-limited rather than source-limited 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 long-term mean instream input concentration C mean [mg L −1 ] (Eq. A1). Additional NO 3 − sources such as the load resulting from NO 3 − release within the stream network and point sources are not considered here (similar to Bertuzzo et al., 2017, andWollheim 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 mixed-land-use 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 in-stream NO 3 − uptake follows first-order 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 · l i /Q i ) and biological (v f ) components (similar to Bertuzzo et al., 2017).
where P i is the wetted perimeter of the cross-section 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 NO 3 − molecules from the water column towards the biofilm at the pelagicbenthic interfaces and the sediments, where the in-stream 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 NO 3 − delivery (Doyle, 2005) (therefore it is not limited to denitrification only). We assume that v f is independent of the in-stream NO 3 − concentration C mean (Pennino et al., 2014;O'Brien et al., 2007) such that the areal uptake rate U = v f ·C mean [mg m −2 d −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 non-linearly for increasing C mean (10 −4 -10 1 mg L −1 ) when considering distinct catchments. However, in Germany, the NO 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 catchment-wide TT as the spatiotemporal median of the sum of all downstream T i (Eq. 2d) for a grid cell in the network ( Out i T i ) (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 first-order reaction constant

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 w +a d < 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 NO 3 − ; (ii) the network-wide 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  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.  Andreadis et al. (2013); Dingman (2007) a For v f , the selected range is an order of magnitude smaller than the one proposed by Marce et al. (2018) as we focus the analysis on the lower v f , where most of the bending happens (Sect. 3.3).

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 sub-catchments for the Selke as well as the Holtemme river system, both part of the Bode, a well-studied 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 tree-shaped river network with N grid cells or nodes.

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 in-stream processes (Rode et al., 2016;Dupas et al., 2017;Yang et al., 2019Yang et al., , 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 catchment-specific 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 landto-stream NO 3 − inputs averaged 1.2 kg N d −1 km −2 , which is similar to the 1.9 kg N d −1 km −2 reported by Winter et al. (2021) for the Selke river (Meisdorf), and it is well within the general 0.001 to 100 kg N d −1 km −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 (Risse-Buhl 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  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.

PAWN sensitivity analysis and correlation analysis
We performed a global sensitivity analysis (GSA) using the moment-independent PAWN method . 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 uninfluen- Table 2. Catchment properties summary: catchment area, median elevation, slope and topographical wetness index (TWI), maximum Strahler stream order, Horton form factor, drainage density, median discharge at the basin mouth over time (Q) with the corresponding runoff (area specific discharge) between brackets, and the coefficient of variation in the discharge in time (CV Q). The latter variable CV Q integrates the frequency of runoff events and the differences in recession constant (so the catchment's "flashiness" in response to rainfall) (Botter et al., 2013). . Germany DEM with the location and outline (shape) of selected catchments, along with their drainage networks (in blue) and outlet location (red triangle). See Table 2 for catchment IDs and properties. tial 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 = 1, . . ., 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|x i (y|x i ∈ I i,k )). Each conditional CDF is obtained by varying all in-puts 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):

ID
where KS( 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 non-influential 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 non-influential. 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 .
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 in-stream 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).

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 network-wide in-stream 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 non-linear synergistic interactions among model parameters and output variables. This non-parametric 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 catchment-wide 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 non-missing 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 Results and discussion

Empirical Curv max
The estimated Curv max for the French observed NO 3 − log(C)-log(Q) data  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 low-frequency C-Q time series is robust. The Spearman rank correlation (ρ = 0.53, p value < 2.2×10 −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 low-frequency time series does not impact the estimated Curv max . Overall, the proposed Curv max metric is suitable to quantify bending in multi-annual, temporally stable log(C)-log(Q) relationships.

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 NO 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, pbias = −0.4 %). 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 sum- Table 3. Classes containing low, medium and high values for response variables v f (uptake velocity), L r.perc (percentage of load removed) and Da (Damköhler number) are used for the CART training and testing samples. Similar classes are obtained for model output Curv max . These classes stem from distributing the non-missing simulation data over 20 percentiles and selecting the percentiles corresponding to low, medium and high literature values, with the respective percentile number (1-20) indicated in brackets. This was done to ensure that for one variable, each class contains the same number of simulation data points. The training sample for constructing the CART model was then allocated 80 % of this data and the test sample 20 %.

Variable
Units Low Medium High References mer and coincide with low simulated NO 3 − concentrations at the catchment outlet. The conservative NO 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 (< −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 mid-stream 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 first-order grid cells (average 24.1 kg N yr −1 ) that represent 55 % of the river network, followed by second-and third-order 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) in-creases with stream order from 1329 kg N yr −1 on average in a first-order grid cell to 5128 and 42 124 kg N yr −1 in secondorder and third-order 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, pbias = −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, andYang 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 third-order 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 NO 3 − concentrations (0.21 mg N L −1 ) and load (L = CQ), which are not represented in our model simulation (the lowest simulated NO 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   Table C1). Three representative grid cells covering low (A), intermediate (B) and high (C) total drainage areas show the incoming land-to-stream load as L in.ls (Eq. 3), the incoming load from upstream as L in.up (Eq. 3), the absolute removed load as L r (Eq. 5) and the outgoing load as L i (Eq. 4) in the log(L)-log(Q) space. The load removed as a percentage of the incoming load is presented on the secondary axis. Note that the corresponding Curv max values for these grid cells are calculated from the log(C)-log(Q) relationships rather than log(L)-log(Q). The insets show the distribution of Curv max and L r for each grid cell within a certain stream order. In Fig. 2a the observed and simulated concentrations are compared at point A.
for the entire simulated time period to better represent an effective long-term network-wide 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). Land-to-stream 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 non-linearly with discharge (Curv max = 0). The onset of a bent log(C)-log(Q) pattern (Curv max = −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., 2008Wollheim et al., , 2017Doyle, 2005;Basu et al., 2011). This decreased NO 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 surface-to-volume 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 land-to-stream 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 NO 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 lower-order to a higher-order 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 NO 3 − C-Q relationships that evolve from the decreasing removal efficiency at higher discharges.

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 land-to 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 land-to-stream 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 reaction-driven 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 transport-driven 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 catchment-specific conditions but rather to explore how NO 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 NO 3 − log(C)-log(Q) relationships in the French catchments (80 % of the values between −0.41 and −0.067; Fig. B4) . 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 Table 4. The 10th, 50th and 90th quantiles of model outputs Curv max , percentage of load removed, (L r.perc ) Damköhler number (Da), regression slope of the log(C)-log(Q) relationship at the basin outlet (b out ), the median concentration at the basin outlet (C out ), the median water velocity (v) and median river network travel times (TT) for each of the 13 German 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 reaction-driven 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 second-largest 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 vari-ability within the catchments (due to the effect of the chosen input parameter set).

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 sub-catchments 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 (ρ = −0.04) but shows a negative correlation with L r.perc (ρ = −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 (ρ = −0.28) such that higher bending coincides with more positive b out . The variable v is additionally strongly negatively correlated with L r.perc (ρ = −0.87), so a high percentage of load removed occurs at low velocities. Da on the other hand is pos-itively 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 (ρ = −0.71) than by TT (ρ = 0.48). Finally C out is negatively correlated with L r.perc (ρ = −0.82) and Da (ρ = −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 r,i ∼ 1/(Q 1−a w −a d ). 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 discharge-controlled transport processes, and no bending would be observed (Curv max = 0). Although increasing v f does correlate with decreasing concentrations (ρ = −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 network-wide velocity v (L r.perc and v, ρ = −0.82). This confirms that dis- Figure 4. The KS max sensitivity index for each of the model input parameters and each of the 13 simulated catchments. The input parameters related to the channel geometry (a d , a w , K d and K w ), land-to-stream loading (b and C mean ) and biogeochemistry (v f ) are shown together with two variables derived from some of the input parameters: the median velocity v and the Damköhler number Da. Each boxplot displays 15 000 bootstrapped estimates of KS max for each of the 13 simulated catchments. charge 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 % NO 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 (ρ = −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 transport-driven system with limited NO 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, ∼ −0.02) is thereby linked to low L r.perc (median 4.8 %), while low Curv max (more bending, ∼ −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 (Curv max < −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 high-bending, low-uptake 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 NO 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 Figure 5. Correlation matrix for the model parameter inputs: channel depth and width exponents (a d , a w ) and coefficients (K d , K w ); slope of the land-to-stream loading (b); concentration of the land-to-stream load (C mean ); uptake velocity (v f ); and the outputs bending of the log(C)-log(Q) relationship at the catchment outlet (Curv max ), effective stream velocity (v), first-order uptake constant (k), travel time (TT), Damköhler number (Da), daily percentage of load removed (L r.perc ), and slope of the log(C)-log(Q) relationship and median concentration at the outlet (b out and C out ). The Spearman rank correlation coefficients (ρ) are given for each combination. 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 (Curv max > −0.03) and high removal (L r.perc > 63.0 %) are even rarer (0.1 % of all simulations) compared to concurrent less bending (Curv max > −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 high-removal 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.

Predicting in-stream processing with Curv max
To determine if observed C-Q bending at the catchment outlet (here Curv max ) can be utilized to quantify in-stream 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 value < 2.2 × 10 −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 in-stream 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 land-tostream 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 high-resolution satellite pictures. The CART models could help obtain an initial probability of NO 3 − removal efficiency in a river network, especially in a context where network-wide 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.

Conclusions
In this study, we explore how low-frequency NO 3 − log(C)log(Q) relationships, observed at a basin outlet, can display bending as a result of network-scale in-stream uptake processes. We established a parsimonious grid-based river network model for 13 distinct German catchments and investigated the influence of in-stream loading, transport and uptake parameters on the bending of log(C)-log(Q) relationships. Based on our exploratory analysis we conclude the following.
-Noisy, multi-annual and low-frequency NO 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 multi-annual 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 land-tostream 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 NO 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, ρ = −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 geomor- Figure 7. CART decision trees for the response variables L r.perc (accuracy = 0.66), Da (accuracy = 0.51) and v f (accuracy = 0.40). The trees are read from top to bottom, where the binary splits are followed to arrive at a histogram, illustrating the probability of a test sample having low, medium or high values (Table 3) for a certain response variable. The variables at the binary splits differ per response variable and consist of the median stream velocity (v; m d −1 ) and Curv max for response variables L r.perc and v f , while width coefficient K w and depth exponent a d are identified additionally for Da. The prediction metrics for each of these classes and response variables are stated in Table C3.
phological 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 (ρ = −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 (Curv max < −0.51) and low removal (L r.perc < 5.2 %, 0.9 % of all simulations) or low bending (Curv max > −0.03) with high removal (L r.perc > 63.0 %, 0.1 % of all simulations) exist that are imposed by respective higher and lower network-wide 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 vari-ables 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 NO 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 low-frequency C and Q observations.

Appendix A
Calculation of c is performed as follows (Jawitz and Mitchell, 2011): with µ q = mean(log Q t.sp · a i ) µ c = log mean C − σ 2 c 2 σ c = b 2 · σ 2 q σ q = var(log Q t.sp · a i ).
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: Appendix B Figure B1. Conceptual figure explaining Curv max . The upper panel shows the log(C)-log(Q) relationship for Selke Meisdorf with the smoothed spline fits to these data for different degrees of freedom (df). The corresponding colors in the lower plot show then the local curvature values for these fitted smoothed splines. Also the log(Q) is indicated, for which the largest local curvature was found for each of the smoothed splines. Curv max is then calculated as the largest local curvature value for a degree of freedom of 5 in this specific case as it is the largest degree of freedom that has the largest local curvature within Q m ±0.05, with Q m equal to the log(Q) of the largest local curvature for df = 3. Note that when Curv max < 0, the curve is concave, while for Curv max > 0 the curve is convex.  In the "At-a-station" panel the changes in velocity, depth and width with Q for grid cell B (Fig. 3) in the middle of the Selke network are evaluated. The "Downstream" panel, which considers all the network grid cells, shows the channel characteristics width, depth and velocity for a time t, with a Q of 0.70 m 3 s −1 at the outlet. The values for the parameters a w , K w , a d and K d are the same for both scenarios (Table C1).    Table C1. Figure B7. Two specific cases of Fig. 6: (a) when Curv max is low (high bending) and L r.perc is low as opposed to high and (b) when Curv max is high (low bending) and L r.perc is high as opposed to low.
Data availability. The data we used were published by third parties and referenced in the text where appropriate, with complete access data in the reference list.