Hortonian Scaling of Coupled Hydrological and Biogeochemical Responses Across an Intensively Managed River Basin

Structural and functional attributes across fractal river networks have been characterized by well-established and consistent hierarchical, Hortonian scaling patterns. In most of the global river basins, spatial patterns of human settlements also conform to similar hierarchical scaling. However, emergent spatial hierarchical patterns and scaling of heterogeneous anthropogenic nutrient loads over a river basin are less known. As a case study, we examined here a large intensely managed river basin in Germany (Weser River; 46K km2; 8M population). Archived data for point-/diffuse-sources of total Phosphorus (Ptot) input loads were combined with numerical and analytical model simulations of coupled hydrological and biogeochemical processes for in-stream Ptot removal at the network scale. We find that Ptot input loads scale exponentially over stream-orders, with the larger scaling constant for point-source loads from urban agglomerations compared to those for diffuse-source contributions from agricultural and forested areas. These differences in scaling patterns result from hierarchical self-organization of human settlements, and the associated clustering of large-scale, altered land-cover. Fraction of Ptot loads removed through in-stream biogeochemical processes also manifests Hortonian scaling, consistent with predictions of an analytical model. Our analyses show that while smaller streams are more efficient in Ptot removal, in larger streams the magnitude of Ptot loads removed is higher. These trends are consistent with inverse scaling of nutrient removal rate constant with mean discharge, and downstream clustering of larger cumulative input loads. Analyses of six nested sub-basins within the Weser River Basin also reveal similar scaling patterns. Our findings are useful for projecting likely water-quality spatial patterns in similar river basins in Germany, and Central Europe. Extensions and generalizations require further examination of diverse basins with archetype spatial heterogeneities in anthropogenic pressures and hydroclimatic settings.


INTRODUCTION
Rivers are dendritic and directed drainage networks, where flows from all upstream contributing areas converge to the river basin outlet. Rivers draining heterogeneous landscapes exhibit universal fractal topography of self-organized, bifurcating networks (Rodríguez-Iturbe and Rinaldo, 2001;Paik and Kumar, 2008). Self-similar structures of river networks have been characterized mainly using two approaches. One is a consistent hierarchical scaling over Horton-Strahler stream-orders (Horton, 1945;Strahler, 1957). The hierarchical (Hortonian) scaling is embedded in the mean values of geomorphological or topologic metrics in drainage networks (Schumm, 1956;Peckham, 1995;Yang and Paik, 2017). The other is a power-law scaling between two geomorphological or geometric attributes in river basins or between one variable and its exceedance probability (Hack, 1957;Rodríguez-Iturbe et al., 1992;Dodds and Rothman, 2000).
Hortonian scaling relationships have been extensively found for diverse descriptors in complex river basin systems, by taking the advantages of simple but systematic and quantitative framework of stream-ordering. The discovery of the Hortonian scaling laws for hydrologic and channel hydraulic geometry variables in river networks allowed to model flow dynamics in streams of which practical datasets were not available (Mantilla et al., 2006;Gupta and Mesa, 2014;Gupta et al., 2015). In addition, spatial patterns of human-related attributes embedded in river basins were characterized through ordering scheme across diverse scales. Fang et al. (2018) presented preferential downward clustering pattern for self-organized human settlements in ∼70% large global river basins. Miyamoto et al. (2011) examined scaling of distinct human-managed land use over 109 large river basins in Japan. Yang et al. (2019a) analyzed archived data for urban settlements and the associated wastewater treatment plants (WWTPs) to examine spatial scaling in three large river basins in Germany. Such findings of selforganized, scale-invariant patterns in human-related activities imply that the associated scaling signatures in spatial patterns are likely to occur for nutrient loads generated from urban agglomerations as well as less populated areas (e.g., agricultural and semi-natural/forested areas).
Nonetheless, modeling biogeochemical dynamics coupled with functional indicators of hydrological conditions (e.g., discharge and stage) in river network pathways, to date, relies on the key assumption that nutrient inputs are spatially homogeneous or unstructured, i.e., loads per unit-area are spatially invariant (e.g., Bertuzzo et al., 2017;Helton et al., 2018). Our primary focus here is to investigate in heterogeneous drainage basins: (1) how the linkages between hydrological and biogeochemical filtering of input anthropogenic nutrient loadings contribute to scaling of corresponding response variables; and (2) if scaling patterns at the entire basin-scale are also evident at sub-basin scales. Given the observed streamorder scaling of hydrological conditions and urban settlements, we hypothesized that anthropogenic nutrient loads, as well as outcomes of biogeochemical processes would also scale across stream-order hierarchies. For the coupled processes, we hypothesized that non-linear biogeochemical buffering would dampen the Hortonian scaling patterns. To test these hypotheses, we simulated in-stream nutrient removal, using broadly-applied structural and functional scaling features over an entire river network, for total phosphorus (P tot ) which is the nutrient of our interest here. We explicitly included the effect of spatially structured point-source nutrient loads from WWTPs, and diffuse-source nutrient loads from distinct land covers. The largest national river basin in Germany (Weser River Basin, ∼46K km 2 ; 8M population) served as the case study of a humandominated river basin, with available archived data needed for model simulations of outcomes of the coupled hydrological and biogeochemical processes.

