## Abstract

As water erodes a landscape, streams form and channellize the surficial flow. In time, streams become highly ramified networks that can extend over a continent. Here, we combine physical reasoning, mathematical analysis and field observations to understand a basic feature of network growth: the bifurcation of a growing stream. We suggest a deterministic bifurcation rule arising from a relationship between the position of the tip in the network and the local shape of the water table. Next, we show that, when a stream bifurcates, competition between the stream and branches selects a special bifurcation angle *α*=2*π*/5. We confirm this prediction by measuring several thousand bifurcation angles in a kilometre-scale network fed by groundwater. In addition to providing insight into the growth of river networks, this result presents river networks as a physical manifestation of a classical mathematical problem: interface growth in a harmonic field. In the final sections, we combine these results to develop and explore a one-parameter model of network growth. The model predicts the development of logarithmic spirals. We find similar features in the kilometre-scale network.

## 1. Introduction

As groundwater emerges to the surface through springs, it often forms a network of streams [1–5]. Here, we investigate the relationship between the growth and geometry of such *seepage networks*. Our attention to the growth of networks fed by groundwater is motivated by the relative simplicity of groundwater flow. As was first recognized by Henry Darcy, groundwater moves through the subsurface according to the diffusion equation [6]. This broad mathematical context for subsurface hydrology was exploited by scientists and engineers over the next century to find exact solutions for groundwater flow in complex geometries. Pre-eminent among these researchers was the applied mathematician Pelageya Polubarinova-Kochina, who demonstrated the power of functions of a complex variable in understanding the flow of groundwater [7]. The role of subsurface flow in shaping landscapes was described most intuitively by Thomas Dunne, who recognized that the flow of groundwater shaping a network of streams is itself shaped by the geometry of the network [1]. Here, we develop a quantitative understanding of this interplay between geometry and growth by relating seepage networks to a classical mathematical problem: interface growth in a harmonic field [8–13].

In this paper, we focus on a particular seepage network on the Florida Panhandle [14–16]. An elevation map of a section of this network is shown in figure 1*a*. The growth of this network is driven by the flow of groundwater through a layer of nearly homogeneous unconsolidated sand [14]. Originally, groundwater emerged in springs at the base of a steep bluff. However, over the last several hundred thousand years, water flowing through the springs has eroded the landscape, causing the springs to migrate forward and new springs to form, eventually developing a highly ramified network [15]. A modern spring is shown along with the junction of two streams in the network in figure 1*b*,*c*, respectively.

This paper is about the dynamics of stream growth. We begin by considering the different length scales and time scales over which groundwater shapes a growing network of springs. Because the growth of streams is slow compared with the flow of groundwater, we separate the growth of streams into two problems. In §3, we find the relationship between the shape of a stationary stream and the flow of groundwater into it. We then describe the growth of the valley in response to this flow in §4. In this section, we derive our main result: the dynamics of spring bifurcation are recorded in the bifurcation angle *α*=2*π*/5. In §5, we confirm this prediction in the Florida network. Finally, in §§6 and 7, we combine the results of the previous sections to grow seepage networks numerically. We find that as seepage networks grow, successive bifurcations produce logarithmic spirals in the network. Similar features appear in the Florida network.

Our theoretical and empirical analyses of the bifurcation angle were recently reported in the study of Devauchelle *et al*. [17]. The following five sections of this paper provide much additional detail that allows us to embed the bifurcation problem into more general issues of network growth examined in subsequent sections. Because our use of complex analysis, although classical in physics, is new in geomorphology, we also hope that this more detailed presentation renders our earlier analysis workable to a broader readership.

## 2. The scales of groundwater flow

Groundwater flows along subsurface pressure gradients [6,7]. When the principal pressure gradients are hydrostatic, the Dupuit approximation relates the height *h* of the water table relative to an impervious substratum to the flux *q*=−*kh*∇*h*, where *k* is the hydraulic conductivity [18,19]. From conservation of mass, the height of the water table changes in response to the divergence of groundwater flux and rainfall; thus
2.1where *t* is time and *P* is the precipitation rate.

The typical length scale of groundwater flow, which we will call drainage length, is , where is the area drained by the network and *l* is the total length of streams. This length is the inverse of the so-called *drainage density*. In the Florida network, . Because the width of streams , streams are effectively one dimension.

