the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Power law between the apparent drainage density and the pruning area
Soohyun Yang
Kwanghun Choi
Download
 Final revised paper (published on 18 Jul 2024)
 Supplement to the final revised paper
 Preprint (discussion started on 23 Nov 2023)
 Supplement to the preprint
Interactive discussion
Status: closed

RC1: 'Comment on hess2023271', Anonymous Referee #1, 10 Jan 2024
General:
The authors analyzed the relationship between the pruning area and the drainage density in a variety of delineated river networks in the US. The theoretical derivations as well as the analysis of the delineated networks suggest that a power law represents this relationship.
This finding is very interesting as it provides a new insight into the organization of river networks. Conceptually the law expresses bifurcation and reach length ratios of the Horton laws, but its application is simpler as it requires no ordering of streams. The background of scaling laws in drainage networks is well elaborated and the theoretical derivations are also plausible. The whole paper has a strong focus on numbers (the exponents of the power laws) and is therefore a bit difficult to read fluently if the concepts are not well known to the reader and I have therefore two minor comments:
You mention that the theoretical values of epsilon and h should be 0.5 but they seem to be below 0.5 (epsilon) and above 0.5 (h). You then mention that the theoretical values are related to dimensional consistency and that the found dimensional inconsistency is related to fractal dimension of topography. Although I understand your reasoning it is a bit difficult to follow. What is the fractal dimension of topography? I think it could help to try to underly these concepts additionally with physical meaning, e.g. a stream is sinusoidal because total loss of energy is minimized, similarly one could argue that a network tends to be compacted in hierarchy because.... It is noteworthy to include a sentence on how fractals and flow connectivity are interrelated and what this physically means (energy).
Another minor comment I have regarding the range of application of the found power law. There seem to be different minimum pruning areas for the law to be applicable. How can this be justified, and could there be a companion power law for smaller pruning areas? If so or if not, why?
Specific:
L22: Is there some physical explanation for this?
L26: Isn´t this the question? How can the topography be “given” if these apparent power laws can be found? There seems to be some organization, so I would suggest to mention this here.
L47: some concepts might need more explanation, what is valley transmissivity in this context?
L85: I think some simpler variable for cumulative length might be better for readability.
L106: Some additional explanation for this equation would help the reader.
L108: What is meant by such a tradeoff?
L171: Is the analysis sensitive to the method used to derive flow direction?
L175: Reach or stream?
L188: You mention that you applied a linear fit. This refers to the logarithms? As one can also see in Fig. 2a this fit is only valid for a certain interval of values, therefore how did you derive the lower cutoff value?
L225: Fig. 2a: Different A0 values have been found, what might the range between the minimum A0 and the maximum A0 be related to?
Citation: https://doi.org/10.5194/hess2023271RC1 
AC1: 'Reply on RC1', Kyungrock Paik, 16 Feb 2024
We thank Referee 1 for constructive comments. Our response to each comment is listed below. We will reflect any change addressed below into the revised manuscript which is requested in the subsequent step.
Comment 1: You mention that the theoretical values of epsilon and h should be 0.5 but they seem to be below 0.5 (epsilon) and above 0.5 (h). You then mention that the theoretical values are related to dimensional consistency and that the found dimensional inconsistency is related to fractal dimension of topography. Although I understand your reasoning it is a bit difficult to follow. What is the fractal dimension of topography? I think it could help to try to underly these concepts additionally with physical meaning, e.g. a stream is sinusoidal because total loss of energy is minimized, similarly one could argue that a network tends to be compacted in hierarchy because.... It is noteworthy to include a sentence on how fractals and flow connectivity are interrelated and what this physically means (energy).
Response 1: Fractal dimensions of river networks vary between 1 and 2, as explained in L240241 and L260 of the original manuscript. The positive noninteger dimension indicates that, in a planar view, a single stream in a given river network is more winding than a straight line, which has a dimension of 1, and simultaneously the entire layout of the network does not completely fill a given surface, which has a dimension of 2.
Regarding the comment to explain the fractal features of river networks based on their physical means, we fully agree with your idea and the associated reason. Particularly, we wish if we could clearly elucidate “which physical mechanism results in which fractal feature”, like the example statement you provided. However, it must be acknowledged that many hypotheses have been proposed to unravel which physical processes induce river network structures to be fractal, resulting in persistent debates on the topic (refer to L302306 of the original manuscript). Thereby, it is hardly to explain the causality based on a onetoone approach. Nonetheless, to do our best within the scope of this study, we clarified the underlying mechanical principle out of proposed hypotheses so far (refer to L306308 of the original manuscript).
Overall, to incorporate your constructive comment in the revised manuscript, we will improve the original L62 as follows. It aims to provide additional introductory explanation on the fractal dimension in the context of this research:
“This has brought the introduction of the fractal dimension (Mandelbrot, 1977), whose values for river networks range between 1 and 2 (e.g., Feder, 1988). Note that further detailed explanations are provided in Sect. 4.”
Comment 2: Another minor comment I have regarding the range of application of the found power law. There seem to be different minimum pruning areas for the law to be applicable. How can this be justified, and could there be a companion power law for smaller pruning areas? If so or if not, why?
Response 2: You correctly understood that individual catchments analyzed in this study have distinct values of minimum pruning area which is eventually equivalent to a source area A_{o}. For a given catchment, A_{o }is a constant determined as the median of channel forming areas extracted from the National Hydrography Dataset Plus Version 2 (NHDPlusV2). The applied A_{o} values in this study can be justified because NHDPlusV2 represents the bluelines of river networks spanning the contiguous US (refer to L174177 of the original manuscript).
For our analysis, we selected the median of all channel forming areas by following this logic: the median represents 50% frequency of total cases of channel forming areas (i.e., half of total case). As a typical starting point, we wanted to comply with the criterion for at least half of all cases to define river networks. The threshold 50% can therefore be directly justified as a neutral choice to set spatially constant source area for a given catchment.
Specific comments:
Comment: L22: Is there some physical explanation for this?
Response: Yes. The precipitation effectiveness index (PE index) is expressed as the sum of the ratio of total monthly precipitation to total monthly evaporation over a year (Thornthwaite, 1931). Hence, the PE index indicates how much moisture is available for plant growth. We will add it in the revised manuscript as follows:
“The variation of ρ among catchments is associated with the climate condition (Melton, 1957; Madduma Bandara, 1974; Wang and Wu, 2013), which can be represented by measures such as the precipitation effectiveness index. This index is defined as the sum of the ratio of total monthly precipitation to total monthly evaporation over a year (Thornthwaite, 1931), with a higher value indicating more moisture available for plant growth.”
Comment: L26: Isn´t this the question? How can the topography be “given” if these apparent power laws can be found? There seems to be some organization, so I would suggest to mention this here.
Response: Thank you for pointing it out. It is general phenomenon that the drainage density varies with the source area for any river network. We will correct it in the revised manuscript as follows:
“On another note, the ‘rate’ at which L_{T} varies with A_{0} is determined by a given topography.”
Comment: L47: some concepts might need more explanation, what is valley transmissivity in this context?
Response: Thank you for pointing it out. We will add the physical expression and dimension of valley transmissivity in the revised manuscript as follows:
“γ is the exponent of a hypothetical power function between A and valley transmissivity T. Note that T is calculated as the product of subsurface crosssectional area and conductivity, which in turn is expressed in units of cubic length per time (Prancevic and Kirchner, 2019).”
Comment: L85: I think some simpler variable for cumulative length might be better for readability.
Response: Thank you for the advice. We will replace the original notation ‘Θ_{w}’ of the cumulative mean length with a new notion. The main body of the revised manuscript will be edited accordingly as follows:
Comment: L106: Some additional explanation for this equation would help the reader.
Response: Thank you for the suggestion. To resolve your concern, we will add supporting explanation on the concept of Eq. (7) in the revised manuscript.
Comment: L108: What is meant by such a tradeoff?
Response: In the context, a tradeoff means that mechanisms and processes forming river networks harmonize with each other under limited space, manifesting a relation between two scaling exponents, h and ε, as h + ε = 1. As we stated in the manuscript, the two scaling exponents, h and ε, are within the respective narrow ranges of h = 0.5 ~ 0.7 and ε = 0.4 ~ 0.46. Furthermore, the interdependence between the two exponents is found to be h + ε = 1. Given that, the tradeoff between h and ε can be understood by comparing two catchments (X and Y): if catchment X has a higher h value than catchment Y, then the ε value of X is smaller than that of Y, provided both exponents are within their respective narrow ranges.
To deliver the meaning of the tradeoff more clearly, the original L108 will be edited in the revised manuscript, as follows:
“h + ε = 1 (Maritan et al., 1996), which suggests a tradeoff between the two exponents by balancing each other within their respective ranges to form the catchment boundary within a confined space.”
Comment: L171: Is the analysis sensitive to the method used to derive flow direction?
Response: We expect no effect of flow direction algorithms on the analysis result. Our response is based on comprehensive references on diverse algorithms to extract flow directions. It is inevitable because no previous study on the ρ_{0}A_{p} relationship has been conducted by using the other method of extracting either single or multiple flow direction.
In general, when producing more accurate and smoother geometric patterns, algorithms to define multiple flow direction (i.e., flow of a cell drains into all downslope neighboring cells) are likely to be superior to single flow direction algorithm (i.e., drainage of a cell occurs in the steepest downslope direction) (Quinn et al., 1991; CostaCabral and Burges, 1994; Pan et al., 2004). However, hydrologic responses of a given catchment are almost no different between single and multiple flow directions (Wolock and McCabe Jr., 1995). This suggests that the accuracy of flow paths at the local scale is compromised as the scope expands to encompass the whole catchment.
Furthermore, we do strongly anticipate the insensitivity of scaling features to a way to define flow direction, based on the study of Paik (2011). The reference demonstrates that power law relationships in Hack’s law and the areaexceedance probability distribution are manifested in river networks extracted by the other single flow direction algorithm for determining flow direction – the Global Deterministic 8 method newly developed by Paik (2008).
Comment: L175: Reach or stream?
Response: Thank you for pointing it out. To further clarify, the following sentence will replace the original L175 in the revised manuscript:
“In NHDPlusV2, a channel forming area is given for stream channels at the most upstream points of individual flow paths in each river network.”
Comment: L188: You mention that you applied a linear fit. This refers to the logarithms? As one can also see in Fig. 2a this fit is only valid for a certain interval of values, therefore how did you derive the lower cutoff value?
Response: It seems that you were confused by the part ‘linfit’ in Matlab’s nlinfit function. The function is designed not for fitting a linear regression but for a nonlinear regression for a given dataset. Indeed, the nlinfit function aims to identify the bestfitting parameters by minimizing the sum of the squares of the residuals for a defined nonlinear model. In this study, we applied the nlinfit function to estimate the best parameters fitting Eqs. (17) and (18), respectively.
To clarify it, we will describe additional explanation in the revised manuscript as follows:
“To estimate the bestfitting parameters, we employed Matlab’s nlinfit function which is designed for nonlinear regression for a given dataset. The objective of the function is to minimize the sum of the squares of the residuals for a defined nonlinear model.”
Comment: L225: Fig. 2a: Different A_{0} values have been found, what might the range between the minimum A_{0} and the maximum A_{0} be related to?
Response 12: Source area A_{0} is intimately influenced by diverse conditions characterizing a catchment, such as hydrological, climatic, geomorphological, and geological conditions. Particularly, hydroclimatic conditions are significantly related to the source area and its corresponding drainage density. This relationship is well justified by the Abraham curve (Abrahams, 1984), which shows that drainage density of a river network decreases in arid regions and increases in humid regions. Thereby, A_{0} variability in our studied catchments indicates that this study encompassed a sufficient variety of environmental conditions affecting the formation of river networks. To demonstrate it in detail, we will conduct further analysis to identify potential relation between the statistics of source area for a given catchment and hydroclimatic condition such as annual precipitation.
Citation: https://doi.org/10.5194/hess2023271AC1