METHODS AND MATERIALS
Study Area: Weser River

General Information
The Weser River Basin, the largest national river basin, is enclosed entirely within Germany (∼46K km 2 , ∼10% of total area). The basin is located in the humid-temperate climatic zone, with mean precipitation of ∼780 mm/yr (600∼1,100 mm/yr) in excess of mean evapotranspiration of ∼505 mm/yr (400∼600 mm/yr), yielding a mean aridity index of ∼0.9 (IWRM-net, 2010; Hirt et al., 2012;Zink et al., 2017). According to the long-term (70 years : 1950-2019) daily monitoring records at the Intchede gauging station (BfG, 2019), located at the most downstream without tidal effect, annual mean volumetric discharge (Q), specific river discharge (q), and the coefficient of variation (CV) for both Q and q were: 317 m 3 /s, 265 mm/yr, and 0.73, respectively. Assuming spatially homogeneous hydroclimatic forcing, steady-state mean river discharge is estimated as the product between the given q and drainage area for each river streams and summed over the entire river network. Transient river flow conditions are outside of the scope of this paper. The nutrient uptake velocity (v f ) for P tot was estimated as 0.17 m/day at the whole basin-scale (Yang et al., 2019b).
The Weser River Basin was delineated using 100 m x 100 m grid DEMs (Zink et al., 2017) with a source area of 10 km 2 to form channels. The extracted river network consists of ∼1,700 stream segments with a total length of ∼12,000 km (here, a segment means a channel portion between two successive confluences). The mainstem (∼700-km) stretches from a small tributary in south-eastern part to the basin outlet with the final 7th order Horton-Strahler stream ( Figure 1A) (Horton, 1945;Strahler, 1957). Structure of the network manifests fractal geometry as commonly found for diverse rivers in the world (Yang et al., 2019a).