There are three relevant time scales for the flow of groundwater. First, the typical time *τ*_{p} between perturbations to the water table due to storms is approximately 10 days. Once rain is stored as groundwater, it returns to the surface on the second time scale *τ*_{f}. The time scale for the flow of groundwater out into the network is days, where we have taken *k*∼10^{−7} m s^{−1} [16]. Because *τ*_{p}≪*τ*_{f}, the water table does not fully drain between storms. Rather, the decay of individual storms is so long that only the average rainfall rate influences the flow of groundwater. Finally, the flow of water through streams erodes the landscape, causing springs to grow forward at a speed *v*∼1 mm y^{−1} [15]. The typical time scale of the network growth is therefore year. Thus, the landscape evolves quasi-statically with respect to the flow of groundwater, which depends only on the average rainfall rate. This separation of time scales allows the growth of streams fed by groundwater to be divided into two problems. First, we ask how the flow of groundwater into a stationary spring depends on the position of the spring in the network. We then investigate the dynamics of spring growth in response to groundwater flow.

## 3. Groundwater flow into a static network

On time scales short compared to the time scale of network growth, the curvature of the water table balances rainfall. Similarly, on time scales long compared to the frequency of storms, the water table feels only the average effect of many independent storms. In this limit, equation (2.1) becomes
3.1where *P*_{0} is the average rainfall rate.

To predict how groundwater flows, equation (3.1) must be combined with appropriate boundary conditions on the streams. Because the sand is extremely porous, any water on the surface will immediately seep into the pore space. Consequently, water flowing on the surface is only possible if the height of the topography is lower than the height of water table. A spring is the first point where the topography falls below the water table. Streams form where the shear imposed by the flowing water is sufficient to incise the topography below the water table. We consider only the case in which water flows in a network of streams as opposed to forming swamps and pools. We note that the dynamics of sediment transport often lead to a channellized surficial flow [20]. The water table therefore intersects the network at the elevation of the streams. In the Florida network, the median slopes of stream beds is *S*∼10^{−2}≪1. Neglecting the small variation in height, groundwater flows into the network at a constant elevation, which we define as *h*=0. This approximation, in combination with equation (3.1), has been found to correctly predict the flow of groundwater into the Florida network [21].

A useful non-dimensionalization of equation (3.1) can be found by solving for the shape of a water table in one dimension. Figure 2 depicts two parallel streams separated by a distance . Integrating equation (3.1) and requiring that *h*=0 on the streams,
3.2The typical relief of the water table is . Defining the dimensionless potential , equation (3.1) becomes
3.3where ℓ is a length scale used to render the Laplace operator dimensionless. By expressing the balance of water table curvature with rainfall in dimensionless form, equation (3.3) identifies the inverse drainage density as the typical scale over which rainfall shapes the water table. If one considers the shape of the water table immediately around a single stream—as illustrated in figure 3—the rain falling within a distance ℓ of the channels represents a negligible fraction of the total groundwater. In this limit (i.e. ), precipitation can be neglected.

### (a) The complex potential

In this section, we consider the flow of groundwater into a stream. Close to the stream, equation (3.3) can be approximated as
3.4As before, the boundary condition on the stream is *ϕ*=0. However, the boundary condition far from the stream must be chosen to match the full solution of equation (3.3). Thus, without solving for the groundwater flow throughout the network, equation (3.4) is sufficient to determine the functional form of the water table near a stream. However, parameters such as the total water discharge are determined by the position of the stream in the network.

To solve equation (3.4) around streams, we move to the complex plane: rather than referring to points using Cartesian coordinates (*x*,*y*), points are represented as complex numbers *z*=*x*+i*y*, where i^{2}=−1. The advantage of this transformation is that both the real and imaginary parts of a function *Φ*(*z*) are solutions of equation (3.4) provided *Φ*(*z*) can be expressed as a power series [22]. In the language of complex analysis, *Φ*(*z*) must be *analytic*. Thus, to find the shape of the water table around a stream, one needs only to find a function *Φ*(*z*)=*ϕ*(*z*)+i*ψ*(*z*), the real part of which vanishes along the stream.

The *complex potential* *Φ*(*z*) gives the shape of the water table without solving any differential equations. As an example, we consider the shape of the water table near a straight section of stream in a network as depicted in figure 3. Expanding the complex potential
3.5where *a*_{1}, *a*_{2} and *a*_{3} are coefficients chosen to match two boundary condition. First, the elevation of the water table must flow into the stream at the elevation of the stream, thus *ϕ*=0 on the real axis (i.e. *z*=*x*+0i). Furthermore, these coefficients must force *Φ*_{0} to match the solution of the Poisson equation (3.3) far from the stream. Applying the first boundary condition, *ϕ*(*x*)=0 only if each of the coefficients is real.

### (b) Flow into a spring

