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 km²; 8M population). Archived data for point-/diffuse-sources of total Phosphorus (P_tot) input loads were combined with numerical and analytical model simulations of coupled hydrological and biogeochemical processes for in-stream P_tot removal at the network scale. We find that P_tot 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 P_tot 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 P_tot removal, in larger streams the magnitude of P_tot 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 self-organized, 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 stream-order 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²; 8M population) served as the case study of a human-dominated 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², ~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³/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² 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).

FIGURE 1

Figure 1. Maps of human-related attributes distributed over the Weser River Basin of which location within Germany is shown in the most left figure. (A) Distribution of ~8.4 million people as the number of persons per grid with a range from 1 to 25. River network with the final 7th order is overlaid with black lines of which thickness and brightness are proportional to stream-orders. Dashed and solid triangles indicate the outlets of four sub-basins with 5th-order (F, W, L, and A for Fulda, Werra, Leine, and Aller) and two sub-basins with 6th-order (W and L for Weser and Leine), respectively. (B) Distribution of ~845 WWTPs categorized as five class-sizes of which color-codes are given in the legend. (C) Distribution of five types of land-cover of which areas are visualized with different colors [(B,C) were adapted from Yang et al. Copyright (2021), with permission from Elsevier].

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 CLC-class-1, 1,576 tons P/yr for CLC-class-2, and 485 tons P/yr for CLC-class-3 (UBA, 2010, 2017). In turn, P_tot input loads per unit-area Ψ [kg P/yr/km²] 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.

TABLE 1

Table 1. Information for the Weser River Basin and six sub-basins.

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):

$\begin{array}{ll}{L}_{in,i}=L_{PS,i}+L_{DS,i}+{\sum }_j{L}_{out,j}& (1)\end{array}$

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:

$\begin{array}{ll}{L}_{out,i}=L_{in,i}{e}^{-\frac{v_{f,i}}{d_i}{\tau }_i}& (2)\end{array}$

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 reach-scale residence time τ_i is straightforwardly estimated through the geomorphological and hydrological conditions at the reach i :

$\begin{array}{ll}{\tau }_{\text{i}}=w_i{d}_i{l}_i/Q_i& (3)\end{array}$

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:

$\begin{array}{ll}{f}_{R,i}=L_{r,i}/L_{in,i}& (4)\end{array}$

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:

$\begin{array}{ll}{Q}_{\text{i}}=q{A}_i& (5)\end{array}$

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:

$\begin{array}{ll}{w}_{\text{i}}=K_w{Q}_i^{{\alpha }_w}& (6)\end{array}$

