How Do Channel Network Expansions and Retractions Affect the Shape of Hillslope Flowpath Length Distributions?

Abstract

The length of hillslope flowpaths plays a crucial role in various hydrological, geomorphological and ecological processes. However, hillslope lengths change as the flowing network expands or contracts in response to precipitation. Currently, there is limited understanding of how channel network dynamics influence the distribution of hillslope flowpaths, with existing studies primarily focusing on the numerical analysis of specific case studies. In this work, we propose an analytical framework to characterize the general principles underlying variations in the hillslope length distribution as a function of the total length of the drainage network, $$L$$. Our model provides a closed‐form solution for the probability density function of the distance‐to‐channel conditioned on the total flowing network length, $$\omega (\ell ,L)$$. As a benchmark, the model has been applied to 15 headwater catchments with a non‐perennial river network, achieving reasonably good performances. The model captures the change in shape of the hillslope length distribution as the network expands, in particular the increased probability associated with shorter flowpaths. Empirical data and model results suggest the existence of some universal features in the way the hillslope length pdf responds to variations in the flowing network: longer paths are systematically more affected than shorter ones, and the overall fraction of impacted pathways decreases as the channel network lengthens. Besides the theoretical insights, our method offers a simple analytical means to determine the shape of the hillslope length distribution for varying drainage densities, even when the numerical evaluation of $$\omega (\ell ,L)$$ is challenging.

Introduction

Hillslope flowpaths identify the flow paths that connect each point in a watershed to the channel network via subsurface or surface flow. The practical calculation of hillslope flowpaths in catchment hydrology follows a standard procedure based on the identification, in a Digital Terrain Map, of the shortest path connecting any hillslope pixel to the first channelized site of the catchment along flow directions determined on a topographic basis. The spatial distribution of the hillslope flowpath length, hereafter termed hillslope length distribution or hillslope flowpath length distribution, represents the probability density function of the length of the hillslope flowpaths, calculated over all the unchanneled sites of the catchment.

The length of hillslope flowpaths is a key mathematical object that affects important hydrological, geomorphological and ecological processes taking place in river basins. In particular, the hillslope length distribution controls the hydrologic response of rivers, modulating the shape of the unit hydrograph and thus flood magnitude (e.g., Borga et al., 2007; Da Ros & Borga, 1997; Di Lazzaro et al., 2016; Grimaldi et al., 2012; Malagó et al., 2018; Piccolroaz et al., 2016; Volpi et al., 2013; White et al., 2003). The unit hydrograph, in fact, shares key features of the rescaled width function (Rinaldo et al., 1995), which coincides with the hillslope length distribution when the channel celerity is much larger than the hillslope velocities (Botter & Rinaldo, 2003; Di Lazzaro, 2009; D’Odorico & Rigon, 2003). Furthermore, hillslope lengths constrain the age distribution of the discharge at the catchment outlet, with notable implications for water quality and solute mobilization, retention and transport (e.g., Benettin et al., 2013; Benettin et al., 2022; Botter et al., 2010; Sloan et al., 2017; Wigington et al., 2005). In geomorphology, the characteristics of hillslope pathways are known to bear a fundamental impact on the sediment yield delivered to downstream water bodies and the corresponding rate of soil erosion (Cavalli et al., 2013; Ferrier & Perron, 2020; Grieve et al., 2016; Li et al., 2022; Liu et al., 2000; Montgomery & Dietrich, 1992; Willgoose et al., 1992). The length of hillslope flow paths also represents a critical factor that shapes many ecological dynamics of a watershed. By influencing the patterns of water movement and nutrient delivery, water flow paths determine the distribution of species along river corridors and the primary productivity of both aquatic and riparian ecosystems (e.g., Tetzlaff et al., 2007; Ward, 1989; Wiens, 2002).

Provided that any hillslope flow path terminates with a channelized site, hillslope lengths are crucially dependent on the extent and configuration of the river network. However, channel networks are inherently dynamic, that is, they do change their shape and extent depending on the wetness conditions in the contributing catchment (Blyth & Rodda, 1973; Botter et al., 2021; Day, 1978). This circumstance implies that hillslope flow paths do change in time as the network expands and contracts. Limited studies have investigated how the hillslope flowpath length distribution changes in response to variations in the total flowing network (Durighetto et al., 2020; Mutzner et al., 2016; Van Meerveld et al., 2019). Existing analyses revealed that channel network dynamics have a significant impact on the mean and the shape of the hillslope length distribution. The latter, in fact, appears to be highly skewed toward the shortest paths when the catchment is wet, though it becomes more uniform during dry conditions, when there is a higher probability for the longer flow paths. Yet, existing studies considered only specific case studies, in which the hillslope flowpath length distributions were numerically evaluated for different network configurations based on a digital terrain map. Instead, a general theory designed to identify how the hillslope flow path length distribution changes in shape and scale as the flowing network expands and retracts is missing.

In this paper, we propose an analytical framework to characterize how the hillslope length distribution changes as a function of the total flowing network length. The primary hypotheses of this study are as follows: (a) the problem can be addressed in a lumped manner, focusing on the statistics of the hillslope length and removing the explicit dependence on spatial coordinates; and (b) the mechanism by which hillslope lengths decrease as the flowing network expands or contract is unique and universally applicable across all river basins, allowing the same framework to be utilized for any catchment with a dynamic river network. These hypotheses will be tested combining analytical derivations and numerical analyses performed on a series of real world case studies.

The reminder of this paper is organized as follows. Section 2 presents the mathematical framework used to tackle the problem at hand, and derives a general closed-form expression for the hillslope length distribution conditional to the total flowing network length. Section 3 describes all the ancillary methods necessary for the model application and evaluation. Section 4 compares the analytical solution with empirical hillslope length distributions derived from the analysis of real world digital terrain maps. Section 5 discusses the value of the model, the perspective for its practical application and the main limitations of the approach. A set of conclusions closes then the paper.

Theoretical Framework and Model Development

In this section we set up the theoretical framework used to describe the variations in the hillslope length distribution caused by network expansions or retractions, the latter being quantified via changes in the total active network length. The framework is based on the formulation of a partial differential equation dictating the behavior of the hillslope length distribution conditional to the total flowing channel length $$L$$, indicated as $$\omega (\ell \vert L)$$—consistently with the notation frequently utilized in statistical analysis. In particular, $$\ell $$ represents the length of a generic hillslope flow path, defined as the distance from a given point in the catchment to the nearest channel site along the flow direction. Meanwhile, $$\omega $$ denotes the probability density function of $$\ell $$ across all sites within the focus catchment. The function $$\omega (\ell \vert L)$$ contains the family of hillslope length distributions resulting from different spatial configurations of the channel network within the catchment (i.e., for different values of the flowing channel length $$L$$). This means that, when $$L$$ is set as fixed (and equal to e.g., $${L}^{\ast }$$), $$\omega \left(\ell \vert {L}^{\ast }\right)$$ represents a Probability Density Function (PDF) in the $$\ell $$ domain (i.e., $$\int \nolimits_{0}^{\infty }\omega \left(\ell \vert {L}^{\ast }\right)d\ell =1$$). From an empirical perspective, $$\omega (\ell ,L)$$ can be calculated as the fraction of pixels of the DTM that generate a hillslope flow path within the range $$(\ell ,\ell +d\ell )$$ scaled to the quantity $$d\ell $$ (to convert the result into a probability density of a continuous random variable). This ratio needs to be calculated keeping the flowing network length as a constant, given that $$\omega $$ represents a conditional probability associated with a fixed length of the flow network. The calculation is than repeated for different values of $$L$$, since channel network dynamics induce temporal fluctuations in the length of the flowing stream network. When the function $$\omega (\ell \vert L)$$ is known, the temporal changes in the flowing network length in response to seasonal and event-based climatic fluctuations can be translated in to the corresponding temporal variations of the hillslope length distribution, for example, by calculating the time-variant PDFs $$\omega (\ell ,t)=\omega (\ell \vert L(t))$$. The latter equivalence suggests that $$\omega $$ can be implicitly viewed not only as a conditional PDF but also as a time-variant probability density of a stochastic process. Since the primary goal of this paper is to identify the probabilistic structure of $$\omega $$ by solving a Partial Differential Equation (Section 3.1) where $$L$$ is treated as a deterministic independent variable of the problem, we will use the notation $$\omega (\ell ,L)$$—instead of $$\omega (\ell \vert L)$$—to denote the hillslope length distribution throughout the remainder of this manuscript.