Having found the flow of groundwater into the side of a stream, we now turn our attention to a more complicated geometry: the flow of groundwater into a spring. Sufficiently close to a spring, a stream can be approximated as straight. We choose a coordinate system in which the stream lies along the negative imaginary axis and the spring is at the point *z*=0.

To solve for the complex potential around a spring, we make use of a classical observation about functions of a complex variable: the composition of analytic functions is also analytic [22]. Thus, if two functions *Φ*(*z*) and *g*(*z*) are solutions to the Laplace equation, then so is *Φ*(*g*(*z*)). The power of this simple observation is in the interpretation of *g* as a change of coordinates. The function *g*(*z*) maps a point *z* to a new point *w*=*g*(*z*). Figure 4 shows that the coordinate change maps the space around the negative imaginary axis to the upper half plane. This change of coordinates, therefore, reduces groundwater flow around a spring to a simple geometry in which the shape of the water table is easily expressed as equation (3.5). In the language of complex analysis *g*(*z*) is a *conformal map* [22]. Combining this coordinate change with equation (3.5), the complex potential around a spring is
3.6where *a*_{1}, *a*_{2} and *a*_{3} are real coefficients determined by the boundary conditions far from the stream.

Each coefficient in equation (3.6) has a physical interpretation related to the position of the spring in the network. At a small distance *ϵ* from the spring, the slope of the water table is determined by the lowest order term in *Φ*_{1}. The discharge in the stream is . Thus, *a*_{1} is proportional to the flux of water into the spring. Physically, *ϵ*∼*w*∼30 cm is the typical width of the stream near the spring. Because the flux into the stream is assumed to be positive, it follows that *a*_{1}>0. Figure 5*b* shows the flow field for *Φ*_{1}=(i*z*)^{1/2}. If an obstruction (e.g. another stream) is to one side of the channel (figure 5*c*), the flow into the spring is asymmetric. As illustrated by figure 5*d*, this behaviour is captured by the coefficient of the odd term i*z*. The third term in the expansion describes the convergence of groundwater into the spring (figure 5*e*–*h*). As shown in figure 6, when a stream grows to a point in the network where *a*_{3} passes a critical value *a*_{c}, the shape of the water table around the tip undergoes a bifurcation in which water begins to flow from two maxima rather than one. When *a*_{3}>*a*_{c} the largest flux into the spring is directly above the stream. When *a*_{3}<*a*_{c}, the largest fluxes come from two points on either side of the stream. A symmetric water table (*a*_{2}=0) bifurcates when the real part of *Φ*_{1} develops a local minimum directly in front of the stream. To find this point, we first move to polar coordinates *z*=*r* e^{i(θ+π/2)}. In this coordinate system, the real part of *Φ*_{1} is
3.7The bifurcation of the water table occurs when
3.8Recalling that the *a*_{1}>0, the critical value of *a*_{3} is
3.9Thus, whenever *a*_{3}<0, the water table develops a bifurcation. When ∥*a*_{3}∥≪1, the bifurcation is far from the tip. As *a*_{3} becomes more negative, the bifurcation moves closer to the spring. There is an analogous bifurcation for asymmetric bifurcations (i.e. *a*_{2}≠0).

Higher order terms in equation (3.6) provide progressively smaller corrections to the flow. Notably, in this special geometry, one can take rainfall into account by adding the field to *Φ*_{0}(*z*). Because this field is , precipitation does not affect the first three terms in the expansion of the water table.

In equation (3.9), *a*_{c} depends explicitly on the choice of *r*. Consequently, the proposed dynamics that lead to a spring to bifurcate introduce a new length scale. To make clear the importance of this length scale, we abandon notion of a critical moment of the water table (i.e. *a*_{c}) in favour to the equivalent critical length
3.10From equation (3.9), springs bifurcate if −*a*_{1}/9*a*_{3}<*r*_{c}. Practical problems make a direct measurement of this length-scale difficult. In particular, because the time scale for the growth of the network is measured in millennia, it is not possible to observe a bifurcation in nature, and consequently it is difficult to test the functional form of equation (3.9).

There are two strong candidates for this new length scale. According to the first, *r*_{c} is determined by the geometry of the network. For example, *r*_{c} could be of order the width of a valley (approx. 100 m). According to this hypothesis, when the water table develops local maxima within a single valley head, each maximum forms a spring which erodes the valley walls. As each spring grows, the valley bifurcates. Alternatively, *r*_{c} could be determined by the dynamics of network growth. According to this hypothesis, random perturbations within the valley head typically lead to the formation of several springs. When *a*_{3}>0, each of these incipient springs compete for water from a single direction. As small differences in position lead to large differences in flux, one spring outcompetes the rest, leading to an intact valley head. This behaviour is analogous to the mathematical system studied by Carlson & Makarov [10]. When *a*_{3} is ‘substantially’ negative, incipient springs growing towards different maxima compete less and less; causing the valley to bifurcate. According to this interpretation, *r*_{c} is the typical distance incipient springs grow before one outcompetes the rest. Notably, in this second case, the stability of the incipient springs could depend on higher order terms in the expansion of the water table. We shall defer judgement on this point and simply interpret this model as a means to generate intuition for the natural system.

