Rationalizing the Differences Among Hydraulic Relationships Using a Process‐Based Model

The use of power law forms to describe hydraulic geometry is a classic subject with a history of over 70 years. Two distinct forms of power laws have been proposed: at‐a‐station hydraulic geometry (AHG) and downstream hydraulic geometry (DHG). Although the utility of these semiempirical expressions is widely recognized, they remain poorly understood in terms of the mechanisms underlying the differences between AHG and DHG, as well as the variability among different systems. In this study, we attempt to address these basic issues. Two hypotheses are proposed: (a) the different geomorphic relationships represented by AHG and DHG result from the control of lateral adjustment of the bank and flow turbulence over short and long timescales, respectively; and (b) the systematic variability of the AHG and DHG exponents is related to the description of the frictional resistance. These two hypotheses are embedded in our theoretical models and lead to explicit functional forms for AHG and DHG. The verification of our hypotheses is based on a large data set consisting of over 550 b‐f‐m exponents and 120 power law hydraulic relations. The analysis highlights the role of uncertainties in data acquisition and theoretical/statistical explanations. In addition, the theoretical expressions of AHG also provide an explanation of at‐many‐stations hydraulic geometry (AMHG) in a physical sense. Overall, our work provides new insights into the fundamental theory of power laws and hydraulic geometry.