It should be noted that the framework developed in this paper is quite general and applies to any type of dynamic channel network, including dynamically fragmented networks that are frequently observed in small headwater catchments. However, for the sake of simplicity, we shall assume hereafter that network expansion starts from the catchment outlet, and consists of upstream movements of the channel head(s)—while network retraction consists of downstream movements of the channel head(s). Accordingly, all the river networks (theoretical and empirical) considered here are assumed to be the result of a suitable network extractions performed by applying a threshold on the contributing area (i.e., all the sites with contributing area larger than a given threshold are assumed to be channelized). Under the above assumption, the sites having a distance $$\ell $$ from a fixed and predefined channel network can be split into two parts, as indicated in Figure 1: (a) the subset of locations that belong to the contributing catchment of the network head(s) (red shading in Figure 1); and (b) the subset of points that do not belong to the contributing area of any network head (green shading in Figure 1). The distinction is key here, as the sites belonging to group A change their distance from the channel network when the heads move in response to network expansions or retractions, while the sites belonging to group B do not.

[IMAGE OMITTED. SEE PDF]

Crucially the proportion of sites whose distance-to-channel is impacted by the network dynamics changes with the network length $$L$$. In fact, when the network is short the distance from the channel changes in response to network expansions/contractions for the majority of the catchment locations. However, when the network is fully expanded, the hillslope length of many sites remains unaffected by network dynamics. Similarly, the number of sites in the catchment with a distance to the channel equal to $$\ell $$ varies with $$\ell $$. Consequently, the number of sites affected or unaffected by network dynamics, which are represented by the lengths of the red and green lines in Figure 1 also depends on the value of $$\ell $$. The fraction of sites with a hillslope flowpath length equal to $$\ell $$ belonging to group A (the impacted sites), scaled with the total number of sites having a distance-to-channel equal to $$\ell $$ (when the network length is $$L$$) is indicated as impact coefficient, and is denoted hereafter as $$\alpha (\ell ,L)$$ to highlight the fact that $$\alpha $$ is a function of both $$\ell $$ and $$L$$. Mathematically, the definition of the impact coefficient can be summarized as: 1 $$\alpha (\ell ,L)=\frac{\text{sites}\,\in \,\mathrm{A}\,\text{with}\,\mathrm{a}\,\text{hillslope}\,\text{flowpath}\,\text{of}\,\text{length}\,\ell \,(\text{network}\,\text{length}\,=\,L)}{\text{total}\,\text{number}\,\text{of}\,\text{sites}\,\text{with}\,\mathrm{a}\,\text{hillslope}\,\text{flowpath}\,\text{of}\,\text{length}\,\ell \,},$$ which in fact corresponds to the ratio between the length of the red line and the sum of the lengths of the red and dark green lines in Figure 1. According to the definition of $$\alpha $$, the fraction of pixels that are not affected by an increase of the network length (corresponding to group B)) is simply given by $$1-\alpha (\ell ,L)$$. Since it is, by definition, a relative fraction, the coefficient $$\alpha (\ell ,L)$$ must take values in the range [0,1].

The Case of a Single Channel

For the sake of convenience, we shall start describing the case in which the network is not branching, that is, there's a single channel head. Let us start considering what happens to the hillslope flowpath lengths when the network expands by an infinitesimal quantity $$dL$$. Based on our previous definitions we can identify two classes of sites in the catchment: the fraction that is impacted by the network expansion and the fraction that is not impacted by the network expansion. In particular, let us focus on the hillslope flowpath with a length $$\ell $$ before the network gets expanded and let us analyze how they respond to the considered network expansion. For a fraction of the sites originally at a distance $$\ell $$ from the network the hillslope length will not change in response to the network expansion. These sites are the ones that belong to group B—the amount of which is given by $$[1-\alpha (\ell ,L)]\omega (\ell ,L)$$. On the other hand, for a fraction $$\alpha (\ell ,L)$$ of the sites originally at a distance $$\ell $$ (the impacted ones, group A) the distance from the nearest channel is reduced by a quantity equal to the increase of the network length, $$dL$$. For similar reasons, there will be a fraction $$\alpha (\ell +dL,L)$$ of pixels originally characterized by a distance-to-channel equal to $$\ell +dL$$ for which, in response to the network expansion, the hillslope length is reduced by $$dL$$ (from $$\ell +dL$$ to $$\ell $$).

All this reasoning can be condensed in the following difference equation for $$\omega $$, which is essentially a local balance of probability: 2 $$\omega (\ell ,L+dL)=\omega (\ell ,L)+\alpha (\ell +dL,L)\omega (\ell +dL,L)-\alpha (\ell ,L)\omega (\ell ,L)\,.$$ In this equation, the probability density associated to the hillslope length $$\ell $$ when the network is expanded (length equal to $$L+dL$$) is the sum of three terms: the first term on the r.h.s. represents the PDF of $$\ell $$ before the network expansion, the second term represents the additional contribution given by the sites for which the distance changes from $$\ell +dL$$) to $$\ell $$, and the third term represents the loss of probability due to the sites for which, in response to the network expansion, the distance to channel changes from $$\ell $$ to $$\ell -dL$$.

Figure 2 graphically summarizes the mechanism though which the shape of the hillslope length distribution changes as the network expands. For a given $$\omega (\ell ,L)$$ (shown in Figure 2a), the impact coefficient defines, for every value of $$\ell $$, the fraction of $$\omega (\ell ,L)$$ that is affected by the network expansion (group A), as shown in Figures 2b and 2c. As the network length increases, the hillslope length distribution associated to the fraction of impacted sites experiences a leftward translation toward lower values of $$\ell $$ (Figure 2d). At the same time, the hillslope length distribution associated with the other sites (group B) remains unaffected. The combination of these two processes leads to the overall change of shape of $$\omega $$, as represented in Figure 2e.

[IMAGE OMITTED. SEE PDF]

Rearranging Equation 2 and dividing both sides by $$dL$$ yields to: 3 $\begin{align*}\hfill & \frac{\omega (\ell ,L+dL)-\omega (\ell ,L)}{dL}=\hfill \\ \hfill & \hspace{2.77695pt}\hspace{2.77695pt}=\frac{\alpha (\ell +dL,L)\,\omega (\ell +dL,L)-\alpha (\ell ,L)\,\omega (\ell ,L)}{dL}.\hfill \end{align*}$

By taking the limit for $$dL\to {0}^{+}$$ at both sides the following Partial Differential Equation (PDE) is obtained: 4 $$\frac{\partial }{\partial L}\omega (\ell ,L)=\frac{\partial }{\partial \ell }\left[\alpha (\ell ,L)\,\omega (\ell ,L)\right].$$

Equation 4 is a 1D linear first-order partial differential equation in conservation form. It is usually referred to as a convection (or advection) equation with non-constant coefficients, where the convection term corresponds to the impact coefficient $$\alpha $$, which is the focus of the following subsection.

Parametrizing $$\alpha $$