A final determination of which, if either, of these two possibilities is correct must rely on an understanding of how these dynamics are recorded in the current geometry of the network (e.g. the distance a channel grows before bifurcating). Although such an understanding is introduced in §7, a test will be left for future work. In the following three sections, we therefore only assume that the dynamics of network growth occasionally cause springs to bifurcate; however, we do not rely on any specific assumptions regarding the cause.

## 4. Dynamics of network growth

Because springs grow in response to the flow of groundwater, physical reasoning suggests that springs grow in the direction from which groundwater flows. This assumption is the central hypothesis of this paper. Proceeding from this hypothesis, in this section we derive our main result: a prediction of the bifurcation angle of streams. We first find the paths—or *streamlines*—along which groundwater flows into a single spring. As we shall see, the streamline flowing into a spring offers a useful characterization of the competition for groundwater between the spring and the rest of the network. From this characterization, we derive the bifurcation angle that balances competition between the branches and competition with the rest of the network.

### (a) Growth of a stream

Thus far in the analysis, the power of complex functions has been exploited chiefly in the interpretation of *ϕ*(*z*) as the shape of the water table. We now turn our attention to the imaginary part of *Φ*(*z*). The physical interpretation of *ψ*(*z*) requires a classical result of complex analysis. It follows from the Cauchy–Riemann equations [22] that the curve *ψ*(*z*)=const. is orthogonal to the elevation contours of the water table. Thus, the imaginary part of *Φ* is constant along the paths of steepest descent down the water table. Therefore, *ψ*(*z*) gives the path along which groundwater flows [7] and, consequently, the direction a spring grows.

As an example, we return to the flow of groundwater into an infinite straight stream. Neglecting cubic terms in equation (3.5), the complex potential can be expressed as
4.1
4.2The streamlines are curves on which *ψ*(*u*+i*v*) is constant. In particular, as shown in figure 7, the streamline given by *ψ*(*w*_{0})=0 defines the curve
4.3along which groundwater flows into the point *w*_{0}(0)=0. Note that this parametrization of the streamline requires that the parameter *u* has the same sign as the coefficient *a*_{2}.

Despite the simplicity of this example, equation (4.3) is useful in determining the flow of groundwater into streams with more complex geometries. To find the paths along which water flows into a stream with any shape, we only need to identify a change of coordinates that maps the infinite straight stream into the more complex geometry. In the following derivations, we use conformal mappings to relate equation (4.3) to the flow of groundwater into a spring and the resulting growth of the stream.

To find the shape of the streamline flowing into a spring, we must identify a function that maps the real axis onto a straight stream ending in a spring. This mapping is the inverse of the change of coordinates used in the derivation of equation (3.6): *z*=−i*w*^{2}. Because the point *w*=0 is mapped onto the spring at *z*=0, the streamline that flows into the spring is
4.4This streamline is shown in figure 7.

Proceeding from the hypothesis that streams grow in the direction from which groundwater flows, equation (4.4) relates the shape of the water table, through the parameters *a*_{1} and *a*_{2}, to the direction in which the stream grows. As a spring grows forward, its position in the network changes. These changes are reflected in the new values of the three coefficients *a*_{1}, *a*_{2} and *a*_{3} that determine the flow of groundwater into the tip and, therefore, the growth of the tip. In the following section, we investigate how the geometry of a bifurcated stream affects the flow of groundwater into each of the nascent springs, the growth of the springs and the resulting geometry of networks cut by sub-surface flow.

### (b) Dynamics of stream bifurcation