Anthropogenic Nutrient Loads
About 8.4 million people (∼10.2% of Germany's population) reside in the Weser River Basin (Figure 1A), with different sizes of urban agglomerations, with ∼97% connection to WWTPs (EEA, 2017). Out of total ∼8,900 WWTPs in Germany (Büttner, 2020), ∼845 WWTPs are located in the Weser River Basin (Figure 1B). Five class-sizes k are assigned to WWTPs based on the population equivalent (PE) served (https:// www.gesetze-iminternet.de/abwv/anhang_1.html) (EEC, 1991). Smaller WWTPs (k ≤ 3) are clustered toward upstream locations serving small communities, while larger WWTPs (k ≥ 4) are located downstream serving larger cities clustered near the river basin outlet (Yang et al., 2019a). Tertiary treatment technology is mandatory in the two largest class-sizes (k ≥ 4) to meet stringent environmental regulations for nutrient loads discharged to the river network, while a suite of primary and secondary treatment is deployed at the smaller WWTPs (k ≤ 3). Sum of P tot loads discharged from ∼845 WWTPs to the Weser River is ∼957 tons P/yr loads (for the period 2012-2016). Around 71% of the P tot loads from point-sources (∼680 tons P/yr) are generated by ∼260 large WWTPs (k ≥ 4), representing about a third of the total WWTPs located in the Weser River Basin. P tot loads from diffuse-sources to river streams are categorized based on distinct land-covers delineated from the CORINE Land Cover (CLC) 100-m grid data for the year 2012 (https://land. copernicus.eu/pan-european/corine-land-cover). Five main classes in the Level 1 of the CLC data are used to characterize different land-covers. In the Weser River Basin (Figure 1C), agricultural areas account for the largest area (CLC-class-2, ∼59%), followed by forested/semi-natural areas (CLC-class-3, ∼31%), and urban surface areas (CLC-class-1, ∼9%). These three CLC-classes covering ∼99% of total area were considered in estimating total P tot loads from diffuse-sources to streams. Respective total P tot input loads contributions from CLC-classes for the year 2010 were calculated as: 419 tons P/yr for CLCclass-1, 1,576 tons P/yr for CLC-class-2, and 485 tons P/yr for CLC-class-3 (UBA, 2010(UBA, , 2017. In turn, P tot input loads per unit-area Ψ [kg P/yr/km 2 ] are given as Ψ UB =101, Ψ AG =58, and Ψ FO =34 for CLC-class-1 to 3, respectively. The estimated total P tot loads from all diffuse-sources is 2,481 tons P/yr, which is 2.6 times larger than those provided by point-sources WWTPs. Table 1 summarizes information on anthropogenic P tot input loads for the entire Weser River Basin, and for two 6th-order sub-basins and four 5th-order sub-basins along with corresponding drainage areas and basin shape factors.

Modeling Elemental Removal at the Network Scale
Models for element removal at the scale of a whole river network have gained attraction in the recent literature (Basu et al., 2011;Wollheim et al., 2015;Bertuzzo et al., 2017). In particular, Bertuzzo et al. (2017) investigated scaling patterns in the dynamics of terrestrial dissolved organic carbon removal assuming a first-order losses in each stream segment. The generic model framework considers spatially invariant nutrient loads (mass/area) from only diffuse-sources to river reaches. Here, we additionally included spatially heterogeneous point-and diffuse-source nutrient loads, and implemented the model for the Weser River Basin, by relying on available datasets for parameter estimation. With the relevant modification of the original model, core governing equations and their underlying processes are described below.
For a given reach i in a river network, incoming P tot loads, L in,i , to the receiving reach are a summation of three parts: (1) point-source P tot loads, L PS,i , discharged from WWTPs; (2) diffuse-source P tot loads, L DS,i , delivered by directly drained areas of different land-covers; and (3) P tot loads transported from direct upstream reaches j, L out,j (i.e., L out,j = 0 for headwaters): a The symbol of Ω indicates the final stream-order for a given river basin. Sub-basins with the same Ω are independent within the entire boundary of the Weser River Basin. b We quantify elongation of a basin shape using a unit shape factor, Fu = [L T,i /(A T,i ) 0.5 ], where L T,i is the length of the main stem in a given basin i, and A T,i is the total drainage area of the basin i. Higher Fu means more elongated shape of a river basin.
We assume temporally constant P tot loads from both point-and diffuse-sources to receiving river reaches, to be compatible with available dataset for mean annual loads. We acknowledge that intra-annual variabilities of nutrient loads are neglected in such an analysis. Outflowing P tot loads at the most downstream point of each reach i L out,i is estimated as the residual P tot loads after in-stream natural attenuation, described as a first-order decay processes during residence time τ i : The net uptake rate scales inversely with the water depth d i as v f ,i /d i , where v f ,i is the uptake velocity for nutrient loss in the hyporheic and benthic zone (Wollheim et al., 2006;Botter et al., 2010;Basu et al., 2011;Yang et al., 2019b). The reachscale residence time τ i is straightforwardly estimated through the geomorphological and hydrological conditions at the reach i : where w i and l i indicate the width and length of the reach, respectively. Here, spatially homogeneous v f ,i is assumed over the entire river network; compared to order of magnitude changes in discharge across stream-orders, biogeochemical and bathymetry variations exert minor influence (e.g., Basu et al., 2011). The fraction of P tot loads removed in each reach f R,i represents the nutrient removal efficiency, as: where the P tot loads removed in each reach, L r,i , is equivalent to (L in,i -L out,i ). For spatially homogeneous hydroclimatic conditions over the entire river network, represented by mean precipitation and evapotranspiration, the steady-state river discharge, Q i , at a reach is proportional to the contributing drainage area A i : where q is the specific river discharge [L/T] indicating net precipitation. Based on well-known scaling relationships of channel geometry (Leopold and Maddock, 1953), we calculate the width w i and water depth d i for each reach as a power-law scaling of Q i for the same frequency of occurrence, as: where the power-law exponents α w = 0.5 and α d = 0.4 are suitable representatives across many rivers and so are the coefficients K w = 10 [m 1−3α w s α w ] and K d = 0.25 [m 1−3α d s α d ] (Leopold and Maddock, 1953;Knighton, 1998). Anthropogenic alterations to channels and flow diversions are neglected in this analysis (see later discussion in section Limitations and Extensions).

Scaling Analyses
Hortonian scaling was used to identify spatial patterns for several coupled processes and responses across hydrological, biogeochemical, and anthropogenic domains over a river basin. Specific attributes examined in this study are: (1) Diffuse-sources P tot loads from three land-cover types to rivers, i.e., urban surface areas [L UR(P) ], agricultural areas [L AG(P) ], and forested/seminatural areas [L FO(P) ]; (2) Summed P tot loads from both pointand diffuse-sources [L Total(P) ]; and (3) Fraction and magnitude of P tot loads removed (f R and L r , respectively) through in-stream biogeochemical processes.
The Hortonian scaling in river networks represents the scaleinvariance in spatial hierarchy of averaged attributes over streamorders (Horton, 1945;Strahler, 1957). For a given attribute Z, the Hortonian scaling is demonstrated when its mean <Z> varies exponentially over stream-order ω with an exponent h: where h > 0 indicates a scale-free increase over ω, such as for drainage area (Schumm, 1956) and human settlements (Fang et al., 2018;Yang et al., 2019a), while h < 0 denotes a scaleinvariant decrease over ω, as for the number of streams (Horton, 1945). The consistent rate of change for <Z> over increasing ω is quantified as a Hortonian scaling ratio (HSR) for the attribute Z, R Z > 0, as:

Stream-Order Model
Given the Hortonian scaling in each source of P tot loads, we adopted the stream-order model of Miyamoto et al. (2011) to estimate the HSR for total P tot input loads as a multiplicative function: where Σδ x = 1 (x = PS, UR, AG, and FO). Each power-law exponent δ represents the relative contribution to the Hortonian scaling of total P tot input loads from a specific source. That is, δ x > 0 indicates clustering toward downstream, while δ x < 0 suggests upstream clustering.

Hortonian Scaling for Anthropogenic P tot Loads
Total of P tot loads, from both point-and diffuse-sources in the Weser River Basin, followed structural hierarchical scaling over stream-orders ω, with R L_TOT = 2 (R 2 > 0.9, p < 0.05) (Figure 2A). P tot loads from each of point-and diffuse-sources also scaled with ω, but each with distinct rate of change (Figures 2B-E). The largest HSR value was for point-source (R L_PS = 2.6; R 2 > 0.9, p < 0.05) followed by diffuse-sources from urban surface areas (R L_UR = 1.9; R 2 > 0.9, p < 0.05), agricultural areas (R L_AG = 1.6; R 2 ≈ 0.7, p < 0.05), and forested/semi-natural areas (R L_FO = 1.4; R 2 ≈ 0.2, p < 0.2). The decreasing trend for HSRs from more human-dense to less habitable areas was consistent in six sub-basins embedded in the Weser River Basin ( Table 2). The values of HSR for a given source of P tot loads were within a narrow range among independent sub-basins. For example, four 5th-order sub-basins had R L_PS = 3.1 ± 0.7, R L_UR = 2.4 ± 0.3, R L_AG = 1.9 ± 0.2, and R L_FO = 1.7 ± 0.3 (mean ± std.; where std. denotes the standard deviation), highlighting consistency in the scaling patterns despite geographical and environmental differences among the sub-basins.
Employing the stream-order model, unknown values for four δ x were estimated based on the empirical data for P tot input loads over six sub-basins and the Weser River Basin ( Table 2). The resultant δ x set was given as [δ PS , δ UR , δ AG , δ FO ] = [0.24, 0.60, 0.27, −0.20], which yielded accurate predictions of R L_TOT compared to empirically derived R L_TOT (difference < 3%) (Figure 3). The estimated values of δ x were aligned with the distinct patterns in the spatial clustering of different land covers over the river basin. The value of δ x < 0 for the P tot loads from forested/semi-natural areas reflected their upstream clustering patterns, whereas that of δ x > 0 for urban and agricultural P tot loads mirrored their downstream-skewed distributions (see Figure 1C).

Fraction of P tot Loads Removed
The fraction of removed P tot loads (f R ) within a reach, compared to incoming loads, reflected the inverse relation of the instream nutrient uptake rate and river discharge (Figure 4A). Across the river basin, f R ranged over two orders of magnitude, representing ∼56% in-stream loss of total P tot input loads over the entire river network under the annual mean river discharge condition. Spatial distribution of f R was structured by streamorders (Figure 4B), decreasing with higher ω, with R fR = 1.9 (R 2 > 0.9, p < 0.05).
We also analytically examined the Hortonian scaling of f R, since the term of f R is expressed indeed as a power function of Q for suitably short reach-unit (see Appendix A). The analytical expression in Equation (A5) presents that Hortonian scaling of f R can be estimated with the exponent of h fR = (α w -1)h Q , if hierarchical scaling of Q, h Q , and hydraulic geometry scaling of stream width, α w , are given for a river basin of interest. We expect R fR values to range between 1.1 and 4.2, estimated using extensively-found literature values for h Q = 1.1 ∼ 1.8 (Schumm, 1956;Rodríguez-Iturbe and Rinaldo, 2001) and α w = 0.2 ∼ 0.9 (Klein, 1981;Knighton, 1998). Analytically derived R fR ∼1.9 for the entire river basin was almost identical with statistically estimated R fR based on simulation results (difference = 0.8%). For the six sub-basins, there were no discrepancies between the R fR results from analytical and statistical approaches (difference < 5%). The R fR values were within its feasible range based on the literature. Compared to the Hortonian scaling for hydrologic condition R Q , the values of R fR declined incrementally for the entire basin and its sub-basins ( Table 2).
In addition to the reach-scale analysis, the fraction of P tot loads removal estimated at the entire river network revealed a scaling pattern over the Weser River Basin. The fraction of total P tot loads removed in streams (F R ) with the drainage areas smaller than a threshold of interest A increased non-linearly with increasing A, reaching ∼56% removal at the basin outlet with maximum A (Figure 4C). The general pattern of F R over A from our simulation results agreed well with the analytical equation derived by Bertuzzo et al. (2017) where α is a scaling exponent of the power-law relation between nutrient input loads and contributing area for each reach, and A T is a source area to define channel initiation). Note that, by following the original framework, we estimated the values of the scaling exponent α from only diffuse-sources P tot loads, resulting in α = 0.91. Underestimation from using the original analytical equation was attributed to not-including point-sources P tot loads.