The equation derived in the previous section indicates that the impact coefficient is a key component of the problem under consideration. An empirical approach based on the analysis of the DTM can be a useful tool to get some insight on the possible shape of $$\alpha (\ell ,L)$$. This consists in analyzing data derived from several terrain maps in which the fraction of affected and unaffected pixels can be experimentally evaluated for different network configurations (and lengths). Specifically, for a given network configuration (i.e., for a given network length $$L$$), the distance of each hillslope pixel from the network is calculated based on the D8 criteria. Then, for each distinct value of the distance-to-channel $$\ell $$, the number of pixels belonging to the contributing catchment(s) of the source(s) is determined and divided by the total number of pixels at that distance. In this way, the empirical pattern of $$\alpha (\ell ,L)$$ is obtained. The analysis provides very useful insight on the behavior of $$\alpha $$, as the example reported in Figure 3 indicates. In particular, empirically derived impact coefficients suggest the existence of three emerging systematic patterns: (a) in spite of the pronounced scattering, $$\alpha $$ in general increases with $$\ell $$ since, for a given network, the short and the long flow paths are affected differently by network dynamics, with the long pathways generally being more impacted than the short ones; (b) the rate of increase of $$\alpha $$ decreases with $$\ell $$ (i.e., the longest pathways are impacted in a more homogeneous manner as compared to the short paths); (c) on average, $$\alpha $$ decreases with $$L$$ for any value of $$\ell $$, thereby implying that the overall impact of network expansions/contractions on the hillslope flowpath lengths decreases as the network becomes longer.

[IMAGE OMITTED. SEE PDF]

While it's not convenient to solve PDE Equation 4 using empirically derived impact coefficients of the type shown in Figure 3 because of the pronounced scattering of $$\alpha $$, we shall exploit the general features exhibited by the observed $$\alpha $$ to identify simple parametrizations of the impact coefficient that are keen to analytical treatments.

The idea pursued in this paper is to parametrize $$\alpha $$ based on purely geometrical arguments applied to simple and idealized morphological settings. To this aim, as a first approximation let us consider a symmetrical V-shaped catchment of infinite size as shown in Figure 4, where the hillslope pathways are parallel and symmetrical with respect to the central valley, where the channel is located. Moreover, they are assumed to converge toward the valley bottom with a uniform angle. On a planar view, the angle between the drainage directions and the horizontal line, perpendicular to projection of the valley bottom on the horizontal plane, is indicated as $$\delta $$. Without loss of generality, the vertical component of the flow paths is neglected.

[IMAGE OMITTED. SEE PDF]

In this setting, the network expands and retracts within the main channel at the bottom of the valley. Therefore, for every network (of arbitrary length $$L$$) one can identify the points having the same distance to the nearest channel, $$\ell $$, and calculate the fraction of these points which is impacted by network expansions. Figure 4 shows, for a given network configuration, the segment containing the points which are impacted by network expansion (with length $$4\ell \,\sin \left(4{5}^{o}\hspace{2.77695pt}-\hspace{2.77695pt}\delta /2\right)$$ and indicated in red) and that including all the points in which the hillslope path length remains the same when the network gets longer (with length $$2L$$, indicated in green). The total length of the curve is then $$2L+4\ell \,\sin \left(4{5}^{o}\hspace{2.77695pt}-\hspace{2.77695pt}\delta /2\right)$$ and the ratio of the impacted length over the total length which provides an estimate of the impact coefficient $$\alpha $$ can be expressed as: 5 $$\alpha (\ell ,L)=\frac{4\ell \,\sin \left(4{5}^{o}\hspace{2.77695pt}-\hspace{2.77695pt}\delta /2\right)}{2L+4\ell \,\sin \left(4{5}^{o}\hspace{2.77695pt}-\hspace{2.77695pt}\delta /2\right)}=\frac{\ell }{L\frac{1}{2\,\sin \left(4{5}^{o}\hspace{2.77695pt}-\hspace{2.77695pt}\delta /2\right)}\hspace{2.77695pt}+\hspace{2.77695pt}\ell }$$

Introducing the parameter $$a=\frac{1}{2\,\sin \left(4{5}^{o}\hspace{2.77695pt}-\hspace{2.77695pt}\delta /2\right)}$$, the previous equation leads to the following expression for the impact coefficient: 6 $$\alpha (\ell ,L)=\frac{\ell }{a\,L+\ell }\hspace{2.77695pt},$$ which constitutes the basis of our analytical hyperbolic model (in which $$\alpha $$ is an hyperbolic function of $$\ell $$ as Equation 6 indicates). Interestingly this model is able to guarantee the three main properties exhibited by empirically derived impact coefficients (namely, $$\partial \alpha /\partial \ell > 0$$, $${\partial }^{2}\alpha /\partial {\ell }^{2}< 0$$, $$\partial \alpha /\partial L< 0$$). Note also that $$a=1/\sqrt{2}$$ when $$\delta =0$$, while for positive values of $$\delta $$, $$a$$ increases significantly (e.g., by a factor of 2 when $$\delta =4{5}^{o}$$).

For real-world catchments, however, the value of $$a > 0$$ is very likely influenced by other morphological features not included in this simplified formulation, such as the shape of the catchment, the presence of converging hillslope pathways, and others. Therefore, we shall consider $$a$$ as a model parameter to be determined on an experimental basis.

The solution to PDE Equation 4 for a coefficient $$\alpha $$ of the kind given by Equation 6 and for the initial condition $$\omega (L=0,\ell )={\omega }_{0}(\ell )$$ is given by 7 $$\omega (\ell ,L)=\frac{aL+\ell }{\ell }{\left[\frac{\ell }{\ell +(1+a)L}\right]}^{\tfrac{a}{1+a}}\,{\omega }_{0}\left(\ell {\left[1+(1+a)\frac{L}{\ell }\right]}^{\tfrac{1}{1+a}}\right),$$ where $${\omega }_{0}$$ represents the hillslope length PDF when there is no flowing network (i.e., for $$L=0$$). The above equation will be referred to as the hyperbolic model and provides a simple equation for describing the change in the shape of the hillslope length distribution of a basin as the flowing network expands or retracts. The initial condition $${\omega }_{0}(\ell )$$ must be known for making Equation 6 fully explicit, and it can be specified as the empirically derived catchment width function or through any analytical approximation of it.

Generalizing the Single-Channel Model for Branching River Networks

In this section, we aim to generalize the framework to cases in which the network is branching. In particular, we first discuss the case in which the network contains more than one channel head and the network is assumed to branch only in corresponence of its head(s)—that is, no tributaries can be generated laterally from existing streams. In this case, we could repeat the same reasoning provided in the previous sections for a single channel, the main difference being that as the network expands, the infinitesimal flowing length increment $$dL$$ has to be split among all the active branches in the network. For the sake of simplicity, let us assume that all the channel heads active at a given time experience the same upstream expansion, that is, all of them get the same active length increment. This assumption does not necessarily reflect physical reality in all cases. However, it seems to be the most reasonable choice here to study the effect of branching on the PDE governing $$\omega (\ell ,L)$$, especially considering that the principle of maximum entropy suggests there should be no significant differences among the various heads of a river network. Under the above assumption, the classification of the catchment sites into the groups (a) and (b) previously identified (namely, the affected and unaffected sites) still applies. However, in this case, when the total active network length increases by a quantity equal to $$dL$$, the sites belonging to the contributing catchment of the channel heads would experience a decrease of their distance-to-channel equal to $$dL/n(L)$$, $$n(L)$$ being the number of active sources when the network has a length equal to $$L$$. In this case Equation 2 needs to be rewritten as follows: 8 $$\omega (\ell ,L+dL)=[1-\alpha (\ell ,L)]\omega (\ell ,L)+\alpha (\ell +dL/n,L)\omega (\ell +dL/n,L)$$

Rearranging terms in Equation 8, dividing both sides by $$dL$$ and dividing and multiplying the right hand side by $$n$$ gives the following equation: 9 $\begin{align*}\hfill & \frac{\omega (\ell ,L+dL)-\omega (\ell ,L)}{dL}=\hfill \\ \hfill & \hspace{2.77695pt}\hspace{2.77695pt}=\frac{1}{n}\frac{\alpha (\ell +dL/n,L)\,\omega (\ell +dL/n,L)-\alpha (\ell ,L)\,\omega (\ell ,L)}{dL/n}\hfill \end{align*}$

Finally, by taking the limit for $$dL\to {0}^{+}$$ at both sides it's possible to get the following PDE: 10 $$\frac{\partial }{\partial L}\omega (\ell ,L)=\frac{1}{n(L)}\frac{\partial }{\partial \ell }\left[\alpha (\ell ,L)\omega (\ell ,L)\right]$$