$\begin{array}{ll}{d}_{\text{i}}=K_d{Q}_i^{{\alpha }_d}& (7)\end{array}$

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\left[{\mathrm{\text{m}}}^{1-3{\alpha }_w}{\mathrm{\text{s}}}^{{\alpha }_w}\right]$ and K_d = 0.25 [${\text{m}}^{1-3{\alpha }_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/semi-natural areas [L_FO(P)]; (2) Summed P_tot loads from both point- and 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 scale-invariance in spatial hierarchy of averaged attributes over stream-orders (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:

$\begin{array}{ll}\langle Z\rangle ~e^{h\omega }& (8)\end{array}$

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 scale-invariant 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:

$\begin{array}{llll}{R}_Z& =& \{\begin{array}{c}{\langle Z\rangle }_{\omega }/{\langle Z\rangle }_{\omega +1}(forh<0)\\ {\langle Z\rangle }_{\omega +1}/{\langle Z\rangle }_{\omega }(forh>0)\end{array}& (9)\end{array}$

$\begin{array}{ll}or{R}_Z=e^{\left|h\right|}.& (10)\end{array}$

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:

$\begin{array}{ll}{R}_{L\text{\_}TOT}={(R_{L\text{\_}PS})}^{{\delta }_{PS}}{(R_{L\text{\_}UR})}^{{\delta }_{UR}}{(R_{L\text{\_}AG})}^{{\delta }_{AG}}{(R_{L\text{\_}FO})}^{{\delta }_{FO}}& (11)\end{array}$

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.

Results

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² > 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² > 0.9, p < 0.05) followed by diffuse-sources from urban surface areas (R_{L_UR} = 1.9; R² > 0.9, p < 0.05), agricultural areas (R_{L_AG} = 1.6; R² ≈ 0.7, p < 0.05), and forested/semi-natural areas (R_{L_FO} = 1.4; R² ≈ 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.

FIGURE 2

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.

TABLE 2

Table 2. Hortonian scaling ratios for the Weser River Basin and six 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).

FIGURE 3

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.

Hortonian Scaling of Biogeochemical Processes

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 in-stream 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 stream-orders (Figure 4B), decreasing with higher ω, with R_fR = 1.9 (R² > 0.9, p < 0.05).

FIGURE 4

Figure 4. Distributions of the removed P_tot loads along the Weser River network under the mean river discharge. (A) Entire basin-scale distribution for the fraction f_R of P_tot loads removed in each reach (L_r) to input mass flux to the same reach (L_in). (B) Stream-order ω based distribution for f_R. (C) Distribution of the fraction of total P_tot loads removed by upstream catchment F_R at a given drainage area A. (D) The ω-based distribution for L_r. (E) Distribution of removed P_tot loads per unit length in each reach, ζ_r, over A. For visual efficiency, values of gray dots in (B,D) were extracted from the most downstream point of each stream for a given ω. Hortonian scaling equation and the corresponding scaling ratio are in (B,D). Blue dots in (C,E) indicated the model results for each stream reach of the Weser River network. Dashed black line in (C) and (E) was plotted by following the analytical scaling of Equations (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.

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) (i.e., $1-F_R=1-{A_T}^{1-\alpha }{A}^{\alpha -1}$, 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 on average, resulting R_Lr = 1.8 (R² > 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) (i.e., ${\xi }_R\cong (K_Y{K}_w{v}_f/q){A_T}^{1-\alpha }{A}^{\alpha +{\alpha }_w-1}$, 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, 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.

Discussion

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 (Jawitz et al., 2020). 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 6th-order and four 5th-order, of the Weser River. The clustered 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 class-sizes 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 flow-weighted 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 intra-annual 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-treatment-directive-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).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

SY and DB acknowledge the support of the Helmholtz Climate Initiative Project funded by the Helmholtz Association. PSCR acknowledges the financial support of the Lee A. Rieth Endowment in Lyles School of Civil Engineering at Purdue University. EB acknowledges the support of the Ca' Foscari University of Venice, Italy. We would also like to thank the reviewers for their constructive comments and suggestions.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frwa.2021.693056/full#supplementary-material

References

Aissa-Grouz, N., Garnier, J., and Billen, G. (2018). Long trend reduction of phosphorus wastewater loading in the Seine: determination of phosphorus speciation and sorption for modeling algal growth. Environ. Sci. Pollut. R. 25, 23515–23528. doi: 10.1007/s11356-016-7555-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Basu, N. B., Rao, P. S. C., Thompson, S. E., Loukinova, N. V., Donner, S. D., Ye, S., et al. (2011). Spatiotemporal averaging of in-stream solute removal dynamics. Water Resour. Res. 47:W00J06. doi: 10.1029/2010WR010196

CrossRef Full Text | Google Scholar

Bertuzzo, E., Helton, A. M., Hall, R. O., and Battin, T. J. (2017). Scaling of dissolved organic carbon removal in river networks. Adv. Water Resour. 110, 136–146. doi: 10.1016/j.advwatres.2017.10.009

CrossRef Full Text | Google Scholar

BfG (2019). Water Levels–Data From Selected Gauging statIons on German Federal Waterways and Shipping Administration (WSV). Koblenz: German Federal Institute of Hydrology (Bundesanstalt für Gewässerkunde). Available online at: https://www.bafg.de/EN/06_Info_Service/01_WaterLevels/waterlevels.html (accessed December 2020).

Biggs, J., von Fumetti, S., and Kelly-Quinn, M. (2017). The importance of small waterbodies for biodiversity and ecosystem services: implications for policy makers. Hydrobiologia 793, 3–39. doi: 10.1007/s10750-016-3007-0

CrossRef Full Text | Google Scholar

Botter, G., Basso, S., Rodriguez-Iturbe, I., and Rinaldo, A. (2013). Resilience of river flow regimes. Proc. Natl. Acad. Sci. U.S.A. 110, 12925–12930. doi: 10.1073/pnas.1311920110

PubMed Abstract | CrossRef Full Text | Google Scholar

Botter, G., Basu, N. B., Zanardo, S., Rao, P. S. C., and Rinaldo, A. (2010). Stochastic modeling of nutrient losses in streams: interactions of climatic, hydrologic, and biogeochemical controls. Water Resour. Res. 46:WR008758. doi: 10.1029/2009WR008758

CrossRef Full Text | Google Scholar

Büttner, O. (2020). DE-WWTP-Data Collection of Wastewater Treatment Plants of Germany (Status 2015, Metadata). HydroShare.

Dodds, P. S., and Rothman, D. H. (2000). Geometry of river networks. I. Scaling, fluctuations, and deviations. Phys. Rev. E 63:016115. doi: 10.1103/PhysRevE.63.016115

PubMed Abstract | CrossRef Full Text | Google Scholar

Doyle, M. W. (2005). Incorporating hydrologic variability into nutrient spiraling. J. Geophys. Res. Biogeosci. 110:G01003. doi: 10.1029/2005JG000015

CrossRef Full Text | Google Scholar

EEA (2017). Urban wast water treatment. Copenhagen: European Environment Agency. Available online at: https://www.eea.europa.eu/data-and-maps/indicators/urban-waste-water-treatment/urban-waste-water-treatment-assessment-4 (accessed December 2020).

EEC (1991). Council Directive of 21 May 1991 Concerning Urban Waste Water Treatment (91/271/EEC). Official Journal of the European Communities. Available online at: http://ec.europa.eu/environment/water/water-urbanwaste/index_en.html (accessed December 2020).

Google Scholar

Ensign, S. H., and Doyle, M. W. (2005). In-channel transient storage and associated nutrient retention: evidence from experimental manipulations. Limnol. Oceanogr. 50, 1740–1751. doi: 10.4319/lo.2005.50.6.1740

CrossRef Full Text | Google Scholar

Fang, Y., Ceola, S., Paik, K., McGrath, G. S., Rao, P. S. C., Montanari, A., et al. (2018). Globally universal fractal pattern of human settlements in river networks. Earth's Future 6, 1134–1145. doi: 10.1029/2017EF000746

CrossRef Full Text | Google Scholar

Ferrazzi, M., and Botter, G. (2019). Contrasting signatures of distinct human water uses in regulated flow regimes. Environ. Res. Commun. 1:071003. doi: 10.1088/2515-7620/ab3324

CrossRef Full Text | Google Scholar

Gall, H. E., Park, J., Harman, C. J., Jawitz, J. W., and Rao, P. S. C. (2013). Landscape filtering of hydrologic and biogeochemical responses in managed catchments. Landsc. Ecol. 28, 651–664. doi: 10.1007/s10980-012-9829-x

CrossRef Full Text | Google Scholar

Gupta, V. K., Ayalew, T. B., Mantilla, R., and Krajewski, W. F. (2015). Classical and generalized Horton laws for peak flows in rainfall-runoff events. Chaos Interdiscip. J. Nonlinear Sci. 25:075408. doi: 10.1063/1.4922177

PubMed Abstract | CrossRef Full Text | Google Scholar

Gupta, V. K., and Mesa, O. J. (2014). Horton laws for hydraulic–geometric variables and their scaling exponents in self-similar Tokunaga river networks. Nonlinear Process. Geophys. 21, 1007–1025. doi: 10.5194/npg-21-1007-2014

CrossRef Full Text | Google Scholar

Hack, J. T. (1957). Studies of longitudinal stream profiles in Virginia and Maryland. United States Geol. Surv. Prof. Pap. 294-B, 45–97. doi: 10.3133/pp294B

CrossRef Full Text | Google Scholar

Helton, A. M., Hall Jr, R. O., and Bertuzzo, E. (2018). How network structure can affect nitrogen removal by streams. Freshw. Biol., 63, 128–140. doi: 10.1111/fwb.12990

CrossRef Full Text | Google Scholar

Hirt, U., Kreins, P., Kuhn, U., Mahnkopf, J., Venohr, M., and Wendland, F. (2012). Management options to reduce future nitrogen emissions into rivers: a case study of the Weser river basin, Germany. Agric. Water Manag. 115, 118–131. doi: 10.1016/j.agwat.2012.08.005

CrossRef Full Text | Google Scholar

Horton, R. E. (1945). Erosional development of streams and their drainage basins; hydrophysical approach to quantitative morphology. Geol. Soc. Am. Bull. 56, 275–370. doi: 10.1130/0016-7606(1945)56[275:EDOSAT]2.0.CO;2

CrossRef Full Text | Google Scholar

IWRM-net (2010). Innovative Instruments and Institutions in Implementing the Water Framework Directive: The German Case Study. Research Report of Integrated Water Resource Management (IWRM), No. 2.2. IWRM.

Jawitz, J. W., Desormeaux, A. M., Annable, M. D., Borchardt, D., and Dobberfuhl, D. (2020). Disaggregating landscape-scale nitrogen attenuation along hydrological flow paths. J. Geophys. Res. Biogeosci. 125:e2019JG005229. doi: 10.1029/2019JG005229

CrossRef Full Text | Google Scholar

Jensen, C. K., McGuire, K. J., McLaughlin, D. L., and Scott, D. T. (2019). Quantifying spatiotemporal variation in headwater stream length using flow intermittency sensors. Environ. Monit. Assess. 191:226. doi: 10.1007/s10661-019-7373-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, D. H., Rao, P. S. C., Kim, D., and Park, J. (2016). 1/f noise analyses of urbanization effects on streamflow characteristics. Hydrol. Process. 30, 1651–1664. doi: 10.1002/hyp.10727

CrossRef Full Text | Google Scholar

Klein, M. (1981). Drainage area and the variation of channel geometry downstream. Earth Surf. Process. Landf. 6, 589–593. doi: 10.1002/esp.3290060608

CrossRef Full Text | Google Scholar

Knighton, D. (1998). Fluvial Forms and Processes. London: Edward Arnold.

Google Scholar

Kumar, R., Heße, F., Rao, P. S. C., Musolff, A., Jawitz, J. W., Sarrazin, F., et al. (2020). Strong hydroclimatic controls on vulnerability to subsurface nitrate contamination across Europe. Nat. Commun. 11:6302. doi: 10.1038/s41467-020-19955-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Leopold, L. B., and Maddock, T. J. (1953). The Hydraulic Geometry of Stream Channels and Some Physiographic Implications. United States Geological Survey Professional Paper 252. Washington, DC: U.S. Government Printing Office.

Google Scholar

Liu, X., Beusen, A. H. W., Van Beek, L. P. H., Mogollón, J. M., Ran, X., and Bouwman, A. F. (2018). Exploring spatiotemporal changes of the Yangtze River (Changjiang) nitrogen and phosphorus sources, retention and export to the East China Sea and Yellow Sea. Water Res. 142, 246–255. doi: 10.1016/j.watres.2018.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Mantilla, R., Gupta, V. K., and Mesa, O. J. (2006). Role of coupled flow dynamics and real network structures on Hortonian scaling of peak flows. J. Hydrol. 322, 155–167. doi: 10.1016/j.jhydrol.2005.03.022

CrossRef Full Text | Google Scholar

Miyamoto, H., Hashimoto, T., and Michioku, K. (2011). Basin-wide distribution of land use and human population: stream order modeling and river basin classification in Japan. Environ. Manage. 47, 885–898. doi: 10.1007/s00267-011-9653-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Mockler, E. M., Deakin, J., Archbold, M., Gill, L., Daly, D., and Bruen, M. (2017). Sources of nitrogen and phosphorus emissions to Irish rivers and coastal waters: estimates from a nutrient load apportionment framework. Sci. Total Environ. 601–602, 326–339. doi: 10.1016/j.scitotenv.2017.05.186

PubMed Abstract | CrossRef Full Text | Google Scholar

Paik, K., and Kumar, P. (2004). Hydraulic geometry and the nonlinearity of the network instantaneous response. Water Resour. Res. 40:WR002821. doi: 10.1029/2003WR002821

CrossRef Full Text | Google Scholar

Paik, K., and Kumar, P. (2008). Emergence of self-similar tree network organization. Complexity 13, 30–37. doi: 10.1002/cplx.20214

CrossRef Full Text | Google Scholar

Peckham, S. D. (1995). New results for self-similar trees with applications to river networks. Water Resour. Res. 31, 1023–1029. doi: 10.1029/94WR03155

CrossRef Full Text | Google Scholar

Ran, L., Wang, S., and Lu, X. X. (2012). Hydraulic geometry change of a large river: a case study of the upper Yellow River. Environ. Earth Sci. 66, 1247–1257. doi: 10.1007/s12665-011-1336-x

CrossRef Full Text | Google Scholar

Rice, J., and Westerhoff, P. (2017). High levels of endocrine pollutants in US streams during low flow due to insufficient wastewater dilution. Nat. Geosci. 10:587. doi: 10.1038/ngeo2984

CrossRef Full Text | Google Scholar

Richardson, J. S. (2019). Biological diversity in headwater streams. Water 11:366. doi: 10.3390/w11020366

CrossRef Full Text | Google Scholar

Rodríguez-Iturbe, I., Ijjász-Vásquez, E. J., Bras, R. L., and Tarboton, D. G. (1992). Power law distributions of discharge mass and energy in river basins. Water Resour. Res. 28, 1089–1093. doi: 10.1029/91WR03033

CrossRef Full Text | Google Scholar

Rodríguez-Iturbe, I., and Rinaldo, A. (2001). Fractal River Basins: Chance and Self-Organization. Cambridge: Cambridge University Press.

Google Scholar

Schumm, S. A. (1956). Evolution of drainage systems and slopes in badlands at Perth Amboy, New Jersey. Geol. Soc. Am. Bull. 67, 597–646 doi: 10.1130/0016-7606(1956)67[597:EODSAS]2.0.CO;2

CrossRef Full Text | Google Scholar

Shin, Y. H., and Julien, P. Y. (2010). Changes in hydraulic geometry of the Hwang River below the Hapcheon Re-regulation Dam, South Korea. Int. J. River Basin Manag. 8, 139–150. doi: 10.1080/15715121003651252

CrossRef Full Text | Google Scholar

Speckhann, G. A., Kreibich, H., and Merz, B. (2021). Inventory of dams in Germany. Earth Syst. Sci. Data 13, 731–740. doi: 10.5194/essd-13-731-2021

CrossRef Full Text | Google Scholar

Strahler, A. N. (1957). Quantitative analysis of watershed geomorphology. EOS Trans. AGU 38, 913–920. doi: 10.1029/TR038i006p00913

CrossRef Full Text | Google Scholar

Thober, S., Kumar, R., Wanders, N., Marx, A., Pan, M., Rakovec, O., et al. (2018). Multi-model ensemble projections of European river floods and high flows at 1.5, 2, and 3 degrees global warming. Environ. Res. Lett. 13:014003. doi: 10.1088/1748-9326/aa9e35

CrossRef Full Text | Google Scholar

UBA (2010). Calculation of Emissions into Rivers in Germany Using the MONERIS Model-Nutrients, Heavy Metals, and Polycyclic Aromatic Hydrocarbons. The German Federal Environmental Agency (UBA). Available online at: https://www.umweltbundesamt.de/en/publikationen/calculation-of-emissions-into-rivers-in-germany (accessed December 2020).

UBA (2017). Waters in Germany: Status and Assessment. The German Federal Environmental Agency (UBA). Available online at: https://www.umweltbundesamt.de/en/publikationen/waters-in-germany (accessed December 2020).

Williams, G. P., and Wolman, M. G. (1984). Downstream Effects of Dams on Alluvial Rivers. Geological Survey Professional Paper. Washington, DC: U.S. Government Printing Office.

Google Scholar

Wohl, E. (2017). The significance of small streams. Front. Earth Sci. 11, 447–456. doi: 10.1007/s11707-017-0647-y

CrossRef Full Text | Google Scholar

Wollheim, W. M., Stewart, R. J., Aiken, G. R., Butler, K. D., Morse, N. B., and Salisbury, J. (2015). Removal of terrestrial DOC in aquatic ecosystems of a temperate river network. Geophys. Res. Lett. 42, 6671–6679. doi: 10.1002/2015GL064647

CrossRef Full Text | Google Scholar

Wollheim, W. M., Vörösmarty, C. J., Peterson, B. J., Seitzinger, S. P., and Hopkinson, C. S. (2006). Relationship between river size and nutrient removal. Geophys. Res. Lett. 33:L06410. doi: 10.1029/2006GL025845

CrossRef Full Text | Google Scholar

Yang, S., Büttner, O., Jawitz, W. J., Kumar, R., Rao, P. S. C., and Borchardt, D. (2019a). Spatial organization of human population and wastewater treatment plants in urbanized river basins. Water Resour. Res. 55, 6138–6152. doi: 10.1029/2018WR024614

CrossRef Full Text | Google Scholar

Yang, S., Büttner, O., Kumar, R., Jäger, C. G., Jawitz, J. W., Rao, P. S. C., et al. (2019b). Spatial patterns of water quality impairments from point source nutrient loads in Germany's largest national River Basin (Weser River). Sci. Total Environ. 697:134145. doi: 10.1016/j.scitotenv.2019.134145

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, S., and Paik, K. (2017). New findings on river network organization: law of eigenarea and relationships among hortonian scaling ratios. Fractals 25:1750029. doi: 10.1142/S0218348X17500293

CrossRef Full Text | Google Scholar

Zink, M., Kumar, R., Cuntz, M., and Samaniego, L. (2017). A high-resolution dataset of water fluxes and states for Germany accounting for parametric uncertainty. Hydrol. Earth Syst. Sci. 21, 1769–1790. doi: 10.5194/hess-21-1769-2017

CrossRef Full Text | Google Scholar

Appendix A

Analytical Derivation for the Hortonian Scaling in the Fraction of P_tot Loads Removed in Each Stream Reach (f_R)

In this appendix, we show an analytical derivation of the Hortonian scaling in f_R,.

Substituting Equation (2) into Equation (3) yields :

$\begin{array}{l}{f}_R=1-e^{(-v_{f}l{K}_W{Q}^{{\alpha }_W-1})}\\ \cong -v_{f}l{K}_W{Q}^{{\alpha }_W-1}=M{Q}^{{\alpha }_W-1}=g(Q),& (\mathrm{A1})\end{array}$

where the approximation originates from the Taylor series of the exponential function up to the first-order term (Ensign and Doyle, 2005). The notations of reach i in f_R and Q are omitted for the generalized function of g(Q) with the constant M = –v_f lK_W, by assuming the difference in l from diagonal or orthogonal flow direction is negligible. The power-law exponent of (α_w – 1) ranges between-1 and 0. Given the well-known Hortonian scaling for the mean value of Q, <Q>, with its scaling exponent h_Q, that is:

$\begin{array}{ll}\langle Q\rangle ~e^{h_Q\omega },& (\mathrm{A2})\end{array}$

we can express the mean of f_R, <g(Q)>, as the Taylor series of the g(Q) by taking until the first term at a real constant μ:

$\begin{array}{ll}\langle g(Q)\rangle \cong g(\mu )+g^{\prime }(\mu )(\langle Q\rangle -\mu ).& (\mathrm{A3})\end{array}$

When μ = <Q>, <g(Q)> scales with <Q> as:

$\begin{array}{ll}\langle g(Q)\rangle \cong g(\langle Q\rangle )=M{\langle Q\rangle }^{{\alpha }_W-1}.& (\mathrm{A4})\end{array}$

Substituting Equation (A2) into Equation (A4) yields an exponential scaling of <f_R> over ω as:

$\begin{array}{ll}\langle f_R\rangle =\langle g(Q)\rangle ~e^{({\alpha }_W-1)h_Q\omega }.& (\mathrm{A5})\end{array}$

Therefore, 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 and hydraulic geometry attributes for stream width are given.