Explaining the climate sensitivity of junction geometry in global river networks

Significance River networks form tree-like branching patterns as they drain landscapes. Research has shown that the angle between joining channels, referred to as the junction angle, is sensitive to climate. We develop a theory based on the principle that river networks self-organize to minimize energy expenditure along their paths and test it using the global dataset of junctions we have extracted. This theory is able to explain much of the observed variation in junction shapes, including their sensitivity to both climate and other factors, such as the drainage area ratio of the two tributaries. The theory does not, however, explain our observation that minor tributaries tend to join large rivers at the outside apex of large-scale bends.

local scales. To account for these discrepancies and to allow for the extraction of fully intact basins, we buffered the edges of 23 each basin tile prior to hydrological processing. For tiles with areas less than 100 km 2 this buffer was fixed at 2 km. For tiles 24 with areas between 100 km 2 and 15625 km 2 the buffer radius (R buf f er ) dependant on the area of the tile (A tile ) as according heads strongly controls the density and topology of the resulting network. While it is possible to apply algorithms to extract 34 channel heads from topographic data (1) they require high pixel resolution and considerable processing time -constraints which 35 made this approach unfeasible in our analysis. Instead we opted to use a simplistic threshold approach, whereby channels are 36 initiated at a threshold drainage area.

37
What constitutes an appropriate drainage area threshold for channel initiation? Using a threshold that is too low may 38 result in the extraction of false junctions where no channels exist. Using a threshold that is too high would result in fewer 39 extracted junctions, sacrificing statistical confidence. The scale at which junction angles are measured is also controlled by the 40 extraction threshold as channel links will be longer when this threshold is higher.

41
In 21 field studies of channel heads (2) the largest recorded initiation threshold area for a channel was 0.8 km 2 with the of the three segments, a line is fitted that minimises the sum of squared perpendicular distances from the data points to the 83 regression line). These bearings are then used to calculate the junction and bending angles using standard vector calculations.

84
The spacing between junctions is determined by the channel extraction threshold: smaller thresholds will have greater network 85 density and thus more closely spaced junctions and shorter channel segments. We use channel segments between junctions 86 rather than segments made up only of channel pixels near the junctions because we are interested in capturing the structure of 87 how channels drain the entire landscape rather than focusing on the hydraulics near the confluence of tributaries. Any junction 88 with a channel segment composed of less than 8 pixels is rejected. This avoids the inclusion of junctions dominated by local 89 hydraulics and, more importantly, provides sufficient degrees of freedom for angular calculations from the linear regression 90 to avoid the production of systematic values. As mentioned previously, we also eliminate junctions where any three of the 91 confluent links have a gradient less than a threshold values to minimise errors from poorly located channels in low relief terrain.

92
The criteria that a junction must include three segments with lengths greater than or equal to 8 pixels means that in large 93 scale images of the channel network, such as an Figure 1, some locations that appear to have a junction are not added to the 94 dataset.
We illustrate examples of such conditions in Figure S2. In the optimal junction model the parameter γ prescribes the relative cost weighting of channel links. When γ = 1 there is 97 no energy-saving advantage for two links to join and form a junction, and hence no optimal junction angle solution. When 98 γ < 1 there is a cost minimising advantage for the two links to join and form a junction. As γ decreases the optimal angle of α 99 becomes wider and the optimal bending angles narrow as the energy-saving advantage of lengthening the resultant channel  Implicit in our calculated optimal junction geometries is the assumption that the energy scaling exponent γ is constant within 106 all the links of a model junction.

107
Our results show that apparent values of γ decrease with increasing drainage area (as is consistent with hydrological evidence 108 for a downstream decrease in c (9-15)). Therefore predictions made by employing the assumption of constant γ cannot be 109 expected to perfectly match reality when the confluent tributaries are disparately sized. We mathematically investigated 110 whether this assumption could be responsible for the discrepancies we observe between theory and reality.

111
Taking the case where the energy scaling exponent γ varies with drainage area, each link (i) with a unique drainage area 112 value (Ai) will have a γ value independent of the other links' γ values such that the weighting of each link in Zamir's (16) 113 optimal junction equations becomes A γ i i : At a junction where two links (i1 and i2) join to form a resultant link (i0), the drainage area is additive according to a 115 simple bifurcation rule: In real rivers, A increases modestly from the upstream end of a link to the downstream end as areas proximal to the link 117 drain into it. For simplicity we disregard these local increases in A and assume the A of any junction link conforms to Eq. 5.

118
The drainage area ratio (AR) of a junction is given by: Using the bifurcation rule (5) and an AR (6) the optimal junction equations (2, 3 and 4) can be simplified to (17): When AR = 1 we find: When AR = 1 optimal junction geometry is sensitive only to γ0, the γ value of the resultant channel. When AR = 1, optimal 122 junction geometry is sensitive to the γ values tributary and resultant channel links.