Which appears to be different from Equation 4, derived in the case of a single channel, as the convection factor is multiplied by the term $$1/n< 1$$. This would suggest that the convection term, that drives the speed at which the hillslope length distribution, changes in response to variations in network length and gets smaller and smaller as the network branches. However, in what follows we will demonstrate that this decrease in the convection term is compensated by a proportional increase in the impact coefficient $$\alpha $$, which is due to the fact that when the network gets branched its ability to impact a larger proportion of sites in the catchment increases as compared to the case of a single channel, because the network penetrates the landscape in a more efficient manner. To demonstrate the latter ansatz, we shall refer again to the simplified morphological setting introduced in Section 3.2, a symmetrical V-shaped catchment with flow directions perpendicular to the bottom valley ($$\delta ={0}^{o}$$), where a single bifurcation is however introduced along the main channel. In this simplified setting, the bifurcation is assumed to give rise to two symmetrical sub-valleys. Figure 5 shows the reference morphological scheme for branching networks, and gives some hints on how to compute on a purely geometrical basis the affected and unaffected lengths in this case (i.e., when the network has a single branching node).

[IMAGE OMITTED. SEE PDF]

For reasonably small values of the hillslope flowpath length $$\ell $$, the impact coefficient can be calculated based on the length of the segments $$\overline{AB}$$, $$\overline{BC}$$ and $$\overline{HD}$$ indicated in Figure 5b. The measures of these three segments can be expressed as: $$\overline{AB}=2\ell \,\sin \left(\pi /8\right)=\ell \sqrt{2-\sqrt{2}}$$, $$\overline{BC}=2\ell \,\cos \left(\pi /8\right)=\ell \sqrt{2+\sqrt{2}}$$ and $$\overline{HD}=\ell \sqrt{2}$$. The impacted length is then equal to $$2\left(\overline{AB}+\overline{BC}\right)=2\ell \left[\sqrt{2-\sqrt{2}}+\sqrt{2+\sqrt{2}}\right]$$, while the length of the segment containing the points which are not impacted by he network expansion is $$2L-2\,\overline{HD}=2L-2\ell \sqrt{2}$$. An estimate of the related impact coefficient is then given by the following equation: 11 $$\alpha (\ell ,L)=\frac{2\ell \left[\sqrt{2-\sqrt{2}}+\sqrt{2+\sqrt{2}}\right]}{2L-2\ell \sqrt{2}+2\ell \left[\sqrt{2-\sqrt{2}}+\sqrt{2+\sqrt{2}}\right]}$$

Which can be approximated as follows: 12 $$\alpha (\ell ,L)\approx \frac{2.1\,\ell }{L\frac{1}{1.2}\hspace{2.77695pt}+\hspace{2.77695pt}\ell }\approx 2\frac{\ell }{L/\sqrt{2}+\ell }$$

The impact coefficient for a network with two active branches is, therefore, approximately twice the impact coefficient in a V-shaped catchment with a single channel (Equation 5 calculated with $$\delta ={0}^{o}$$). This analysis can be easily extended to the case of multiple heads, showing that when $$n$$ sources are active the impact coefficient $$\alpha $$ is about $$n$$-times the hyperbolic coefficient derived for a single channel as given by Equation 6. Consequently, the factor $$1/n$$ which modulates the convection term in Equation 10 cancels out when we consider the proper shape of the impact coefficient in a network that branches in correspondence of the channel heads. Therefore, under the assumptions made the PDE determining the evolution of the hillslope length distribution as a function of the total network length $$L$$ is approximately the same as the one derived in Section 3.1 for a single channel, that is, Equation 4 with a convection term $$\alpha $$ expressed by Equation 6.

The general case of branched river networks, where new branches can emerge both at the channel head(s) and laterally from existing streams, presents significant challenges. In this scenario, the group A of pixels is no longer composed solely of those contributing to a channel head; it must also include all pixels that contribute to flow at sites within the active network where new branches are expected to be originated. While our approach could potentially address this complexity, providing a general analytical description of the case in which new branches emerge from existing streams seems infeasible. Therefore, we address this generalized case using a more practical approach, as detailed below. We proceed with the ansatz that the same model can be applied to describe all the hillslope length distributions, regardless of the degree (and type) of branching of the channel network. All the theoretical arguments available, in fact, suggest that the process by which the hillslope length distribution changes in response to network expansions and retractions does not significantly depend on the degree of branching of the network. We then verify our ansatz a posteriori, showing that the performance of our model does not correlate well with the branching features of the channel network (as shown in the Supporting Information S1).

Ancillary Material and Methods

Study Catchments

The hillslope length distributions were evaluated for 15 headwater catchments in Italy, Switzerland and North America, which are listed in Table 1. These proof-of-concept applications were used to test the performance of the analytical model and provide suitable initial conditions for the function $$\omega (\ell ,L=0)$$ as per Equation 7. The catchments were chosen for three main reasons: (a) empirical information on the temporal variations of the flowing length and the spatial configuration of the channel network is available; (b) the channel network is mostly temporary (i.e., the minimum observed drainage density is zero in most cases); (c) they have a relatively small contributing area $$A$$, a circumstance which makes the model application relatively cheap from the numerical point of view. The computation of $${\omega }_{e}(\ell ,L)$$ was based on a digital terrain map (DTM) of the area. While the actual spatial distribution of the flowing network in the selected case studies showed in most case complex and fragmented patterns (Botter et al., 2021), the active channel network was identified as the ensemble of pixels with a contributing area larger than a given threshold. For any given configuration of the channel network, the distribution of the distance-to-channel of all the hillslope pixels in each catchment was calculated along the water flowpaths, which were identified based on a standard D8 algorithm (O’Callaghan & Mark, 1984; Tarboton, 1997). The distribution of the hillslope flow path was calculated for different values of the threshold contributing area, so that the hillslope flowpath length distribution was obtained for different values of the total network length $$L$$, which correspond to different values of the drainage density $${D}_{d}=L/A$$. The considered networks are characterized by drainage density varying in between zero (when the outlet is the only site being channelized and the whole catchment is considered as a single hillslope) and a maximal drainage density corresponding to the empirically measured maximum drainage density (or a reasonable approximation of it). The relevant information of the study catchments is shown in Table 1. For more details the reader is referred to the quoted literature.

Table 1 Study Catchments and Their Relevant Properties

Catchment	Label	Area [pixel]	Area ($$\mathrm{k}{\mathrm{m}}^{2}$$)	Pixel ($$\mathrm{m}$$)	$${D}_{d,\max }$$ ($$\mathrm{k}{\mathrm{m}}^{-1}$$)	$${\ell }_{\max }$$ ($$\mathrm{k}\mathrm{m}$$)	References
Borsoia	BRS	551,461	13.79	5	0.56	7.86	–
Coweeta 12	CWT12	124,989	0.12	1	6.07	0.82	Jensen et al. (2017)
Coweeta 33	CWT33	306,933	0.31	1	3.40	1.51	Jensen et al. (2017)
Coweeta 40	CWT40	407,453	0.41	1	5.10	1.56	Jensen et al. (2017)
Fernow 14	FNW14	15,300	0.14	3	2.96	0.63	Jensen et al. (2017)
Fernow 16	FNW16	17,612	0.16	3	2.42	0.68	Jensen et al. (2017)
Fernow 37	FNW37	41,109	0.37	3	2.32	1.10	Jensen et al. (2017)
Focobon	FCB	21,043	8.42	20	2.24	6	–
Montecalvello	MNT	149,888	3.75	5	3.69	3.82	Durighetto et al. (2023)
Poverty 25	PVT25	249,506	0.25	1	5.11	1.08	Jensen et al. (2017)
Poverty 35	PVT35	426,424	0.35	0.90	5.32	1.53	Jensen et al. (2017)
Turbolo east	TURE	25,343	0.63	5	7.78	1.81	Senatore et al. (2021)
Turbolo west	TURW	18,367	0.46	5	3.50	1.38	Botter et al. (2021)
South Fork Potts	SFP	81,102	0.73	3	2.00	1.70	Jensen et al. (2017)
Rietholzbach	RTH	780,311	3.12	2	3.51	3.71	Durighetto et al. (2025)