Magnitude of P tot Loads Removed
The removed P tot loads in each reach (L r ) increased over stream-orders (ω) (Figure 4D). The distribution of L r for the entire river basin was hierarchically structured FIGURE 2 | Stream-order based distributions of P tot loads for the Weser River Basin. (A) Total P tot loads from both point-and diffuse-sources L Total(P) . (B) Reported P tot loads discharged from WWTPs L PS(P) . (C-E) Estimated P tot loads delivered from urban surfaces L UR(P) , agricultural areas L AG(P) , and forested/semi-natural areas L FO(P) , respectively. For a given stream-order ω, each gray dot indicates the estimated P tot loads from directly draining areas into a stream of the ω-th order. Mean of gray dots within the same ω are depicted as red bars, and their linear trends over ω are fitted on the semi-log plot. Hortonian scaling equations with exponent h, and the resultant Hortonian scaling ratios are embedded in each plot. A small letter and a number in the right upper parenthesis indicate a range of the p-value and R 2 value for each regression line fitted to estimate Hortonian scaling trends. For the p-value, three categories are given as: no mark for p < 0.1, (a) for 0.05 ≤ p < 0.1, and (b) p ≈ 0.2. For the R 2 value, five classifications are denoted as: no mark for R 2 > 0.9, (1) for 0.8 < R 2 ≤ 0.9, (2) for 0.7 < R 2 ≤ 0.8, (3) for 0.6 < R 2 ≤ 0.7, and (4) for R 2 ≤ 0.6. on average, resulting R Lr = 1.8 (R 2 > 0.9, p < 0.05). The Hortonian scaling of L r was also found for the six sub-basins with a narrow range of 1.9 ± 0.1 (mean ± std.; Table 2). The increasing trend of L r over ω in the Weser River Basin reflected the spatial distribution of P tot input loads from both point-and diffuse-sources, and their downstream aggregation with converging flows. Hence, despite a much less efficient removal of P tot loads (i.e., lower f R ) in larger streams, the greater magnitude of P tot input loads yielded the positive relation between L r and ω.
Length-normalized L r , i.e., the removed P tot loads per unit length ζ R , was lower in smaller drainage areas A and increased for larger A (Figure 4E). The overall trend of our model simulation was aligned with the analytical expression derived by Bertuzzo et al. (2017) where K Y is the nutrient input loads per unit area, and each of K w and α w are the coefficient and exponent in the power-law function of the hydraulic geometry for river width). As presented in Figure 4C, the analytical scaling here also showed underestimation for the simulated results, which is consistent with the original model assumption that neglected point-sources. Nevertheless, FIGURE 3 | Comparison of the Hortonian scaling ratio for total P tot loads input from point-and diffuse-sources, R L_TOT . Data-driven R L_TOT on the x-axis was originated from the data analyses, as presented in Figure 2A. Estimated R L_TOT on the y-axis was the result of the multiplicative model for R L_TOT as a power function of four Hortonian scaling ratios for P tot input loads from individual point-and diffuse-sources. In addition to the Weser River Basin (blue square), its six sub-basins with its corresponding final stream-order Ω = 5 or 6 (orange circle; purple triangle, respectively) are marked with the acronym of each name, as shown in Figure 1A.
for the first time our study corroborated that the analytical scaling derived from model simulations reproduces the scaling of nutrient load removal in a large river basin characterized by strong spatial heterogeneity.