. Park (1977) and Rhodes (1977Rhodes ( , 1987 collated existing observations to develop a b-f-m ternary diagram consisting of over 400 AHG exponents and 150 DHG exponents (Figure 1). The AHG and DHG exponents occupy different parts of the ternary diagram with limited overlap. This implies that AHG and DHG are intrinsically representative of different conditions.
Both the AHG and DHG b-f-m exponents exhibit highly scattered distributions in the ternary diagram ( Figure 1) and no clear trend can be observed. Studies usually averaged the fitted b-f-m exponents from different systems (Ferguson, 1986;Leopold & Maddock, 1953), or introduced additional parameters (e.g., sediment grain size) attempting to construct unified power law relationships (Millar, 2005;Parker et al., 2007). Some additional insights were provided by Jung et al. (2013) and Kim and Paik (2018), who explored the variation in suspended sediment concentration for systems characterized by b d > f d and b d < f d , respectively. Efforts have also been made to deduce the analytical solutions of b-f-m exponents through hydraulic equations (Gleason, 2015). Studies on AHG generally concentrate on the variability of the hydraulic variables with the short-term rise and fall of water level (Dingman, 2007;Ferguson, 1986). For example, Ferguson (1986) preset a parabolic channel cross section and deduced that (b s , f s , m s ) = (0.23, 0.46, 0.31), with b s = f s /2. On the other hand, studies on DHG primarily emphasize the long-term morphodynamic equilibrium of the channel profiles and the associated bankfull geometry and kinematics (Cao & Knight, 1996;Finnegan et al., 2005;Huang et al., 2002;Langbein, 1963;Millar, 2005;Savenije, 2003).
In contrast to the extensive attention for b-f-m exponents, researchers have shown much less interests in the a-c-k coefficients of AHG or DHG (Equation 1). This might be attributed to the large variability in the fitted values of the a-c-k coefficients, ranging from 10 0 to 10 3 (Gleason & Smith, 2014;Rhoads, 1991;Stewardson, 2005). A notable progress was provided by Gleason and Smith (2014), who plotted a series of fitted AHG parameters along a river in a semilog coefficient-exponent diagram and observed that b s , f s , and m s indicate strong linear decrease with log s a , log s c , and log s k , respectively. Gleason and Smith (2014) termed this finding as at-many-stations hydraulic geometry (denoted as "AMHG" hereafter) and further developed a statistical tool to evaluate the discharge based on channel width data, which can be derived from remote sensing images. Barber and Gleason (2018) confirmed that AMHG is a common phenomenon among 191 natural fluvial systems in the continental United States.
Much of our understanding of hydraulic geometry has so far stemmed from fitting power laws to field measurements. However, such data-driven analyses do not address fundamental questions like: "Why do power functions represent these geomorphic relationships?," "Do b-f-m exponents vary systematically between different systems?," "Why do AHG and DHG exponents exhibit different distributions in the b-f-m ternary diagram?" and "How are the a-c-k coefficients physically related to b-f-m exponents?" (Gleason & Smith, 2014;Park, 1977;Parker et al., 2007;Rhodes, 1977Rhodes, , 1987Richards, 1973). Existing theoretical approaches often preassume the existence of power functions of the discharge (Cao & Knight, 1996;Ferguson, 1986;Finnegan et al., 2005;Huang et al., 2002;Langbein, 1963;Millar, 2005;Savenije, 2003) and/or constrain the channel cross section to be a particular shape, such as a rectangle (e.g., Cao & Knight, 1996;Huang et al., 2002;Langbein, 1963), trapezium (e.g., Eaton et al., 2004) or parabola (e.g., Dingman, 2007). Therefore, these are still room to develop a general analytical approach to address these underlying questions.
In this study, we attempt to gain some explanatory insights into the use of the power law in hydraulic geometry. Two hypotheses are proposed in terms of the hydrodynamic and sedimentary processes that may be responsible for the power law relationships exhibited by AHG and DHG (see Section 2 for a detailed description). The theoretical models are developed attempting to derive the functional forms of AHG and XU ET AL. 10.1029/2020WR029430 2 of 21 Figure 1. b-f-m ternary diagram of the at-a-station hydraulic geometry (blue dots) and downstream hydraulic geometry (red dots) exponents, redrawn from Park (1977) and Rhodes (1977Rhodes ( , 1987. DHG. Moreover, in addition to the data sets of b-f-m exponents shown in Figure 1, we collate up to 120 sets of complete power law relationships with both b-f-m exponents and a-c-k coefficients (see Tables 1 and 2). The derived functional forms are compared with the data sets, as a verification of these hypotheses and with the previous, related work by Jung et al. (2013). Our objective is to develop a more robust explanation that brings new insights into these geomorphic/hydraulic relationships.

Hypothesis 1: Fast and Slow Processes
The evolution of a channel cross section depends on two key processes: (a) erosion/deposition of stream bed driven by the flow turbulence, and (b) lateral adjustment of the bank (e.g., bank erosion/collapse) induced by gravity. The balance between these two processes has been extensively investigated in order to develop analytical solutions of the equilibrium channel profile (e.g., Ikeda & Izumi, 1990, 1991Izumi et al., 1997;Parker, 1978aParker, , 1978b. However, it was usually neglected that the turbulent flow tends to shape the channel profile slowly (De Vries, 1975), while the width adjustment related to bank erosion/collapse can be much faster (Nagata et al., 2000;Zhao et al., 2020). For instance, Simon (1992) observed that the rate of the width adjustment can reach up to 50 m/year. This is consistent with the short and long sampling frequencies of parameters used in AHG and DHG. This motivated our hypothesis that the different geomorphic relationships represented by AHG and DHG result from the control of lateral adjustment of the bank and flow turbulence, respectively.

Hypothesis 2: Ratio m/f and Frictional Resistance
Eliminating Q in H = cQ f and U = kQ m (Equation 1), a power law relationship can be deduced between the velocity and the water depth characterized by the ratio m/f (i.e.,  / m f U H ). Rhodes (1977Rhodes ( , 1987 and Ferguson (1986) used this ratio seeking to decipher the scatter in the b-f-m ternary diagram, but only a few values were investigated (e.g., m/f = 1/2, m/f = 2/3, and m/f = 1). The large data set available enable us to notice a continuous change of b-f-m exponents with m/f ( Figure 2). The AHG data approximately distribute within the interval 0.2 < m/f < 2.2 and the DHG data mostly concentrate in the interval 0 < m/f < 1.1. Polynomial regression is employed to fit the curvilinear trends of the AHG and DHG exponents (the thin lines in Figure 2). The R-squared values ( 2 R ) of the AHG trendlines are 8%, 73%, and 78% for exponents b s , f s , and m s , and are 17%, 49%, and 86% for the DHG counterparts. The relatively high R 2 values for f and m are largely determined by their implicit correlation to m/f. However, these two data sets do show different trends under the same operation, which highlights how their systematic variability with m/f differ.
The frictional resistance can be expressed in a quadratic form, say    2 D C U , where ρ is the density of water and C D is the drag coefficient. The value of C D is determined by the bed roughness length Z 0 and the water depth H, and can be evaluated using a power law, say z H , where α is a constant coefficient. It thus can be inferred that the frictional resistance is characterized by the ratio U/H q , which coincides with the geomorphic relationship  / m f U H . This elicits our second hypothesis: the systematic variability of the AHG and DHG exponents physically depends on the frictional resistance characterized by the ratio m/f (denoted as q hereafter). This hypothesis also resonates with the suggestions of several researchers (Ferguson, 1986;Leopold et al., 1964;Richards, 1973) who had attempted to link the hydraulic geometry to the frictional resistance.

Theoretical Models
An idealized channel cross section is developed to investigate the two hypotheses proposed above. The bed shear stress over the cross section is calculated based on a general momentum and continuity equations in a Cartesian coordinate system (as e.g., Yu & Knight, 1998): where x and y are the longitudinal and the lateral coordinates (m), τ is the bed shear stress (Pa), ρ is the density of water (kg/m 3 ), g is the gravity (m 2 /s), h is the water depth (m), S is the dimensionless hydraulic slope,  yx is the turbulent shear stress (Pa), and  1 is a Boolean variables, which allows to include/exclude the turbulence term. The momentum balance described by Equation 2 is valid for wide and shallow channels, where the lateral bed slope is small. Another commonly used version has been derived based on an orthogonal curvilinear coordinate system, in which the channel cross section is divided by normals perpendicular to the bed (e.g., D'Alpaos et al., 2006;Fagherazzi & Furbish, 2001;Lundgren & Jonsson, 1964;Pizzuto, 1990). The Cartesian coordinate system employed in this study facilitates the derivations of the theoretical forms.
The bed shear stress is linked to the velocity and the water depth using a general Manning-Strickler formula (Parker et al., 2007): where u is the velocity (m/s), n is Manning's coefficient, and q is corresponding to the ratio m/f. The turbulent diffusion is evaluated using an eddy viscosity closure (Falconer, 1993): where v e is the depth-averaged eddy viscosity coefficient (m 2 /s) expressed as    * e hu (Fischer, 1973), in which is a dimensionless coefficient and is the friction velocity (m/s). The evolution of the channel profile is governed by the mass conservation law: where t is time (s), Q ed is the volumetric erosion-deposition rate (m/s), Q bl is the volumetric lateral bed-slope sediment flux (m 2 /s), and  2 is a Boolean variable, which is used to include/exclude the bed-slope term. The local erosion/deposition is linearly related to the bed shear stress (Mehta, 1986): where M is a reference erosion/deposition rate (m/s), and  e is the threshold shear stress (Pa). The lateral bed-slope sediment flux Q bl , without loss of generality, is determined by the local bed shear stress and the bed slope (Parker, 1984): where N is a reference sediment flux (m 2 /s) and  b is the threshold shear stress for the lateral movement of sediments (Pa). This formula was developed to model the downslope bedload sediment flux of sand (e.g., Ikeda, 1989), while Talmon and Wiesemann (2006) extended its application to total load transport. The morphodynamic evolution governed by Equations 2-8 can be simulated numerically (as e.g., Pizzuto, 1990;Wobus et al., 2006;Xu et al., 2019), which helps to determine the convergence of the ideal systems. A detailed description of the numerical approach is provided by Xu et al. (2019) and the code is available through the https://doi.org/10.5281/zenodo.3825909.
Note that Hypothesis one is invoked into the system through the Boolean values  1 and  2 , while the interplay between the fast and slow processes are not explicitly embodied in the simulations or derivations. The strategy is to simulate the ideal equilibrium profiles fully determined by the short-term and long-term processes, respectively, and then to compare the resultant geomorphic relationships with the AHG and DHG data sets, respectively. Hence, two theoretical models are developed with (a)   1 0 and   2 1, namely, only bed-slope sediment transport is considered for the lateral adjustment of the profile (denoted as "BS" XU ET AL. 10.1029/2020WR029430 6 of 21  (2007), Dudley (2004), Harvey (1975), Jowett (1998), Lewis (1969), Mackey et al. (1998), Singh (2014), Stall and Yang (1970), and Stewardsons (2005). AHG, at-a-station hydraulic geometry.

of 21
Authors  Table 2 Values

of b-f-m Exponents and a-c-k Coefficients of DHG Data Collated From Previous Research
hereafter); and (b)   1 1 and   2 0, which emphasizes the adjustment of the channel profile to the flow turbulence (denoted as "FT" hereafter). Note that such either or cases do not exist in natural systems. An additional simulation with   1 1 and   2 1 was also performed to demonstrate the functionality of the model. Hypothesis 2 is embedded in the Manning-Strickler friction formula (Equation 4). The value of the exponent q in the friction formula is given by the ratio m/f calculated based on the data sets of b-f-m exponents.

Numerical Simulations
The equilibrium profiles resulting from the numerical simulations of BS and FT models are shown in Figure 3. The BS profile is composed of a vertical bank and a curved channel portion (see the blue line in Figure 3a). The bed shear stress peaks at the channel axis and thus the convergence of the lateral sediment flux is compensated (see the red line in Figure 3a). The channel edge is characterized by the threshold bed shear stress for the lateral movement of sediments (i.e.,    b ). The FT model results in a smoother shape without sharp demarcation between the channel and the bank (see the blue line in Figure 3b). The whole profile is characterized by an evenly distributed bed shear stress equal to the threshold value  e (see the red line in Figure 3b). This implies a balance of the hydraulic pressure and the lateral turbulent diffusion throughout the channel profile (Xu et al., 2020). In particular, the bed resistance is fully compensated by the lateral turbulent diffusion at the channel edge, where the hydraulic pressure gSh tends to 0. The simulation with both processes included (i.e.,   1 1 and   2 1) produces an intermediate equilibrium profile (Figure 3c), similar in shape to previous analytical and numerical approaches (e.g., Ikeda & Izumi, 1990;Jung et al., 2013;Parker, 1978a;Pizzuto, 1990  Note. Allen et al. (1994); Andrén (1994); Andrews (1984); Bomhof et al. (2015); Bray (1973); Dudley (2004); Elliott and Cartier (1986); Ellis and Church (2005); Emmett (1972Emmett ( , 1975; Fahnestock (1963); Griffiths (1980); Hey and Thorne (1986); King (2004); Leopold and Maddock (1953)

AHG Parameters (BS Model)
Substituting Equations 2, 7, and 8 into Equation 6   The thin lines indicate the trends of the scattered data fitted using polynomial regression. Notice that the AHG data approximately distribute within the interval, 0.2 < q < 2.5, while the DHG data concentrate in a smaller interval, 0 < q < 1.1. Hence, panel a is wider than panel b for clarity. The DHG has been redrawn from Xu et al. (2020).
and the solution gives a parabolic profile (Figure 3a): ,

DHG Parameters (FT Model)
Since   2 0 in Equation 6 for the FT model, the equilibrium state is characterized by the threshold shear stress τ e over the cross section (see the red line in Figure 3b). Substituting the equilibrium condition τ = τ e in Equation 2 with   1 1, the governing equation for the equilibrium profile achieved in the FT model is obtained: where   this system converges to a critical solution characterized by the largest width-to-depth ratio, due to eddy diffusion of momentum. Then an implicit relationship between the bed slope and the water depth can be derived following Xu et al. (2019): The functional forms of the width B, the cross-sectional area A and the discharge Q can further be derived by integration, as used for the BS model (Section 3.2): where ( Similarly, eliminating h 0 in Equation 22, the geomorphic scales predicted by the FT model also inherently satisfy power law relationships, but with different functional forms of b-f-m exponents and a-c-k coefficients: 24) 4. Discussion

b-f-m Exponents and a-c-K Coefficients
As shown in Section 3, the BS and FT models lead to explicit functional forms of b-f-m exponents and a-c-k coefficients for AHG (Equations 15 and 16) and DHG (Equations 23 and 24). The next step is to compare the theoretical expressions to parameters extracted from field measurements, including more than 550 sets of b-f-m exponents (Park, 1977;Rhodes, 1977Rhodes, , 1987 and up to 120 sets of complete power law relationships giving both b-f-m exponents and a-c-k coefficients (see Tables 1 and 2). The behavior of the hybrid model (i.e.,   1 1 and   2 1) depends on the relative influence of the turbulence and the bed-slope term in the FT and BS models respectively. This aspect is illustrated in the supporting information. The comparison is conducted in a semilog coefficient-exponent diagram following Gleason and Smith (2014)'s approach (Figure 4). The complete power law data sets are plotted in a-b, c-f and k-m tuples as scattered points, and the predicted ranges are presented as boxes. The boxes enclose most of the corresponding data points, which implies that the theoretical models predict close-to-reality width-to-depth ratio. The shapes of the boxes are different between BS and FT, while their predictions in the ranges of the a-c-k coefficients are similar to each other, both indicating that c ≈ k ≪ a. The points outside the boxes are mainly attributed to their b-f-m exponents being higher or lower than the predicted ranges, as constrained by Equations 15 and 23. For example, the BS model predicts that   0.125 0.333 s b with 0 < q < 2.5, while some of the data points indicate a value larger than 0.5 or smaller than 0.1 (blue circles in Figure 4a). Note that a larger ranges parameter would result in larger boxes.
XU ET AL.

Dynamics and Timescales
The predictions of the BS and FT models both exhibit striking consistencies with the large data sets in terms of the systematic variabilities of b-f-m exponents as well as the ranges of the a-c-k coefficients (Figures 2  and 4). Moreover, a notable improvement is that all these geomorphic relationships are derived based on empirical-theoretical formulations of the underlying processes (e.g., bed resistance, eddy turbulence), rather than a priori determining B, H or U as power functions of Q, or restricting the channel profile to a specific shape. These advantages enable us to provide an explanation on fundamental issues related to the hydraulic geometry.
The ratio m/f expresses how strongly the mean velocity varies with the water depth (i.e.,  / m f U H ) and thereby is characterized by the frictional resistance. Hence, the control of m/f on b-f-m exponents indicated in both the theoretical models and the data sets (Figure 2) highlights the role of frictional resistance in the systematic variability of both AHG and DHG. This is not surprising, since friction participates in almost all processes that governs the evolution of the channel profile. For instance, in addition to the bottom resistance (Equation 4), friction also affects the lateral momentum exchange resulting from eddy viscosity (   * e hu in Equation 5), erosion/deposition (τ, in Equation 7), and bed-slope sediment flux (τ, in Equation 8).
The different relationships exhibited by AHG and DHG also provide an example that reflects the duality of morphodynamic systems as change can be driven by different processes operating over different, short XU ET AL.   Tables 1 and 2. The ranges of the parameters are given as: 0 < q < 2.5 for the BS model and 0 < q < 1.1 for the FT model; Manning's coefficient, 0.004 s/m 1/3 < n < 0.012 s/m 1/3 ; dimensionless eddy viscosity, 0.1 < Λ < 5; threshold bed shear stresses, 0.2 Pa < τ e < 0.5 Pa and 0.1 Pa < τ b < 0.4 Pa; and the ratio of reference erosion/deposition rate and bed-slope sediment flux, 1 < N/M < 50. Note that larger ranges of the model parameters can result in bigger boxes. and long, timescales. In this context, Baar et al. (2019) suggested that different calibrations of bed-slope sediment flux are needed to model systems over different spatiotemporal scales. Our theoretical approach suggests that AHG and DHG are governed by lateral bed-slope sediment transport (BS) and FT, respectively (Figures 2 and 4). These two processes can reach an absolute balance with steady boundary conditions as proved in many theoretical and numerical studies (e.g., Parker, 1978b;Pizzuto, 1990, see also Figure 3c). However, natural variability in the discharge tends to continually perturb the system from such an equilibrium configuration. At the same time, BS and FT induce a response to the changing conditions at different speeds and their preferred hydraulic relation can be manifested over short and long timescales.

Uncertainties in Power Law Relationships
The data sets of AHG or DHG do not exactly conform to the predictions of the theoretical models. The data points of b-f-m exponents still show considerable spread around the theoretical lines ( Figure 2). Also, the predicted ranges of the a-c-k coefficients do not enclose all the data points (Figure 4). These inconsistencies may be related to the simplifications in the theoretical models, including neglecting the spatiotemporal inhomogeneity of bed materials, separating the intertwined processes (i.e., BS and FT), and the assumptions of stationary equilibria (Figure 3), which is rarely achieved in natural systems. On the other hand, the data sets themselves may also contain uncertainties. The relatively large differences between the fitted trends and the theoretical lines of AHG exponents, within the interval 0 < m s /f s < 0.2, may be attributed to the paucity of data points ( Figure 2b). As pointed out by Richard (1997), the AHG data are usually measured from relatively stable channel cross sections, where vertical changes in water level are fast compared to changes in width. While this type of channel cross sections improves the accuracy of the flow measurements, it might also lead to an underestimation of the exponent b s for AHG. Moreover, we are unable to confirm that the collated field data are measured and processed to a uniform standard, or whether the definitions of the width and the depth are consistent, as well as the sampling frequency for the time-varying parameters. These uncertainties should be considered when interpreting the results.
Moving beyond the parameters and into the underlying controversy of the existence of power law relationships. The ability of power law relationships to capture the observations for so many systems around the world has generated many attempts to invoke various forms of "theoretical" explanation. But the data scatter has led to skepticism as to the validity of such relationships and the adequacy of any such theoretical explanation has also been questioned (Bates, 1990;Richards, 1973;Teleki, 1972). In fact, as long as the descriptions of the underlying processes are semiempirical, it is impossible to theoretically demonstrate the power law relationships in a strict sense. The present derivation is no exception. For example, if a logarithmic expression of friction is used (e.g., Colebrook-White), it can be immediately deduced that the geomorphic relationships are not power laws (Julien & Wargadalam, 1995). In fact, the present derivation demonstrates that power laws provide a suitable approximation of AHG and DHG using common expressions of the underlying processes (e.g., frictional resistance and eddy viscosity closure). Notice that we do not a priori restrict the forms of B, H or U as power functions of Q.

Revisiting AMHG: Can It Be Physically Explained by the BS Model?
Gleason and Smith (2014) plotted a series of fitted AHG parameters along a river using a semilog coefficient-exponent diagram and observed that s b , s f , and s m generally indicate a strong linear decrease with log s a , log s c , and log s k , respectively. They termed this result as AMHG. Gleason and Wang (2015) mathematically explained that such linear trends result from imposing each AHG as a power law. Specifically, suppose we have K sets of geomorphic relationships,  b i i B a Q with i = 1,2,3, … , K. Taking logarithm on both sides gives: which represents K straight lines in a log-log Q-B plot (Figure 5a). Gleason and Wang (2015) proposed that such a family of AHG yield a perfect AMHG, if all the AHG curves exactly intersect the same point, say XU ET AL.

10.1029/2020WR029430
14 of 21  0 B B and  0 Q Q (e.g., red cross in Figure 5a). In this case, all the a i -b i tuples lie on a decreasing line in the semilog coefficient-exponent diagram given by (Figure 5b): Gleason and Wang (2015) then suggested that a river is characterized by a strong AMHG, if the intersections of the AHG curves are concentrated around a "centroid." Notice that this explanation essentially is based on the superposition of power law fits to the data and does not exclude cases where the Q-B data points are poorly represented by a power law. To illustrate this point, we randomly generate 20,000 sets of Q-B sequences within an interval centered on   in Figure 5a). However, the corresponding a-b tuples (gray circles in Figure 5b) still show a strong decrease trends as the "perfect" AMHG suggests (compare the black dashed line and the red line in Figure 5b).
Notice that a large data set (e.g., 500 Q-B sequences here) is needed to mimic the AMHG trend with random data and the fitted a-b tuples exhibit a divergent distribution (see gray circles in Figure 5b). On the contrary, the measured Q-B sequences can result in more concentrated patterns and even a small data set can exhibit a strong AMHG trend (e.g., less than 20 Q-B sequences as reported by Gleason & Smith, 2014). Gleason and Wang (2015) supposed that this might be attributed to the hydraulic self-similarity of the cross sections along a river, which promotes the concentration of intersections of the AHG curves (i.e., Equation 26).
Recently, Brinkerhoff et al. (2019) attempted to confirm the physical existence of AMHG, based on Dingman (2007)'s theoretical model of AHG. Here our concern is more specifically on whether such "self-similarity" can be physically explained by the BS model.
The AMHG data (i.e., a series of AHG parameters along a river) of the six rivers reported by Gleason and Smith (2014) are compared with the BS model ( Figure 6). The consistency in the ranges of a-b, c-f, and k-m tuples (Figure 6a) resembles the comparison presented in Section 4.1 (Figure 4a). In terms of the trends, the BS model can produce any curve within the predicted range through calibration. Here, we configure the BS model taking advantage of Equation 15, in which variations of the b-f-m exponents only depend on q. Then each river model is calibrated for hydrodynamic and sedimentary parameters n, M, N,  e and  b , which are assumed to be constant along the river. The resultant a-b, c-f, and k-m curves approximately show linearly decreasing trends, similar to the distribution of the data points ( Figure 6). This demonstrates that the BS model, as configured, is able to reproduce the trends exhibited by AMHG. Notice that here we do not pursue a best fit to the observed data points. The parameters of the model are picked using a simple routine that randomly adjusts parameters (within the ranges derived in Section 4.1) until a minimum is found. In fact, these parameters may vary spatially in natural systems. For example, bed materials are known to vary in a spatially coherent way (Blom et al., 2016). A better elaboration of such coherent variability may further improve the explanatory power of the proposed model.

Limitations and Suggestions for Future Research
As highlighted in Sections 4.3 and 4.4, uncertainties are implicit in many aspects in the research of hydraulic geometry, including the inhomogeneity of natural systems, the errors in data acquisition, and the usage of semiempirical expressions for underlying processes. These limitations can be misleading when determining the parameters of the hydraulic relations and even the legitimacy of the power laws. Therefore, future research should be devoted to eliminating or reducing such uncertainties. Specifically, we consider that the following aspects are worthy of in-depth investigation: 1. The Cartesian coordinate system is employed in the description of the momentum balance (Equation 2). This approach is valid when the lateral bed slopes of the channel cross sections are generally small. However, for the almost vertical bank formed in the BS profile (e.g., Figure 3a), considering only the vertical profiles of the velocity leads to the underestimation of the bed shear stress on the vertical bank and hence the channel width. A more accurate momentum balance can be achieved using an orthogonal curvilinear coordinate system (e.g., Fagherazzi & Furbish, 2001;Lundgren & Jonsson, 1964;Pizzuto, 1990). 2. We assume a specific friction formula (Equation 4), an eddy viscosity closure assuming an identical turbulent transfer (Equation 5), and a constant threshold shear stress (Equations 7 and 8). Recent works have revealed that a river cross section with heterogeneous bed and bank material deviates significantly from such ideal threshold configuration (Dunne & Jerolmack, 2020;Francalanci et al., 2020). The present derivation in fact is a simplified abstraction useful to emphasize the overall trends of fluvial hydraulic geometry. Future models may improve the predictions with more elaborated description of such heterogeneity. 3. The derived theoretical expressions involve several parameters (q, n, Λ,  e ,  b , M, and N) calibrated using short-term laboratory experiments or field measurements (Section 4.1). New approach is needed to verify the validity of such calibrations over short and long timescales. 4. The analyses of AMHG imply that the parameters used to model the underlying physical processes (e.g., n,  e ,  b , M, and N) may exhibit a spatial coherence as supposed by Gleason and Wang (2015) (see also XU ET AL.

10.1029/2020WR029430
16 of 21 Section 4.4). This conjecture needs to be confirmed through a careful inspection of the sediment properties along a river, as well as the associated variability of these parameters. 5. Although timescale is highlighted, the specific timescales characterizing BS or FT are not given, nor is the timescale embedded in the derived theoretical expressions. Hence, the variability in the hydraulic relations caused by the interaction of processes operating over different timescales needs further attention. 6. The above suggestions urge us to increase the detail and number of variables measured when testing power law relations. Also, developing standards for the definition and measurement of variables will reduce this source of uncertainty. 7. We stress that the power law formed AHG and DHG are approximations of natural fluvial channels, and that the power law fitting does erase many details relevant to realistic channel geometries. Developing XU ET AL.

10.1029/2020WR029430
17 of 21 Figure 6. Semilog coefficient-exponent diagram of the at-many-station hydraulic geometry data extracted from Gleason and Smith (2014), and the predictions of bed-slope model. Panel (a) shows the predicted ranges (boxes) and the integrated data of all four rivers (points, as Figure 4a). Panels (b)-(e) compare the predicted trends (lines) to the four rivers (points) separately using different parameters. Note that the length of each line is limited by the predicted range of b-f-m exponents (Equation 23).
new theoretical descriptions beyond the power law paradigm is a challenging and valuable direction for future work.

Conclusions
Theoretical expressions of AHG and DHG have been derived without any a priori assumptions on the form of the relationship or the cross-sectional shape. The resulting expressions therefore have dimensionally consistent and physically meaningful coefficients. The functional forms of b-f-m exponents and a-c-k coefficients have been verified using over 550 b-f-m exponents and 120 power law hydraulic relations from rivers around the world. In addition, the predicted AHG a-b, c-f and k-m tuples agree well with the at-manystations hydraulic geometry (AMHG) data exhibiting linearly decreasing trends. The main findings of this study are summarized as follows: 1. The trends of the AHG and DHG exponents are revealed and they both depend on the frictional resistance. 2. The large variability of the a-c-k coefficients results from the number of physical and geometric variables present in the theoretical expressions. 3. The different behaviors of AHG and DHG are attributed to the bed-slope sediment transport and the flow turbulence, respectively. 4. Finally, the analysis highlights the uncertainties in the power law formed hydraulic geometry possibly resulting from data acquisition and/or theoretical assumptions.

Data Availability Statement
The code developed in this study are available through the https://doi.org/10.5281/zenodo.3825909. Data sets for this research compiled from the published papers can be downloaded from https://coastalhub. science/data.