Degree of Branching of the Network

The study catchments of Table 1 have different sizes and morphological characteristics (e.g., slope, elongation). Therefore, they represent a good benchmark for testing the hyperbolic model on a broad range of conditions. In particular, the way in which the network expands—starting from a single channel up to the fully expanded and branching configuration—differs significantly across the selected study catchments. Likewise, the rate at which the network branches does not remain constant during the network growth and is expected to change from catchment to catchment. For instance, some networks start branching at very short lengths but with a low branching rate, while other network can display a single channel for high values of the total flowing length and then they can branch relatively fast (see Figure 6). To properly assess the variability of the model performances under a variety of conditions, we have quantified the heterogeneity of the branching degree of the different networks considered in this paper. Figure 6 shows how the number of branches of the network ($${n}_{r}$$) increases as the flowing network gets longer, for a random selection of case studies. The plot indicates that the degree of branching of the considered networks is highly heterogeneous, and the branching rate changes considerably with $$L$$ and depends on the specific case study. The normalized branching indexes used in the y-axis of the Figure and its non-normalized version were used to analyze the correlation between the model performance and the degree of branching of the network (if any).

[IMAGE OMITTED. SEE PDF]

Model Evaluation

The performance of the model is evaluated through a comparison between the analytical solution given by Equation 7 and the empirical hillslope length distributions obtained from the analysis of the DTM. For every value of the channel network length $$L$$, the difference between the hyperbolic modeled $$\omega (\ell ,L)$$ and its empirical counterpart $${\omega }_{e}(\ell ,L)$$ is quantified by a scalar integral error, $$\varepsilon $$, defined as (Botter et al., 2008): 13 $$\varepsilon (L)=\frac{1}{2}\int \nolimits_{0}^{+\infty }\vert \omega (\ell ,L)-{\omega }_{e}(\ell ,L)\vert d\ell ,$$ where the factor $$1/2$$ is used to normalize the error in the interval $$[0,1]$$ (i.e., two completely disjointed PDFs would originate a unit integral error $$\varepsilon =1$$, while if, e.g., $$\varepsilon =0.1$$ then the modeled and the empirical hillslope length distributions do overlap for $$90\%$$ of their area). The value of the integral in Equation 13 is dependent on the flowing network length $$L$$, so the error of the hyperbolic model given in Equation 13 should be seen as a function of $$L$$. An overall estimate of the model performance, over the whole range of network lengths, can be calculated by averaging $$\varepsilon (L)$$ over $$L$$ in the range from $$L=0$$ (no flowing network) to $${L}_{\max }$$ (the maximum extent of the flowing network in each study catchment, as derived from observational studies): 14 $$\bar{\varepsilon }=\frac{1}{{L}_{\max }}\int \nolimits_{0}^{{L}_{\max }}\varepsilon (L)\,dL$$

The average error $$\bar{\varepsilon }$$, as given by Equation 14, is used in this paper to quantify the overall model performance across the range of network length observed in each study catchment. Accordingly, $$\bar{\varepsilon }$$ is minimized to select the best value of the parameter $$a$$ of the hyperbolic model, $${a}_{\mathit{opt}}$$. The errors defined by the indexes $$\varepsilon (L)$$ and $$\bar{\varepsilon }$$, however, do not provide any information about how the shape of the hillslope length distribution is captured by the model. Therefore, a visual comparison of $$\omega (\ell ,L)$$ and $${\omega }_{e}(\ell ,L)$$ will be also provided for some reference value of the total flowing length $$L$$. This visual comparison will be used to check if the modeled PDFs capture the key features of the empirical hillslope length distributions. We believe that the definition of other and more formal performance metrics is not necessary here, as the main purpose of the paper is to describe the analytical model and propose a series of proof-of-concept applications.

Results

Figures 7–9 show a comparison between modeled and empirical hillslope length distributions for some reference values of the drainage density $${D}_{d}$$. The Figures show three case studies with various degrees of performance (based on the error function defined by Equation 14): (a) the Fernow 37, for which the model error is the lowest; (b) the Coweeta 33, for which the model error is the highest and (c) the Poverty 25, for which the model error is similar to the mean of the errors across all the study catchments. The case studies not included in this Figures are reported in the SI, and a summary of the performance for all the considered catchments is given in Table 2.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Table 2 Summary of the Model Performances for the Various Catchments Introduced in Table 1

Catchment	$$a={a}_{\mathit{opt}}$$	$$a=1$$	$$\bar{{n}_{r}}$$	$$\bar{B}$$
$$\bar{\varepsilon }$$ (%)	$${\varepsilon }_{\max }$$ (%)	$$\bar{\varepsilon }$$ (%)	$${\varepsilon }_{\max }$$ (%)
Borsoia	17.2	19.6	20.6	25.1	67.18	4.37
Coweeta 12	9.6	14.2	9.6	14.9	1.86	4.81
Coweeta 33	20.6	31.1	23.7	34.4	3.29	7.17
Coweeta 40	12.6	17.1	20.0	34.2	3.17	4.42
Fernow 14	12.9	23.9	13.5	25.2	1.04	7.28
Fernow 16	14.4	23.3	15.2	25.0	1.13	10.54
Fernow 37	6.9	11.8	7.6	13.0	2.44	6.32
Focobon	23.0	30.8	23.2	31.2	16.82	2.92
Montecalvello	11.4	21.4	12.7	21.2	15.96	3.40
Poverty 25	8.4	11.9	15.0	26.9	3.13	4.01
Poverty 35	13.2	18.5	19.9	31.7	2.60	2.90
Turbolo east	17.5	23.6	18.9	27.3	6.52	4.74
Turbolo west	11.7	17.7	15.7	23.1	3.66	3.92
South Fork Potts	7.3	14.7	8.8	16.1	2.13	3.43
Rietholzbach	10.0	17.8	10.3	18.7	15.81	4.10

In the Figures, the empirical hillslope length distributions are compared to the model for two values of the parameter $$a$$, namely the default value $$a=1$$ and the optimized value $${a}_{\mathit{opt}}$$ that results from the minimization of the average error as defined by Equation 14. The empirical curves have been smoothed to remove the spurious fluctuations induced by the discretization of the catchment domain in the DTM. Nevertheless, the original version of $${\omega }_{e}$$ are shown in the background as light colors. In the comparison, all the modeled distributions have been obtained by imposing a smoothed version of the empirical hillslope length distribution for $$L=0$$ as initial condition. In all cases, the hyperbolic model correctly captures the increasing likelihood of shorter $$\ell $$ values emerging in the data as drainage density increases and the stream network expands. Overall, the hyperbolic model seems to be able to capture the global change in the shape of the hillslope length distribution due to the dynamics of the channel network within the catchment, and the transition toward a monotonically decreasing PDF observed when the drainage density increases.

Figure 10 shows a synthesis of the model performances. In particular, we show the patterns of the error functions $$\varepsilon (L)$$ defined in Equation 13 obtained with the hyperbolic model (Equation 7) for the three reference study catchments (the same shown in the previous figures). The two curves shown on each plot refer to different values of the parameter $$a$$: the calibrated value $${a}_{\mathit{opt}}$$, and the default value $$a=1$$. The curve corresponding to the optimal value $${a}_{\mathit{opt}}$$ should be seen as the best performance achievable with the hyperbolic model when attempting to model the family of hillslope length distributions obtained for the whole range of drainage densities investigated in the paper. For this reason, $${a}_{\mathit{opt}}$$ is not necessarily the best-performing value for each single configuration and length of the active network. In all cases, the initial condition $${\omega }_{0}$$ imposed in the hyperbolic model in Equation 7 consists in the hillslope length distribution at $$L=0$$ computed from the DTM—thereby implying that the error is by definition equal to zero for $$L=0$$. For the sake of representation, the empirical distributions have been smoothed before the computation of the error in order to remove the spurious fluctuations produced by the discretization of the catchment area into squared cells. This operation allowed us to focus on the general modifications of the hillslope length distribution produced by changes in the flowing length $$L$$. The errors obtained with the calibrated value of $$a$$ and the ones given by the default choice $$a=1$$ are very close to each other in most cases (with the exceptions of the catchments Coweeta 12, Coweeta 40, Poverty 25 and Poverty 35), with average errors $$\bar{\varepsilon }$$ in the range $$5\%\mbox{--}25\%$$ depending on the case study. In general, the error is higher for larger values of the drainage density, and the maximum values of $$\varepsilon $$ seldom exceed $$30\%$$. These results suggest that the performance of the model is by and large satisfactory. In particular, the hyperbolic model with $$a=1$$ can be used to get a reasonable first order approximation of the hillslope length distribution without any calibration (beside the specification of the distribution for a fixed and known value of $$L$$). Nevertheless, calibrating $$a$$ can improve the model performance most of the times, especially for relatively high values of the drainage density.