AC1: 'Reply on RC1', Kyungrock Paik, 16 Feb 2024

RC2: 'Comment on hess2023271', Massimiliano Schiavo, 15 Jan 2024
The comment was uploaded in the form of a supplement: https://hess.copernicus.org/preprints/hess2023271/hess2023271RC2supplement.pdf

AC2: 'Reply on RC2', Kyungrock Paik, 16 Feb 2024
We thank the Referee 2 for constructive comments. Our response to individual comments is listed below. We will reflect any change addressed below into the revised manuscript which is requested in the subsequent step.
Comment 1: L. 2025: the literaturebased explanation of the relationship between drainage density and climatic conditions of catchments can be elaborated a little bit more, e.g. providing further explanation about how L_{T} increases at the decrease of A_{0}, and vice versa.
Response 1: Thank you for your constructive comment. To clarify the intermediate mechanism, we will substitute the original L2125 with the following sentences in the revised manuscript:
“The variation of ρ among catchments is associated with the climate condition (Melton, 1957; Madduma Bandara, 1974; Wang and Wu, 2013), which can be represented by measures such as the precipitation effectiveness index. This index is defined as the sum of the ratio of total monthly precipitation to total monthly evaporation over a year (Thornthwaite, 1931), with a higher value indicating more moisture available for plant growth. A_{0 }reduces as the catchment becomes wetter, water accumulates more readily in the soils of lowgradient areas, and saturated areas expand accordingly. This mechanism leads to the enlargement of the stream network (greater L_{T}). Conversely, when the catchment gets drier, A_{0} increases, which in turn results in the contraction of the stream network (Godsey and Kirchner, 2014; Hooshyar et al., 2015; Durighetto et al., 2020).”
Comment 2: General comment: you could provide a 2panel Figure with sketched two kinds of basins for a prompt (visual) appraisal of the involved variables, in particular A_{0} and A_{p}. They can be drawings of two different catchments you investigated, or just two exemplary (not real) ones, or again a sketch representing two catchments one upstream and one downstream of the same river.
Response 2: Thank you very much for the constructive comment. In Sect. 2.2 of the revised manuscript, we will add a new figure that conceptually explains the core variables studied in this work (attached). This will be designated as Figure 1 in the revised manuscript.
Comment 3: Figure 1 and lines above: the usage of 5 climatic regions in the US is forgotten in the rest of the manuscript. It seems that the networks, which are then investigated in their PL behaviors, are not related to the climatic region they are located anymore after Figure 1. Thus, there are two choices: 1) you can make the US Figure smaller and surround it (e.g. in enlargements) with as many Figures as the catchments you provide in Figure S1, and referring to climatic areas to generally frame these catchments and their overall climatic (hence precipitation? Please specify) conditions; 2) (I suggest this one) you can make a more impactful use of these climatic regions. For example, you can make a boxplot chart by classifying all the apparent drainage densities by climatic region. See, as an example, the classification of salinity values for depositional environments done by Schiavo et al. (2023) (Figure 5). If you choose to go down for keeping the subdivision in climatic zones of the US territory and therefore this influence on networks’ structure, you should provide a classification of all the required exponents for all 5 climatic areas (see the following points).
Response 3: Thank you very much for your considerate comment. After deep consideration, we have decided to incorporate the climatic feature merely as one of explanatory variables for the study areas. (i.e., the first option you suggested). The primary reason for this decision is that a climatedependent classification of the analyzed power law exponents does not yield impactful findings. In other words, it demonstrates that river networks formed in different climate zones surprisingly reveal consistent attributes in describing fractal structures.
To provide more suitable guiding information at an earlier stage, we propose a new figure that presents the layouts and locations of the studied14 river networks, by merging Figure S1 with adjusted Figure 1 from the original manuscript (attached). This new figure will be designated as Figure 2 in the revised manuscript.
Comment 4: Figure S1: do all the DEM have the same spatial resolution? If yes, it’s ok; if not, you should homogenize the meshes before employing any routing algorithm. Rigon et al. (1996) and Maritan et al. (2002) correctly underlined the multiscaling problem when treating DEMs at different resolutions.
Response 4: Yes. All DEMs analyzed have the same resolution.
Comment 5: Figure 2a: it is not clear to me which are the values of the 3 thresholds employed for discretizing the xaxis, i.e. the pruning area variable; they seem to be about 0.02, 0.5, and 3 km^{2}, respectively. How did you choose them? Could you report the numerical values and motivations (if any) you employed for these choices?
Response 5: We’d like to remind that only one threshold of A_{p} was applied in the xaxis to characterize the ρ_{0}A_{p} relationship. The threshold A_{p} value is defined as a source area A_{0} for a given catchment. The value of A_{0 }is a constant determined as the median of channel forming areas extracted from the National Hydrography Dataset Plus Version 2 (NHDPlusV2). All relevant descriptions are given to Table 1, L174177 and L206207 of the original manuscript.
Despite the existing explanation in the manuscript, we fully anticipate that the original Figure 2a can confuse readers as you addressed those questions. To resolve the expected issue, we will revise the figure (attached), keeping only the vertical line for the median A_{0} (i.e., solid black line). It will be designated as Figure 3a in the revised manuscript. Note that we will insert a new threshold of drainage area which distinguishes hillslope region and unchannelized valley networks (i.e., dashed black line). Its relevant explanation is given in the next Response 6.
Comment 6: Figure 2a: it would be nice to further discuss the PL relationship of each trait of the 3 portions you identified, i.e. values of apparent density for A_{p}<0.02, 0.02<A_{p}<0.3, and A_{p}>0.3 (see e.g. D'odorico and Rigon, 2003). In other words, to offer results for the first, middle, and tail traits of the powerlaw relationship proposed in Figures 2a and 2b. This aspect would be very interesting to widen the investigation to the pruning areaapparent density relationships not only in the whole branch, but in each upstream, medium, and downstream trait. It also would provide a stronger confidence interval for estimating PL’exponents. Then, you’d come up with 4 exponents for each n=1,14 network, each exponent referring to a different portion of the network (like η1, η2, and η3), and the 4th for the “whole” one (η, you already did this).
Response 6: Linking to our Response 5, your Comment 6 is interpreted to inquire the existence of a power function even in the range of A_{p} < A_{0}. It is a very interesting and insightful topic to study further because the power function is highly likely to corroborate the selfsimilarity in flowpaths that can potentially appear in nonfluvial regions. Indeed, there are previous studies demonstrating that rill channels and valley networks show scaling relationships and geometric branching patterns found from river networks at equilibrium (e.g., Raff et al., 2004; Gangodagamage et al., 2011; Saybold et al., 2018).
Interestingly, the insight on the companion powerlaw had been also addressed during our work. It led us to apply a power function to an extent including unchannelized valley networks. The extent is defined by a lower limit A_{T}, i.e., a turnover area in the drainage arealocal slope curve and an upper limit A_{0}, i.e., a channel forming area (Montgomery & FoufoulaGeorgiou, 1993; Perera & Willgoose, 1998; McNamara et al. 2006). Given that, the companion powerlaw covering both nonfluvial and fluvial regions is expressed as : ρ_{a} = αA_{p}^{}^{ηv}, for A_{T} < A_{p} < A_{0}. In other words, the newly introduced powerlaw exponent η_{v} is calculated from Zone 2 illustrated in the new Fig. 3a (see Response 5).
Although preliminary analyses on the topic had been conducted, we had not made further progress in that direction after careful consideration at that time. The decision was primarily aimed at avoiding deviation from the main novel contribution of this study: the powerlaw relation between the apparent drainage density and the pruning area for river networks. Nevertheless, to acknowledge the referee’s insightful comment, we completed the analyses during the preparation of this Author Comment letter. New results on η_{v }and A_{T} will be summarized as a Table, which will be inserted as Table 1 of the revised manuscript.
Moreover, we will add a new Figure 3c in the revised manuscript to compare the scaling exponents η_{v}_{ }and η together for all our study sites (attached). They will be categorized by dominant climate conditions based on the KöppenGeiger climate classifications. Note that our study catchments are located in temperate or continental climate conditions.
The new analyses show that the A_{p}ρ_{a} dynamics in the extent covering unchannelized valley networks (i.e., Zone 2) reveal the scaling feature, i.e., ρ_{a} = αA_{p}^{}^{ηv}, as similarly found from the extent only for river networks (i.e., Zone 3). The powerlaw exponent η_{v} (0.59 ± 0.06 with R^{2} > 0.94 and MSE < 0.2) is greater than η (0.45 ± 0.02) for all 14 studied catchments. The finding suggests that the geometric hierarchy of networks formed in Zone 2 is more compact than that in Zone 3. The distinct difference is attributable to more bifurcation in flow paths of valley networks than that of river networks, related to geomorphic regime transitions (e.g., Gangodagamage et al., 2011).
Besides the new contents in Table 1 and Fig. 3c, we will add further discussion on scientific insights of the powerlaw exponent η_{v }in the perspectives of fractal dimension and relevance to dominant climate conditions.
Comment 7: Figure 2b: please comment on the impact of the threefold classification of PLs upon the variability of apparent density under varying pruning areas. These three kinds of exponents allow you to better investigate the behavior of the plot you give in this figure. Indeed, the 1st trait is constant, and the other two are sloping. Please comment adequately, also referring to each network portion’s total area and branching structure.
Response 7: Thank you for the careful comment. In fact, the threshold classification is not useful in describing details of the original Figure 2b (i.e., the normalized ρ_{0}A_{p }distribution). Because it is aimed to clearly demonstrate the power law behavior captured in the trunk part of the original Figure 2a.
To clarify the issue, we will add the explanation in the caption of the new Figure 3 stated in our Response 6 as follows:
“The xaxis range corresponds to the extent of Zone 3 marked in Fig. 3a.”
Comment 8: Figure 3: can you give the correlation coefficient of all those exponents obtained upon the (24) and those with the (25) to quantify the correlation discrepancies? Moreover, could you provide the other 3 panels of the same Figure by plotting η1, η2, and η3 in the same way (and with them correlation coefficients)?
Response 8: Thank you for your constructive comment. We will add the following sentence to provide the correlation coefficient between η values from Eqs. (24) and (25) in Sect. 4 of the revised manuscript:
“Although a high correlation coefficient of 0.95 emerges in the pairs of η values from Eqs. (24) and (25), the two mathematical expressions for η result a contrasting trend when compared with observed η values.”
Regarding the suggested extra panels, we recognize its importance and potential support when additional analyses will be carried out to identify the powerlaw features found in different extents of the ρ_{0}A_{p} distribution. As stated in our Response 6, the suggested aspect falls outside the scope of the current study after our careful consideration.
Comment 9: Flow routing method: the choice of the D8 flow method is not adequately supported. Indeed, it has been probably (implicitly?) chosen concerning the DInf or other ones because D8 guarantees the maximum energy (local) dissipation (e.g. in Schiavo et al. (2022)) by always following the steepest descent. Please clarify this point, eventually referring to Schiavo et al. (2022), for a complete thermodynamic framework (please elaborate on these concepts a little bit) of the processes you are investigating.
Response 9: Thank you for the valuable comment. Firstly, we’d like to note that the D8 method for allocating flow direction was selected not by us, but by the research team formulating the NHDplusV2 dataset. Unfortunately, any reasoning behind the choice of the D8 method is not justified in the NHDplusV2 User Guide (McKay et al., 2012). Given that, we believe that providing an exact answer to your comment is beyond the scope of our responsibility.
Nonetheless, here we’d like to response as best as we can based on our knowledge of flow direction extraction algorithms. In general, when producing more accurate and smoother geometric patterns, algorithms to define multiple flow direction (i.e., flow of a cell drains into all downslope neighboring cells) are likely to be superior to single flow direction algorithm (i.e., drainage of a cell occurs in the steepest downslope direction) (Quinn et al., 1991; CostaCabral and Burges, 1994; Pan et al., 2004). However, hydrologic responses of a given catchment are almost no different between single and multiple flow directions (Wolock and McCabe Jr., 1995). This suggests that the accuracy of flow paths at the local scale is compromised as the scope expands to encompass the whole catchment.
Furthermore, we do strongly anticipate the insensitivity of scaling features to a way to define flow direction, based on the study of Paik (2011). The reference demonstrates that power law relationships in Hack’s law and the areaexceedance probability distribution are manifested in river networks extracted by the other single flow direction algorithm for determining flow direction – the Global Deterministic 8 method newly developed by Paik (2008).
To clarify the suitability in the D8 method application, we will advance the original L171172 as the following sentences in the revised manuscript :
“Flow direction for each cell was assigned using the Deterministic 8 method (O'callaghan and Mark, 1984), which chosen by the NHDplusV2 research team (McKay et al., 2012). The flow direction extraction algorithm is underpinned by the principles of maximizing energy dissipation in surface water flow and minimizing power in groundwater flow (Schiavo et al., 2022). The individual cells in DEM were postprocessed to discard depression or sink cell.”

AC2: 'Reply on RC2', Kyungrock Paik, 16 Feb 2024