123
And when AR = 0 we find analytically that optimal junction geometry becomes independent of γ with optimal angles of 124 α = 90 • , β1 = 180 • , and β2 = 90 • . As the larger bending angle is predicted to approach 180 • irrespective of the link γ values, 125 the assumed constancy of γ cannot be producing the discrepancy between observed and predicted junction angles when the 126 drainage area ratio approaches zero. 4. An alternative derivation of optimal junction angle equations using tributary slope ratio 128 The area ratio, AR, is linked to the slope ratio, SR, through the empirical slope-area scaling law (i.e., S = ksA −θ ). It is 129 therefore theoretically possible to derive versions of Roy's (17) optimal junction angle equations (Eqs. 7, 8 and 9) formulated 130 in terms of SR: In the main body of this paper we have presented Roy's (17) original equations using AR. We prefer the use of AR for three 134 reasons. Firstly, arriving at optimal predictions using a version of the optimal junction angle equations incorporating the 135 slope ratio requires the assumption of a reference θ value to transform the bifurcation rule (A0 = A1 + A2) so that S0 can 136 be formulated. Introducing this third variable parameter significantly complicates the process of finding optimal junction 137 configurations. Secondly, drainage area (A) data from real rivers is much less noisy than gradient (S) data. Finally we believe A 138 represents a more fundamental parameter than S in the context of network geometry as S is generally thought to be controlled  When measuring junction geometry it is necessary to approximate each channel link in the channel network as a straight line.

142
In our analysis this is achieved by performing an orthogonal linear regression through the points that compose each channel 143 link. The sinuosity and curvature of real river channels varies, so the 'fit' of the vector line to the real channel link also varies. with gradients exceeding 0.02 (Fig. S9) our dataset are therefore consistent with the known variation of river planform sinuosity with channel slope.

157
Median R 2 values decrease slightly from ∼ 0.78 to ∼ 0.72 as link length increases from 0 to 5 km (Fig. S10). The great 158 majority of links in our dataset fall within this length range. For the rare links longer than 5 km R 2 values increase with length.

159
For links of 20 km or more R 2 values are ∼ 1.

160
Pertinent to our analysis of junction geometry, the R 2 values we observe in our dataset, although noisy, are in general 161 reassuringly high. This demonstrates that, on average, the vectors we have fitted to channel links are reasonable approximations. 162 6. Verification of Roy's optimal geometry equations 163 We have verified Roy's (17) optimal geometry equations with a brute-fore iterative search for the optimal configurations, where 164 'optimal' is defined as the minimum possible value of energy expenditure (E) calculated by: Where wi is the weighting factor of link (i). Calculating E for any junction configuration requires as input the starting [20] Using the rule of cosines the angle γ can be calculated from sides a, b and c: And it follows that: The line O1O2 bisects line L3 and is therefore equal to twice the altitude of O1CO2 with respect to base O1O2. Finally, −2( a 2 sin(β 1 ) )( b 2 sin(β 2 ) ) sin(β1 + β2 + arccos( a 2 +b 2 −c 2 2ab )) ( a 2 sin(β 1 ) ) 2 + ( b 2 sin(β 2 ) ) 2 + 2( a 2 sin(β 1 ) )( b 2 sin(β 2 ) ) cos(β1 + β2 + arccos( a 2 +b 2 −c 2 2ab )) [24] This last step can be repeated to also find L1 and L2.
193  We provide examples of junction selection or absence that can be obscured when looking at the dataset at regional scales. The junctions from this image have not yet had the gradient threshold applied (junctions where any one of the three segments has a gradient < 0.0001 is removed) but junctions with short segments (8 pixels or less) are removed. Panel A shows the same area as in Figure 1 in the main text. Panel B shows a junction (1) where all three segments are greater than 8 pixels but because of flat terrain one segment runs parallel to another and this parallel segment cannot be seen at the scale of panel A. Such junctions are frequently eliminated by the gradient threshold. Panel C appears to have 3 junctions but two of these are removed (junctions 2 and 3) because both have a segment shorter than 8 pixels. The other apparent junction (4) is actually not a true junction due to parallel flow routing -once again this detail is not visible at the scale of Panel A. In panel D both junctions (5 and 6) are also removed due to short link lengths.      S8. Junction angles (α) and bending angles (β1, β2) sorted by climate aridity zone and binned by drainage area ratio (A R ). Black line is the median value, the darker coloured shading represents the 25th to 75th data percentiles, and the lighter coloured shading represents the 5th to 95th data percentiles. The histogram along the right hand side of the plots is a count of the number of junctions in each of the 20 bins. . This plot shows the relationship of the R 2 value from the orthogonal linear regression performed on each channel link (effecively a measure of 'fit') and the average link gradient. Black line is the median value of each data bin, the darker grey shading represents the 25th to 75th data percentiles, and the lighter grey shading represents the 5th to 95th data percentiles. The histogram at the top of the plot is a count of the number of channel links in each bin. . This plot shows the relationship between the R 2 value from the orthogonal linear regression performed on each channel link (effecively a measure of 'fit') and the link length. Black line is the median value, the darker grey shading represents the 25th to 75th data percentiles, and the lighter grey shading represents the 5th to 95th data percentiles. The histogram at the top of the plot is a count of the number of channel links in each bin.