[IMAGE OMITTED. SEE PDF]

The evidence that the error functions flatten or decrease for the largest values of the drainage density—jointly with the poor correlation between the mean error and the mean branching degree of the network—indicates that the hyperbolic model can be used both for single-channels and highly branched networks, that is, the hillslope length distribution is not impacted by the degree of branching of the network but only depends on the total flowing length $$L$$.

Discussion

Given the complexity of the process examined and the simplicity of the analytical formulation presented in the paper, the model performance across the case studies analyzed is deemed satisfactory, with a mean error of the order of $$10\%$$ (i.e., the modeled and the observed hillslope length distributions do overlap for $$90\,\%$$ of their area). The error slightly increases when the parameter $$a$$ is not calibrated and set equal to the default value of 1, though remaining below $$20\%$$ in most cases. In particular, the model explains why, when the network is long enough, the mode of $$\omega $$ approaches zero, and the hillslope length PDF becomes nearly monotonically decreasing (i.e., the higher the hillslope length, the lower the probability density), as noted by other empirical studies (Durighetto et al., 2020; Mutzner et al., 2016; Van Meerveld et al., 2019). From the mathematical viewpoint, this change of shape of the hillslope length PDF in response to an increase of the network length $$L$$ is primarily driven by the power-law term of Equation 7, which has a monotonically decreasing behavior and tends to infinity for $$\ell \to 0$$. From the physical view point, this term reflects how flowpaths of various lengths are impacted in a different manner by network expansions and contractions (see below).

One of the key features of the model is the lumped nature of the approach, in which the spatial characteristics of the catchment (morphology, patterns of contributing area, configuration of the active network, etc) are surrogated by the variables $$\ell $$ and $$L$$, and the corresponding impact coefficient $$\alpha $$. The impact coefficient, in particular, is responsible for quantifying to what extent flowpaths with different length are impacted by the underlying network dynamics. The lumped formulation proposed in this paper is not designed to capture the response of each flow path to temporal variations in the flowing length, but can only capture the patterns emerging at the catchment scale when key statistics among all the active hillslope flowpaths are considered. From this perspective, the good performance of the model indicates that while the effect of network expansion and contraction of different hillslope flowpaths could be quite heterogeneous in space (some flowpaths are heavily impacted by changes in the flowing length, others are not), systematic patterns still arise at the basin scale, with the longest flow paths that are consistently more impacted by network dynamics than the shorter ones. This suggests that there's a stronger tendency for longer flow paths to shift toward shorter lengths, resulting in a sizable change in the shape of $$\omega $$ rather than a mere leftward shift of the probability density function. Additionally, the overall impact diminishes as network length increases, indicating that the rate at which the hillslope length distribution adapts to network expansion/contraction decreases with longer flowing networks. The presence of these common characteristics in the shape of the impact coefficient allowed us to effectively describe the system focusing on the statistical properties of the flowpath lengths in the $$(\ell ,L)$$ domain.

The main model deficiencies emerging from the comparison between the empirical and modeled hillslope length distributions are the following: (a) the model in some cases has some spurious fluctuations which generate multi modal PDFs that are not observed in their empirical counterparts; (b) the model tends to overestimate the probability associated to the shortest paths (with a length close to zero). This suggests that eliminating the spatial dimension of the problem (e.g., neglecting the spatial configuration of the active network to focus only on its total length) and parametrizing the impact coefficient using the hyperbolic model imply in some cases a misrepresentation of important features of $$\omega (\ell ,L)$$ (e.g., the modal probability density or the precise position of the mode). The primary cause of the spurious fluctuations in the modeled distributions is likely the purely convective nature of the PDE Equation 4, which lacks an inherent smoothing effect. If the catchment's width function $${\omega }_{0}(\ell )$$ contains significant peaks or irregularities, these features are likely to be preserved (with only minor changes in their shape) in the distributions $$\omega (\ell ,L)$$ as the river length $$L$$ increases. In fact, there is often a noticeable correspondence between the number of spurious fluctuations in the width function and those observed in the modeled distributions, particularly at high values of drainage density. More research is needed to identify how these limitations can be overcome within the analytical framework proposed in this paper.

The poor correlation of the model performance with key morphometric indexes and the degree of branching of the network (see Supporting Information S1) hints at the existence of a general mechanism that regulates the changes in the hillslope length distribution induced by variations of the flowing length. This mechanism seems to be independent on the spatial structure of the river network and the morphology of the basin. Accordingly, the standard hyperbolic model with $$a=1$$ generally performs quite well without requiring calibration. This indicates that as the active network expands or contracts, the spatial configuration of the flowing network evolves, influencing the statistics of hillslope lengths in a manner i.e., largely consistent across all river basins, regardless of the specific shape of the network or the morphometric characteristics of the landscape. The alteration in the hillslope length distribution is modulated by the power-law terms in Equation 7, which enhance the probability of the shortest flow paths as $$L$$ increases, rendering $$\omega (\ell ,L)$$ nearly independent of the initial distribution $${\omega }_{0}$$. The presence of these power-law terms in the analytical solution for $$\omega $$ challenges the use of an exponentially decreasing function to approximate the hillslope length distribution, as suggested by previous studies (e.g., Rodríguez-Iturbe & Valdés, 1979; Skaugen & Onof, 2014; Tucker et al., 2001). However, we lack definitive evidence to claim that this discrepancy reflects a flaw in the analytical model, particularly regarding the representation of the convective term $$\alpha $$ in Equation 4. Rather, this may suggest that the hillslope length distribution does not exhibit an exponential form but is instead monotonically decreasing, allowing it to be effectively approximated by an exponential function in some cases.

A potential limitation of the analysis presented in this paper is that only stream networks extracted with a threshold on the contributing area were considered. Thus, the effect of network disconnections and heterogeneous spatial patterns of node persistency (Durighetto et al., 2022) on the hillslope length distribution still remains to be explored. This important task is the object of ongoing research.

In conclusion, we emphasize that improving our ability to predict or analytically characterize how hillslope lengths vary with the active drainage density has significant implications for both hydrology and geomorphology. The pervasive nature of network dynamics poses a call for hydrology to integrate the expansion and contraction of flowing networks into existing rainfall-runoff models. When a catchment is wet, not only do water velocities in surface and subsurface domains increase, but the distance to the channel also decreases significantly, enhancing the capacity of river basins to rapidly transport water and mass downstream (Skaugen et al., 2023). Consequently, lumped, IUH-based formulations that employ a time-invariant transfer function linked to the rescaled width function (i.e., the function $$\omega \left({u}_{h}\,t\,,L\right)$$, where $$t$$ is time, $$L$$ the (constant) network length and $${u}_{h}$$ is a suitable water velocity in the hillslope) may offer a biased representation of critical hydrological processes due to an inaccurate estimation of the relevant hillslope length scales. This is especially true during high-flow conditions, when the network is extended and hillslope flow paths are shorter, or during droughts, when the network contracts and distances to the channel increase significantly. From this perspective, this paper represents a seminal step toward rethinking classical hydrological theories based on static networks, such as the geomorphological theory of the hydrological response and the width function-based unit hydrograph. A more precise characterization of the distribution of hillslope lengths is also crucial for geomorphological sciences, particularly for quantifying sediment transport from the uplands to flowing channels and for developing hillslope-channel connectivity metrics that reflect the temporal changes in hillslope lengths, as influenced by variations in catchment water storage (Cavalli et al., 2013; Heckmann et al., 2018). The closed-form solutions presented here can boost the development of analytical tools linking hillslope length statistics with streamflow patterns, and facilitate the practical estimation of flowpath length distributions in situations where numerical analysis may be challenging. e.g., in large basins, numerically evaluating the function $$\omega (\ell ,L)$$ can be computationally demanding. Additionally, there are circumstances in which the resolution of the Digital Terrain Map is too coarse to accurately capture the shape of the hillslope length distribution at high drainage densities, as the relative number of hillslope pixels may become comparable to that of channelized pixels. In all these circumstances, our simple analytical model could offer a valuable tool for hydrological and geomorphological studies.