Hydrological and Biogeochemical Filtering
Our model framework can be described using the concept of filtering for hydrological and biogeochemical processes. Hydrological responses observed in river networks can be interpreted as the outcomes of sequential functioning of drainage area acting as a "filter" of external forcing (Gall et al., 2013;Kim et al., 2016;Kumar et al., 2020). Following the underlying mass balance principle, hydroclimatic forcing (precipitation, PR) filtered by evapotranspiration (ET) losses generates mean river discharge (Q) through linear filtering of contributing drainage area (A), resulting Q = A(PR-ET). If both PR and ET are assumed spatially homogeneous over a river basin, as assumed in this study, Q linearly varies with A which is known to scale exponentially with stream-orders, i.e., Hortonian scaling (Schumm, 1956). Hence, spatial heterogeneity of Q is dependent only on heterogeneity of A across stream orders in the river basin. Furthermore, well-known scaling relations of mean Q with water depth, width, and flow velocity along downstream channels indicate the hierarchical patterns in spatial heterogeneity of the hydraulic geometry variables (Leopold and Maddock, 1953;Paik and Kumar, 2004;Gupta and Mesa, 2014).
Biogeochemical responses in river networks can also be interpreted using the "filter" concept just discussed for hydrological responses. The filters might be either natural or engineered attenuation processes, or often both . Nutrient loads applied to the landscape are buffered by crop/tree uptake, and further attenuated by biogeochemical filtering in the vadose and saturated zones, before entering the stream networks as diffuse-sources (Mockler et al., 2017;Liu et al., 2018). Similarly, nutrient loads from urban agglomerations are aggregated and treated at the WWTPs and discharged to river networks as point-sources (Aissa-Grouz et al., 2018;Yang et al., 2019a). Note that the multiple steps of filtering from original sources, such as landscape surface or households, to receiving streams can be flexibly implemented in a simulation model depending on data availability of nutrient loads input and temporal scale of interests for a given study. For example, our model did not need to consider transport pathways of nutrient loads from original sources to streams because of reliable estimated data on nutrient loads directly reaching streams (see section Anthropogenic Nutrient Loads). Once nutrient loads arrive in rivers and are transported along channel flow paths, they are removed through in-stream processing represented as a first-order nutrient uptake process, which scales inversely with stream water depth, which in turn scales with Q (Wollheim et al., 2006;Botter et al., 2010;Basu et al., 2011). Theoretical analysis of Bertuzzo et al. (2017) presents that, at the scale of an entire river network, the cumulative in-stream losses of dissolved organic carbon through the nonlinear in-stream filtering scale with contributing area, A (note: Q scales with A). This suggests that our model framework for phosphorus can be directly extended to nitrogen and carbon attenuation, assuming the first-order rate processes for in-stream removal (Wollheim et al., 2015;Bertuzzo et al., 2017;Helton et al., 2018).
Indeed, our finding of consistent R fR < R Q in the Weser Basin represents that additional non-linear filtering by biogeochemical (in-stream nutrient loss) processes at the scale of river network dampens the variation based on the "mean" values for a given attribute across stream-orders. Our analytical derivation for R fR demonstrates that the damping effect from the coupling complex dynamics is indicated as the multiplication of (α w -1) in the Hortonian scaling exponent of f R . Given the universal features for the Hortonian scaling of Q and the hydraulic geometry relations, the analytical estimation for R fR should be applicable to other river basins, which needs to be examined in future research. A more detailed discussion on R fR is reported below (section Scaling of Input Loads and In-stream Removal).