When a stream bifurcates, the flow of groundwater into the newly formed springs is influenced chiefly by competition between the springs, rather than by the distant network. It follows from equation (3.6) that if the length of each branch is *δ*≪(*a*_{1}/*a*_{2})^{2}, the flow is dominated by the lowest order term in the expansion of the complex potential. Recalling that the first coefficient *a*_{1} only determines the magnitude of the flux of water entering the stream, the position of the stream in the network does not affect the direction from which water enters the newly formed springs. The direction in which nascent springs grow is therefore only influenced by the shape of the bifurcation. Thus, although the bifurcation of the spring is related to changes in *a*_{3}, this coefficient only weakly affects the growth of the springs immediately after the bifurcation.

As shown in figure 8, when the bifurcation angle *α* is small, the springs compete with one another for groundwater, causing the nascent streams to grow apart. When *α* is large, the springs compete strongly with the main stream, causing the streams to grow together. In this section, we identify the special angle at which the springs of a newly bifurcated tip grow forward without curving.

Thus far in the analysis, we have assumed that, sufficiently close to a spring, a stream can be approximated by a straight line. However, when a stream bifurcates, this approximation breaks down. Consequently, equation (3.6) is no longer an adequate description of the flow of groundwater into the spring. To predict the flow of groundwater into a bifurcated stream, we must identify a new change of coordinates in which the shape of the surrounding water table is simply expressed. As shown in figure 9, the function
4.5maps the real axis onto a bifurcated stream. In this new coordinate system, the two springs lie on the real axis at the points *w*=±1. As shown in figure 9, equation (4.5) maps these points onto two short streams which meet at an angle *α*. Thus, the nascent streams lie on the lines
4.6where 0<*r*<1. Notably, the only length scale needed to describe this bifurcation is the length of the nascent streams. Consequently, if the growing streams maintain this initial shape, the length of the streams can be rescaled by time to reveal an unchanging configuration. Such a configuration is called *self-similar*.

To find the direction that the springs grow, we must find the shape of the streamlines flowing into the tips. Because the shapes of streamlines are simply expressed when flowing into an infinite straight stream, we use *g*_{0}(*w*) to relate the flow in this simple geometry to the flow into a bifurcated spring. Shortly after a tip bifurcates, the flow of water into the nascent branch is only affected by the lowest order term in the expansion of the complex potential. Neglecting quadratic terms in equation (3.6),
4.7where the coefficient *a*_{1} has been incorporated into the definition of *Φ*_{2}. As shown in figure 9, the streamlines flowing into the springs at *w*=±1 are mapped from the vertical lines *w*_{±}(*v*)=±1+i*v*. The streamline entering the spring at *w*=+1 follows the curve
4.8in the physical plane, where
4.9The corresponding streamline *z*_{−}(*v*)=*g*(*w*_{−}(*v*)) is obtained by symmetry. Given the shape of the streamline flowing into the spring for an arbitrary value of *α*, we now identify the special bifurcation angle for which the two nascent branches enter the streams without curving. The curvature *κ*_{α}(*v*) is given by (see appendix A)
4.10where ℑ(⋅) signifies the imaginary part. Expanding equation (4.10) around the tip (*v*=0) yields
4.11As the streamline approaches the tip (), the curvature vanishes for
4.12and consequently streams grow without curving. At all other angles between 0 and 2*π*, the curvature is infinite, diverging like 1/*v*. We therefore identify *α*=2*π*/5 with self-similar growth.

The growth of needle-like branches in harmonic fields has long been recognized as a physical system of singular elegance. Historically, this subject has been motivated both by a combination of purely mathematical questions [10,23–25] and problems of diffusion-limited aggregation [26–28] and its application to turbulence and phase transitions [12,29,30]. The prediction of the self-similar bifurcation angle *α*=2*π*/5 is a special case of several previous analyses. Selander [31] investigated the stability of one-dimensional interfaces growing in Laplacian fields and showed that two branches growing off the tip of a branch-cut are stable if their opening angle is 2*π*/5. Carleson & Makarov [10] extended these results to more general settings. Hastings [28] derived the same angle in an analysis of the dielectric breakdown model at the upper critical point, assuming, as we have, that the shape of the bifurcation does not change during the initial steps of growth. Finally, Gubiec & Szymczak [32], following a lucid discussion of needle growth in a harmonic field, predict that the angle between two needles growing from the side of an infinite straight boundary asymptotically approaches *π*/5. To adapt this result to the bifurcation of a spring, one need only map the infinite boundary with the coordinate change *z*=−i*w*^{2}. Doing so doubles the bifurcation angle to 2*π*/5.

## 5. Bifurcations in the Florida network

To test the hypothesis that springs grow in the direction from which groundwater flows, we compare the geometry of the Florida network to model predictions. The comparison of the observed bifurcation angle to the prediction *α*=2*π*/5 is a particularly good test of the model for two reasons. First, as shown in figure 10*a*, myriad springs in the Florida network allow for many independent measurements of the bifurcation angle from which to estimate its characteristic value. Second, the comparison between theory and observation requires no fit parameters.

To measure the bifurcation angles of streams in the Florida network, we must first extract the network of streams from the topography. Streams stand out in the landscape as sharply incised cuts through the otherwise smooth topography. Consequently, to extract the streams automatically from a high resolution (1.3 m horizontal, 0.05 m vertical) elevation map (figure 10), we identify all points where the elevation contours of the topography are sharply curved [16]. In practice, a point is considered to be in a stream if the contour curvature *κ*>10^{−0.9} m^{−1}. Because this algorithm also identifies sharply incised rills on the steep valley walls well above the known height of the stream, an additional threshold is placed on the slope *S* of the topography, *S*<10^{0.1}. Finally, the one-dimensional network passing through these points is extracted by skeletonizing [33] a binary image of stream locations. Using the location of the channels we can identify the different streams and index them by their Horton–Strahler order [34], as shown in figure 10*b*. Having the ordered stream we approximate the growth direction for each channel by a straight line using the reduced major axis [35]. This method takes into account that both *x* and *y* coordinates have similar error bounds. Furthermore, the crude approximation by a straight line averages out random fluctuations and meandering of the stream. We then define *α* as the angle between the two upstream segments of intersecting stream directions.

Figure 11 shows the histogram of bifurcation angles estimated from 4966 bifurcated streams (data available in the electronic supplementary material). The mean bifurcation angle 〈*α*_{m}〉=71.9^{°}±0.8^{°} (95% CI) is not significantly different from the predicted value of 2*π*/5=72^{°}. The median angle of the distribution is *M*_{α}=72.07^{°}. Additionally, we observe that the distribution of bifurcation angles is approximately normal. Indeed, a two sided goodness-of-fit test (Lillifors test) [36] cannot reject this possibility (*p*=0.28).

The excellent agreement between theory and observation leads us to conclude that the mean of the distribution is determined by the self-similar growth in a harmonic field. Notably, the growth law presented in this paper defines only the direction that the channels grow and makes no assumptions about the functional form of how fast the channels grow; consequently, this bifurcation angle is typical of a class of growth laws (i.e. those that grow in the direction from which groundwater flows). Conversely, the variance of the distribution is determined by how perturbations to the bifurcation angle and channel lengths are amplified or diminished by the dynamics of channel growth. Thus, in addition to explaining the mean of the distribution of bifurcation angles, this result identifies the variance as a likely record of the growth law specific to a landscape.

## 6. The motion of a spring

In §4*b*, we found the angle between bifurcated tips by considering the directions the springs grow from an initial configuration. Notably absent from this derivation was a mathematical description of how the position of the springs changes in time. Here, we briefly outline two formalisms by which one can grow streams.

Equation (4.4) gives the path along which groundwater flows into the spring. In the simplest manifestation of growth, a stream grows by extending a short distance along this curve. As illustrated in figure 12, when the spring grows forward along the streamline a small distance d*s* it turns by an angle *β*. The resulting change in the position of the spring is therefore
6.1Expanding equation (4.4) near the tip (*u*=0) and re-expressing *u* in terms of d*s* relates the angle *β* to the shape of the water table through the equation
6.2These two equations provide the basis by which one can grow streams numerically (§7).

A notable complication of this equation is that, as *ds* tends to zero, *β* vanishes, apparently causing springs to grow straight. This peculiarity is clarified by recognizing that the curvature of the streamline *κ*∼*β*/d*s*∼d*s*^{−1/2} diverges, causing the streamline to curve sharply near the tip. Thus, in the limit of vanishing d*s*, both the direction and curvature of the streamline must be taken into account when determining the growth of the spring. Here, we briefly outline an alternative formulation of tip growth that deals more naturally with this complication.

It is useful to consider the time evolution of a bifurcation conceptually. In the framework developed in §4, to find the direction the springs grow we first identify a change of coordinates *g*_{0}(*w*) that maps an infinite straight stream onto the bifurcated stream. Using this mapping, it is straightforward to find the path along which groundwater flows into each spring and thus to grow the streams. The result is a new configuration of springs that is only slightly different from the initial geometry. To grow the springs further, we must identify a new mapping *g*_{δt}(*w*), which maps the infinite straight stream into the new geometry. The key insight is that the growth of the bifurcation is described entirely by the changes in the mapping. Thus, rather than growing the streams directly and being forced to find ever more complicated mappings, it is simpler to evolve the mapping. In this formulation, one uses the growth law for the boundary to derive a differential equation—called a *Loewner equation*—for the evolution of the mapping as the boundary grows [11,29]. There is large literature devoted to the application of the Loewner equation to understanding the growth of needle-like curves in response to both deterministic [10,32] and stochastic [25,29,30,37,38] forcing.

The Loewner equation can be used to describe the motion of an isolated spring growing in response to the complex potential *Φ*_{1}(*z*) given in equation (3.6). As described by Gubiec & Szymczak [32], the Loewner equation for a function *h*(*z*,*t*) that maps the potential around a needle onto the linear potential *Φ*_{ℓ}(*w*)=*w* is
6.3where *v* is the speed the of tip growth. Applying this result to the growth of a stream in a Poisson field is straightforward if one can identify a change of coordinates *w*=*h*(*z*,0) that maps the initial complex potential *Φ*_{1}(*z*) around the stream onto *Φ*_{ℓ}(*w*). Given this mapping, equation (6.3) describes how this initial mapping changes in time, thereby describing the growth of the needles. Finding such a mapping is trivial. Rather than using the function as before, which straightens the stream while contorting the streamlines, we instead use the complex potential itself to define a new coordinate system: *w*=*Φ*_{1}(*z*). This transformation exploits the observation that the orthogonal contours of *ϕ* and *ψ* provide a natural coordinate system in which the field *Φ*_{1} is simply represented.

## 7. Network growth in a Poisson field

Having developed a description of the growth of springs in response to a harmonic field, we now investigate the growth of networks.

The first step in growing a network of streams is to understand how groundwater is partitioned throughout the network. To do so, we solve the Poisson equation (3.1) for the shape of the water table *ϕ*. As shown in figure 13, the network is grown in a square box of area . On three sides of the domain, the groundwater flux vanishes, while on the fourth side the water table reaches a constant elevation *ϕ*=0. The streams are one-dimensional boundaries on which *ϕ*=0. In these simulations, the network grows from a straight channel of length 0.01 on the lower boundary. We solve for the resulting Poisson field using the finite-element method as implemented by the software FreeFem++ [39].

Given the shape of the water table around each spring, the next step is to grow the springs. According to equation (3.9), when *a*_{3}>*a*_{c}, the groundwater flows into the spring from a single direction, causing the spring to grow forward. Assuming the springs grow forward a distance d*s*, equations (6.1) and (6.2) relate the coefficients *a*_{1} and *a*_{2} to the new position of the new spring. Following past work [15,21,40], we assume that the velocity *v* at which a tip grows is proportional to the discharge; thus *v*=*a*_{1}. If *a*_{3}<*a*_{c}, the groundwater no longer flows into the spring from a single direction, causing the spring to bifurcate. When the spring bifurcates, the angle between the nascent springs is *α*=2*π*/5. A network grown in a Poisson field in which *a*_{c}=0 is shown in figure 15. This network grew from a straight incipient channel of length 0.01 on the lower boundary. Physically, this initial condition corresponds to the growth of a network from a randomly nucleated incipient channel.

To estimate the values of *a*_{1}, *a*_{2} and *a*_{3} from the solution *ϕ* of the Poisson equation (3.1), we exploit a useful property of the shape of the water table around a spring. The functions in the expansion of the complex potential in equation (3.6) are an orthogonal basis for *Φ*_{1}(*z*). In this space, the inner product is
7.1where 𝒞 is a circle of radius *ϵ* centred at the spring and *δ*_{mn} is the Kronecker delta. Given this inner product, any coefficient in the expansion can be found by projecting *ϕ* onto the real part of this basis. Thus, the *n*th coefficient in the expansion of the water table is
7.2In practice, *ϵ* is chosen to be small compared with the distance between neighbouring springs, but large enough so that *ϕ* is well resolved by the numerics within 𝒞.

As shown in figure 14, the qualitative shape of the network depends on the choice of *r*_{c}, the length scale that, as described in equation (3.10), determines when a spring bifurcates. When *r*_{c} is small compared with the size of the drainage basin (i.e. *r*_{c}≪1), the bifurcation of springs is not influenced by the distant boundaries. In this limit, *r*_{c} is the typical distance a tip grows before bifurcating. However, when *r*_{c}≫1, the only length scale that influences the growth of the network is the size of the drainage basin. Notably, as the network grows to fill the drainage basin, even this length scale is lost. As a result, the long-time trajectories of the springs are self-similar.

As shown in figure 15, the curve defined by turning right at each bifurcation traces out a spiral. Qualitatively similar structures can be found in the Florida network, one of which is illustrated in figure 16. To characterize this shape, we first approximate each channel segment in this spiral as straight, the length of the *m*th segment being *λ*_{m}. As each bifurcation angle is *α*=2*π*/5, it follows that the position of the *m*th vertex, as represented in the complex plane, is
7.3In the limit *r*_{c}≫1, there are no dynamically relevant lengths; consequently *λ*_{m}/*λ*_{m−1}=*Λ*, for some positive constant *Λ*<1. The dimensionless *screening length* *Λ* is determined by the Poisson equation and the boundary conditions on the network. The positions of vertices are elements of the geometric series
7.4As illustrated in figure 15, the resulting spiral grows towards the point
7.5Moreover, the angle *χ*_{m} between the tangent of the *m*th channel segment and the ray connecting *ζ*_{m} and is
7.6where ℑ(⋅) denotes the imaginary part. Because *χ*_{m} is independent of *m*, the vertices asymptotically lay on a logarithmic spiral [41]. Alternatively, as shown in figure 16, this spiral can be constructed from a rhombic gnomon. Because the bifurcation angle is *α*=2*π*/5, each rhombus is composed of two isosceles triangles with an aspect ratio of the golden ratio . So-called *golden triangles* were of some mystical significance to the Pythagoreans [42]. Thus we find that, in the limit *r*_{c}≫1, streams cut by subsurface flow share a common form with mollusc shells [42] and the flight of hawks [43].

## 8. Discussion and conclusion

In this paper, we have shown that when a spring bifurcates to form two streams, the bifurcation angle between the nascent streams is *α*=2*π*/5. Notably, the derivation of the bifurcation angle only requires that springs grow in the direction from which groundwater flows and that the initial growth of streams be self-similar. Both of these assumptions can be motivated by a sense of parsimony [44]. First, the streamline flowing into the spring provides the only natural preferred direction for growth. Similarly, if the nascent streams were to break self-similarity by curving as they grow, then the initial curvature of the streams introduces a new and unmotivated length scale into the dynamics of stream growth. Alternatively, the second assumption can be relaxed entirely by assuming a particular relationship between the flux of water entering the stream and the rate at which it grows [10]. Thus, the dynamics of spring bifurcation do not depend explicitly on any processes taking place within the stream. The strong agreement between theory and observation leads us to conclude that this hypothesis is accurate.

By way of contrast, we briefly consider the vascular network of animals. According to Murray’s law, veins connect in a way that minimizes viscous dissipation [45,46]. By considering the work required to push a laminar flow through cylindrical veins, Murray found the symmetric bifurcation angle *γ*=75^{°} [47]. In this system, the bifurcation angle results from processes taking place within the network. Finally, although *γ* is not far from the predicted bifurcation angle of streams *α*=72^{°}, this value of *γ* assumes laminar flow through cylinders and is therefore not directly applicable to the confluence of turbulent streams (Re∼10^{4}). When similar reasoning is applied to rivers, the angle is a function of slope *s* and discharge *Q* [48]. In the Florida streams, theory and observation provide the relation *Q*∝*s*^{−2} [16]. Minimizing dissipation within streams would then predict a bifurcation angle of 90^{°} for two streams of equal discharge and slope [48], in disagreement with our measurements.

Although our derivation of the bifurcation angle is only directly applicable to streams fed by groundwater, the underlying reasoning is general. In any network of thin channels that grows by draining an area, there is some path leading into the channel tips. It is reasonable to suspect that channels may generally grow along these paths. To find the bifurcation angle in a more general system, one need only understand the interplay between network geometry and the shape of the path leading into a tip. This observation suggests that the bifurcations of many other networks, such as rivers, ice streams, plant roots and fungal hyphae [49] differ from seepage networks only in small details.

## Funding statement

This work was supported by Department of Energy Grant FG02-99ER15004. H.J.S. was additionally supported by the Swiss National Science Foundation.

## Acknowledgements

We would like to thank The Nature Conservancy for access to the Apalachicola Bluffs and Ravines Preserve, and K. Flournoy and D. Printiss for guidance on the Preserve. We thank Y. Couder, L. Kadanoff, M. Mineev-Weinstein, P. Szymczak, and G. Vasconcelas for helpful discussions and L. Brooks, K. Kebart and the Northwest Florida Water Management District for providing high-resolution topographic maps. We thank one of the reviewers for useful suggestions.

## Appendix A

Let *f*(*v*) be a curve in the complex plane parametrized by the real variable *v*. The curvature *κ* is defined by [50]
A1where *θ* is the tangential angle with respect to the real axis. The derivative of *f* can be expressed as
A2Equivalently,
A3where the overbar signifies the complex conjugate. Taking the logarithm, and then the derivative with respect to *v*, leads to
A4Solving for *θ*′, one obtains
A5where ℑ(⋅) denotes the imaginary part. Inserting equation (A.5) into (A.1), we find
A6

## Footnotes

One contribution of 13 to a Theme Issue ‘Pattern formation in the geosciences’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.