Conclusions

This paper proposes a novel analytical model for the computation of the hillslope flowpath length as a function of the total active drainage density in a catchment. The model consists of a closed-form analytical solution for the PDF of the distance to channel conditional to the total flowing network length, $$\omega (\ell ,L)$$. As a benchmark, the model has been applied to 15 headwater catchments in which the network is non-perennial using the observed width function as a boundary condition for $$\omega (\ell ,L=0)$$. The application suggested the following conclusions:

• The model shows reasonably good performances: the mean relative error is in the range from 5 to 25$$\%$$, and the maximum errors—typically observed for drainage densities in between 0.5 and 3 $$\mathrm{k}{\mathrm{m}}^{-1}$$—remain below $$30\%$$ in most cases. Furthermore, the model is able to capture the change of shape exhibited by the hillslope flowpath length when the network expands, in particular the increase of probability associated to the shortest flowpaths, and the monotonically decreasing pattern exhibited by the function $$\omega (\ell ,L)$$ when the network becomes fully expanded and $$L$$ is sufficiently large.

• Our results indicate that the shape of the hillslope length distribution is modulated by the impact coefficient $$\alpha $$, which quantifies how flowpaths of different length are impacted by network expansion and contraction. In particular, empirical data and model performance evaluation both indicate that the fraction of pathways i.e., impacted by active network changes increases systematically with the length of the pathways (i.e., longer paths are more impacted than shorter paths). Moreover, as the channel network gets longer, the overall fraction of impacted pathways decreases.

• In spite of the observed heterogeneity of important morphometric features across the study catchments, our analysis indicates that the mechanism through which the hillslope length distribution changes in response to variations of the flowing network is somewhat “universal” (i.e., it can be quantified by a single law), as expressed mathematically in the power law terms that appear in Equation 7.

• Besides the theoretical insight, the practical advantages of our analytical model are manyfold. In particular, our framework allows one to estimate how in catchments with non-prennial streams the hillslope length distribution changes (starting from a given initial condition) as a function of the active drainage density using a simple analytical equation. This is possible also in cases in which the numerical calculation of $$\omega (\ell ,L)$$ from a DTM is problematic because it is computationally too expensive (e.g., for large river basins) or inaccurate (e.g., for small basins in which the resolution of the DTM is comparable to the mean hillslope lengths).

• The main limitations of the model include a poor ability to capture the number, the position and the amplitude (i.e., the $$\omega $$ values) of the modes of the hillslope length distribution (especially for high values of the drainage density), and spurious fluctuations that appear in the analytical $$\omega $$ function in some of the case studies. Further research is needed to identify how the model can be improved and extended to cases in which the active river network is dynamically fragmented.

Acknowledgments

This research was supported by the European Community's Horizon 2020 Excellent Science Programme (Grant H2020-EU.1.1.-770999). We thank the associate editor and two anonymous reviewers for the time and effort spent in handling our manuscript. Open access publishing facilitated by Universita degli Studi di Padova, as part of the Wiley - CRUI-CARE agreement.

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Data Availability Statement

The DTM input data of the analyzed catchments and the MATLAB scripts used for the analyses reported in this paper are available at (Cenzon et al., 2024) (Research Data Unipd) via with Creative Commons: Attribution 4.0 license.

References

Benettin, P., Rinaldo, A., & Botter, G. (2013). Kinematics of age mixing in advection‐dispersion models: Age mixing in advection‐dispersion models. Water Resources Research, 49(12), 8539–8551. https://doi.org/10.1002/2013WR014708

Benettin, P., Rodriguez, N. B., Sprenger, M., Kim, M., Klaus, J., Harman, C. J., et al. (2022). Transit time estimation in catchments: Recent developments and future directions. Water Resources Research, 58(11), e2022WR033096. https://doi.org/10.1029/2022WR033096

Blyth, K., & Rodda, J. C. (1973). A stream length study. Water Resources Research, 9(5), 1454–1461. https://doi.org/10.1029/WR009i005p01454

Borga, M., Boscolo, P., Zanon, F., & Sangati, M. (2007). Hydrometeorological analysis of the 29 August 2003 flash flood in the Eastern Italian alps. Journal of Hydrometeorology, 8(5), 1049–1067. https://doi.org/10.1175/JHM593.1

Botter, G., Bertuzzo, E., & Rinaldo, A. (2010). Transport in the hydrologic response: Travel time distributions, soil moisture dynamics, and the old water paradox. Water Resources Research, 46(3), 2009WR008371. https://doi.org/10.1029/2009WR008371

Botter, G., & Rinaldo, A. (2003). Scale effect on geomorphologic and kinematic dispersion. Water Resources Research, 39(10), 2003WR002154. https://doi.org/10.1029/2003WR002154

Botter, G., Vingiani, F., Senatore, A., Jensen, C., Weiler, M., McGuire, K., et al. (2021). Hierarchical climate‐driven dynamics of the active channel length in temporary streams. Scientific Reports, 11(1), 21503. https://doi.org/10.1038/s41598‐021‐00922‐2

Botter, G., Zanardo, S., Porporato, A., Rodriguez‐Iturbe, I., & Rinaldo, A. (2008). Ecohydrological model of flow duration curves and annual minima. Water Resources Research, 44(8), 2008WR006814. https://doi.org/10.1029/2008WR006814

Cavalli, M., Trevisani, S., Comiti, F., & Marchi, L. (2013). Geomorphometric assessment of spatial sediment connectivity in small Alpine catchments. Geomorphology, 188, 31–41. https://doi.org/10.1016/j.geomorph.2012.05.007

Cenzon, A., Botter, G., & Durighetto, N. (2024). Empirical hillslope length distributions and hyperbolic model [Dataset]. Centro di Ateneo per le Biblioteche dell'Università degli Studi di Padova. https://doi.org/10.25430/researchdata.cab.unipd.it.00001433

Da Ros, D., & Borga, M. (1997). Use of digital elevation model data for the derivation of the geomorphological instantaneous unit hydrograph. Hydrological Processes, 11(1), 13–33. https://doi.org/10.1002/(SICI)1099‐1085(199701)11:1⟨13::AID‐HYP400⟩3.0.CO;2‐M

Day, D. G. (1978). Drainage density changes during rainfall. Earth Surface Processes, 3(3), 319–326. https://doi.org/10.1002/esp.3290030310

Di Lazzaro, M. (2009). Regional analysis of storm hydrographs in the rescaled width function framework. Journal of Hydrology, 373(3–4), 352–365. https://doi.org/10.1016/j.jhydrol.2009.04.027

Di Lazzaro, M., Zarlenga, A., & Volpi, E. (2016). Understanding the relative role of dispersion mechanisms across basin scales. Advances in Water Resources, 91, 23–36. https://doi.org/10.1016/j.advwatres.2016.03.003

D'Odorico, P., & Rigon, R. (2003). Hillslope and channel contributions to the hydrologic response. Water Resources Research, 39(5). https://doi.org/10.1029/2002WR001708

Durighetto, N., Bertassello, L. E., & Botter, G. (2022). Eco‐hydrological modelling of channel network Dynamics—Part 1: Stochastic simulation of active stream expansion and retraction. Royal Society Open Science, 9(11). https://doi.org/10.1098/rsos.220944

Durighetto, N., Britt, R., Schirmer, M., & Botter, G. (2025). Variability in the shape of the active length‐streamflow relationship in temporary streams: Insights from an empirical analysis. Water Resources Research, 61. https://doi.org/10.1029/2024WR038732