Scaling of Input Loads and In-stream Removal
Our analyses show that the mean values of nutrient P tot input loads exhibit consistent Hortonian scaling, with larger scaling constant for point-source loads delivered from WWTPs, relative to that for diffuse-sources contribution from agricultural and forested/semi-natural areas (R L_PS ∼ 2.6 vs. R L_AG ∼ 1.6; R L_FO ∼ 1.4). This trend is also noted in the six sub-basins: two 6thorder and four 5th-order, of the Weser River. The clustered  (12) and (10) derived in Bertuzzo et al. (2017), respectively. Following the original framework, the exponent α = 0.91 was fitted from scaling of input mass flux only from diffuse-sources over drainage areas. The other exponent α w = 0.5 came from the hydraulic geometry scaling between width and river discharge. spatial patterns result from self-organized hierarchy of human settlements and the associated different sizes of WWTPs, and human-altered land-covers. These results are consistent with more downstream clustered spatial organization of larger classsizes of WWTPs which discharge P tot loads from larger urban aggregations of populations (Yang et al., 2019a,b). These urban settlements account only for <10% of the basin area, and are embedded largely within agricultural areas (∼59%), while the forested/semi-natural areas (∼31%) are clustered upstream.
Our findings indicate how scalings of these P tot input loads are filtered by coupled hydrological (river discharge) conditions and biogeochemical (in-stream nutrient loss) processes along converging flows in the river network. We identify that the fraction of P tot loads removed (f R = 1 -L out /L in ) scales inversely with stream-orders, with a scaling constant of ∼0.66, reflecting the inverse relationship of removal rate constant and stream discharge (or stage). We also demonstrate close agreement between empirical and analytical estimation of the Hortonian scaling for f R (R fR ). However, the magnitude of P tot loads removed (L R ) increases over stream-orders with a scaling exponent of ∼0.6, mirroring the trend of greater P tot input loads toward downstream over the Weser River Basin.
Our analyses show that while smaller streams (ω ≤ 3) are more efficient in nutrient removal, larger total P tot loads are removed in the reaches of higher-order streams (ω ≥ 5). Our numerical simulation results show that scaling patterns over drainage area for the cumulative fraction of P tot loads removed (F R ) and the P tot loads removed per unit stream length (ζ r ) are well-approximated by the analytical expressions, originally derived by Bertuzzo et al. (2017), for spatially homogeneous, diffuse-source loads proportional only to contribution areas. Thus, the analytical expressions for R fR and ζ r are useful for estimating the hierarchical gradient in the fraction of in-stream nutrient removed and the magnitude of in-stream nutrient losses, respectively, without exact simulation of biogeochemical processes for other catchments.
While the Hortonian scaling based on mean values is indeed interesting, we also recognize an important signature for greater variability manifested in lower-order streams draining highly diverse and small sub-catchments. The low-order streams with small contributing areas are biogeochemically most-reactive and are also critically important for preserving aquatic ecological diversity (Wohl, 2017;Richardson, 2019). Hence, they require greater attention in catchment monitoring and management efforts (Biggs et al., 2017;Jensen et al., 2019). With converging flows within the dendritic drainage network, hydrological and biogeochemical attributes in larger streams represent flowweighted spatial averaging, thereby contributing to decreased variability. A future study on the variability trends is expected to provide constructive insights on how to synthetically interpret the Hortonian scaling and its corresponding variabilities over the entire river networks.