Durighetto, N., Noto, S., Tauro, F., Grimaldi, S., & Botter, G. (2023). Integrating spatially‐and temporally‐heterogeneous data on river network dynamics using graph theory. iScience, 26(8), 107417. https://doi.org/10.1016/j.isci.2023.107417

Durighetto, N., Vingiani, F., Bertassello, L., Camporese, M., & Botter, G. (2020). Intraseasonal drainage network dynamics in a headwater catchment of the Italian alps. Water Resources Research, 56(4), e2019WR025563. https://doi.org/10.1029/2019WR025563

Ferrier, K. L., & Perron, J. T. (2020). The importance of Hillslope scale in responses of chemical erosion rate to changes in tectonics and climate. Journal of Geophysical Research: Earth Surface, 125(9), e2020JF005562. https://doi.org/10.1029/2020JF005562

Grieve, S. W., Mudd, S. M., & Hurst, M. D. (2016). How long is a hillslope? Earth Surface Processes and Landforms, 41(8), 1039–1054. https://doi.org/10.1002/esp.3884

Grimaldi, S., Petroselli, A., & Nardi, F. (2012). A parsimonious geomorphological unit hydrograph for rainfall–runoff modelling in small ungauged basins. Hydrological Sciences Journal, 57(1), 73–83. https://doi.org/10.1080/02626667.2011.636045

Heckmann, T., Cavalli, M., Cerdan, O., Foerster, S., Javaux, M., Lode, E., et al. (2018). Indices of sediment connectivity: Opportunities, challenges and limitations. Earth‐Science Reviews, 187, 77–108. https://doi.org/10.1016/j.earscirev.2018.08.004

Jensen, C. K., McGuire, K. J., & Prince, P. S. (2017). Headwater stream length dynamics across four physiographic provinces of the Appalachian highlands. Hydrological Processes, 31(19), 3350–3363. https://doi.org/10.1002/hyp.11259

Li, L., Nearing, M. A., Heilman, P., Nichols, M. H., Guertin, D., & Williams, C. (2022). Rangeland hillslope lengths: A case study at the walnut Gulch experimental Watershed, Southeastern Arizona. International Soil and Water Conservation Research, 10(4), 597–609. https://doi.org/10.1016/j.iswcr.2022.02.004

Liu, B. Y., Nearing, M. A., Shi, P. J., & Jia, Z. W. (2000). Slope length effects on soil loss for steep slopes. Soil Science Society of America Journal, 64(5), 1759–1763. https://doi.org/10.2136/sssaj2000.6451759x

Malagò, A., Vigiak, O., Bouraoui, F., Pagliero, L., & Franchini, M. (2018). The hillslope length impact on SWAT streamflow prediction in large basins. Journal of Environmental Informatics, 32(2), 82–97. https://doi.org/10.3808/jei.201700384

Montgomery, D. R., & Dietrich, W. E. (1992). Channel initiation and the problem of landscape scale. Science, 255(5046), 826–830. https://doi.org/10.1126/science.255.5046.826

Mutzner, R., Tarolli, P., Sofia, G., Parlange, M. B., & Rinaldo, A. (2016). Field study on drainage densities and rescaled width functions in a high‐altitude alpine catchment. Hydrological Processes, 30(13), 2138–2152. https://doi.org/10.1002/hyp.10783

O'Callaghan, J. F., & Mark, D. M. (1984). The extraction of drainage networks from digital elevation data. Computer Vision, Graphics, and Image Processing, 28. https://doi.org/10.1016/S0734‐189X(84)80011‐0

Piccolroaz, S., Di Lazzaro, M., Zarlenga, A., Majone, B., Bellin, A., & Fiori, A. (2016). HYPERstream: A multi‐scale framework for streamflow routing in large‐scale hydrological model. Hydrology and Earth System Sciences, 20(5), 2047–2061. https://doi.org/10.5194/hess‐20‐2047‐2016

Rinaldo, A., Vogel, G. K., Rigon, R., & Rodriguez‐Iturbe, I. (1995). Can one gauge the shape of a Basin? Water Resources Research, 31(4), 1119–1127. https://doi.org/10.1029/94WR03290

Rodríguez‐Iturbe, I., & Valdés, J. B. (1979). The geomorphologic structure of hydrologic response. Water Resources Research, 15(6), 1409–1420. https://doi.org/10.1029/WR015i006p01409

Senatore, A., Micieli, M., Liotti, A., Durighetto, N., Mendicino, G., & Botter, G. (2021). Monitoring and modeling drainage network contraction and dry Down in Mediterranean headwater catchments. Water Resources Research, 57(6), e2020WR028741. https://doi.org/10.1029/2020WR028741

Skaugen, T., & Onof, C. (2014). A rainfall‐runoff model parameterized from GIS and runoff data. Hydrological Processes, 28(15), 4529–4542. https://doi.org/10.1002/hyp.9968

Skaugen, T., Stavang, A. E., Lawrence, D., & Møen, K. M. (2023). Catchment response times – Understanding runoff dynamics from catchment distances and celerities. Hydrological Sciences Journal, 68(8), 1127–1138. https://doi.org/10.1080/02626667.2023.2201449

Sloan, B. P., Mantilla, R., Fonley, M., & Basu, N. B. (2017). Hydrologic impacts of subsurface drainage from the field to watershed scale. Hydrological Processes, 31(17), 3017–3028. https://doi.org/10.1002/hyp.11218

Tarboton, D. G. (1997). A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resources Research, 33(2), 309–319. https://doi.org/10.1029/96WR03137

Tetzlaff, D., Soulsby, C., Bacon, P. J., Youngson, A. F., Gibbins, C., & Malcolm, I. A. (2007). Connectivity between landscapes and Rivers capes—A unifying theme in integrating hydrology and ecology in catchment science? Hydrological Processes, 21(10), 1385–1389. https://doi.org/10.1002/hyp.6701

Tucker, G. E., Catani, F., Rinaldo, A., & Bras, R. L. (2001). Statistical analysis of drainage density from digital terrain data. Geomorphology, 36(3–4), 187–202. https://doi.org/10.1016/S0169‐555X(00)00056‐8

Van Meerveld, H. J. I., Kirchner, J. W., Vis, M. J. P., Assendelft, R. S., & Seibert, J. (2019). Expansion and contraction of the flowing stream network alter hillslope flowpath lengths and the shape of the travel time distribution. Hydrology and Earth System Sciences, 23(11), 4825–4834. https://doi.org/10.5194/hess‐23‐4825‐2019

Volpi, E., Di Lazzaro, M., & Fiori, A. (2013). Analytical modeling of the hydrologic response under moving rainstorms: Storm–catchment interaction and resonance. Journal of Hydrology, 493, 132–139. https://doi.org/10.1016/j.jhydrol.2013.04.025

Ward, J. V. (1989). The four‐dimensional nature of lotic ecosystems. Journal of the North American Benthological Society, 8(1), 2–8. https://doi.org/10.2307/1467397

White, A. B., Kumar, P., Saco, P. M., Rhoads, B. L., & Yen, B. C. (2003). Changes in hydrologic response due to stream network extension via land drainage activities. JAWRA Journal of the American Water Resources Association, 39(6), 1547–1560. https://doi.org/10.1111/j.1752‐1688.2003.tb04438.x

Wiens, J. A. (2002). Riverine landscapes: Taking landscape ecology into the water. Freshwater Biology, 47(4), 501–515. https://doi.org/10.1046/j.1365‐2427.2002.00887.x

Wigington, P. J., Moser, T. J., & Lindeman, D. R. (2005). Stream network expansion: A riparian water quality factor. Hydrological Processes, 19(8), 1715–1721. https://doi.org/10.1002/hyp.5866

Willgoose, G., Bras, R. L., & Rodriguez‐Iturbe, I. (1992). The relationship between catchment and hillslope properties: Implications of a catchment evolution model. Geomorphology, 5(1–2), 21–37. https://doi.org/10.1016/0169‐555X(92)90056‐T