Comparison to Monitored P Concentration
We examined the scaling of modeled P tot and monitored soluble reactive phosphorus concentrations (Supplementary Figures 1,  2). Although the median values of the modeled and monitored concentrations were different in smaller streams (ω ≤ 3) [at the 1% significance level (Kruskal-Wallis test)], both datasets exhibited weak scaling over ω as supported by larger streams (ω ≥ 4) without significant differences in their distributions. The weak scaling patterns reflect an increased filtering (i.e., signal dampening) from the combined influence of local and cumulative hydrological processes (e.g., dilution and flow aggregation) and in-stream biogeochemical attenuation resulting the residual nutrient loads in a given reach (Supplementary Data).

Limitations and Extensions
Our analyses here are based on mean-annual discharge, and constant mean-annual input P tot loads from diverse sources for one study river basin; as such, they are most useful for long-term projections. We acknowledge that inter-and intraannual dynamics in discharge (and stage) and P tot loads will need to be considered to examine how temporal variability modifies the scaling patterns of mean values we presented herein. Probability density function of Q and its statistical moments (e.g., mean, median, mode, variance, etc.) are useful for characterizing stochastic temporal variabilities in Q. Botter et al. (2013) suggested that hydrological regimes of rivers can be divided into two broad groups: persistent, with low variability in discharge (CV Q < 1), and erratic, with high variability (CV Q > 1). Doyle (2005) proposed the use of a "functionally equivalent" river discharge (Q fed ) for reproducing the effect of in-stream nutrient removal over the full spectrum of temporal variability in Q. His analysis showed that Q fed was equivalent to Q mode for low variability in Q (i.e., persistent hydrological regimes), while their similarity became weaker as the Q-variability increased. Other fixed values of constant Q are also useful for bracketing likely range in hydrological conditions. For example, the 10th and 90th percentile Q can represent, respectively, the low-flow conditions during droughts, and high-flow conditions during floods (Rice and Westerhoff, 2017;Thober et al., 2018). We note that the scaling constant for the fraction P tot loads removed (f R ) will be identical for any temporally constant Q at the same frequency of occurrence; only the intercept will increase (decrease) with decreasing (increasing) mean values of Q over stream-orders. This topic on temporal variabilities in Q and P tot loads over a wide range of climate and land cover characters should be addressed in a future study.
We used the fixed values for the coefficients and exponents of power-law scaling in the hydraulic geometry relations (α w and α d in Equations 6, 7) based on global-scale findings from established literature. Our model framework neglected possibly modified hydrological flow regimes due to anthropogenic alterations to the river network, to focus on our primary aim at identifying the hierarchical scaling characters in biogeochemical responses over heterogeneous landscapes. We do acknowledge that hydraulic structures, such as dams and their operations, do alter river flow regimes (e.g., Williams and Wolman, 1984;Ferrazzi and Botter, 2019). In addition, over the long term, channel modifications also change hydraulic geometry variables (e.g., Shin and Julien, 2010;Ran et al., 2012) over a wide range depending on a river and its floodplain (e.g., Knighton, 1998). In the Weser River Basin, primarily flood-control dams are located at higher altitudes in upstream locations (Speckhann et al., 2021), and their effect is indeed captured in the mean discharge we used on the basis of measured discharge at the basin outlet (Supplementary Data; Supplementary Figure 1). Site specific studies in river basins with significant alterations of river flow regimes and modifications of channels should be conducted in the future, to better understand how the Hortonian scaling patterns found here would be different for the human-altered hydrological and geomorphological attributes.
Despite the limitations mentioned above, our findings can be used to predict water-quality impairments in two other large basins in Germany (Rhine and Elbe), which have similar spatial patterns and scaling features (Yang et al., 2019a). Extensions to other basins would require characterizing the spatial patterns of urban agglomerations with high connectivity to WWTPs, and changes of land-cover and land-use. Spatial patterns for nutrient loads from point-sources WWTPs reflect aggregated wastewater from households served by each WWTP, and centralized treatment which reduces significant nutrient loads before discharging them to streams, in order to comply with environmental regulations. Thus, spatial patterns of urban agglomerations, connectivity to centralized WWTPs, variations in treatment technologies adopted, and enforcement of regulations also influence the scaling patterns of point-sources loads in other river basins. Compared to German basins with multiple small/medium size cities, river basins with mega-cities near the outlet may be characterized by a larger Hortonian scaling exponent. In such basins, agricultural or forested landscapes, and diffuse-sources nutrient loads, might inversely mirror human settlement areas, thus exhibiting weaker scaling than suggested by our analyses here. Extending our analyses to explore scaling patterns in diverse managed river basins, however, remains a future challenge for developing a generalized understanding of the scaling for nutrient input loads and biogeochemical responses in large river networks across wider climate and hydrological regimes.

CONCLUSIONS
The Hortonian scaling is a classical signature representing the consistent hierarchical distribution of natural or anthropogenic attributes over a whole river network. In this study, we examined whether the Hortonian scaling is embedded in the spatially heterogeneous distribution of human-induced nutrient loads (total phosphorus, P tot ) reaching streams, and the biogeochemical processes and responses coupled with hydrological conditions. Our analyses on the Weser River and its six sub-catchments exhibited the existence of the Hortonian scaling patterns in the P tot input loads to receiving streams, and further their distinctions by source-types and land covers. This suggests additional evidences for self-organized distributions of human settlements and the relevant managed landscapes, as previous studies demonstrate. To investigate hierarchical structures in biogeochemical processes and responses over the river network, we used the model framework for element nutrient removal at the scale of an entire river network which consisted of well-established scaling features in hydrological responses, hydraulic geometry variables, and in-stream nutrient attenuation processes. Both the fraction and magnitude of P tot loads removed in each reach showed the Hortonian scaling features, but with opposite changing trends. As stream-order increased, the former decreased while the latter increased, under a given steady-state river discharge. This demonstrated higher efficiency of in-stream nutrient uptake in smaller streams, and greater nutrient input loads toward downstream induced by land cover patterns and aggregation through river flows. The results from the numerical model simulations were in agreement with analytical expressions presented in the current and previous studies. We acknowledge that the typical values of parameters in the hydraulic geometry functions should be carefully evaluated for other study areas, especially where river flow regimes and the associated channel geometry were significantly altered by human-built hydraulic structures. Nonetheless, the scientific significance analytically and numerically found from the large river basin warrants further investigation for generalization of the presented Hortonian scaling features in coupled hydrological and biogeochemical processes over spatially heterogeneous landscapes.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found at: The German WWTPs data for 2,000 PE were obtained from a public source provided by European Environment Agency (https://www.eea.europa.eu/data-and-maps/data/waterbase-uwwtdurban-waste-water-treatmentdirective-5#tab-european-data). The population data was from the Gridded Population of the World, version 4 Adjusted Population Count.

AUTHOR CONTRIBUTIONS
SY and PSCR conceived the study and led the manuscript writing process. SY and EB developed the model and performed the analyses. All authors contributed to discussing the results, identifying the scope of investigations, revising the manuscript, and finalizing the submitted version.

FUNDING
This study was funded by Helmholtz Centre for Environmental Research-UFZ, and by Purdue University through the Lee A. Rieth Endowment in Lyles School of Civil Engineering and the NSF-RIPS project: Resilience simulation for water, power and road networks (Award no. 1441188).