## Abstract

Streams shape landscapes through headward growth and lateral migration. When these streams are primarily fed by groundwater, recent work suggests that their tips advance to maximize the symmetry of the local Laplacian field associated with groundwater flow. We explore the extent to which such forcing is responsible for the lateral migration of streams by studying two features of groundwater-fed streams in Bristol, Florida: their confluence angle near junctions and their curvature. First, we find that, while streams asymptotically form a 72° angle near their tips, they simultaneously exhibit a wide 120° confluence angle within approximately 10 m of their junctions. We show that this wide angle maximizes the symmetry of the groundwater field near the junction. Second, we argue that streams migrate laterally within valleys and present a new spectral analysis method to relate planform curvature to the surrounding groundwater field. Our results suggest that streams migrate laterally in response to fluxes from the surrounding groundwater table, providing evidence of a new mechanism that complements Laplacian growth at their tips.

## 1. Introduction

River networks form captivating geometries, but their present-day position often obscures the rich history of their growth and evolution [1]. Fortunately, the processes that shape this history, such as the migration of river bends [2,3], stream capture [4,5], and confluence movement and development [6,7], often leave behind morphological evidence. However, such processes have not been well studied for groundwater-fed streams. In such landscapes, re-emerging groundwater erodes surface grains by sapping and undermining hill slopes [8,9], forming shallow streams and steep ravines around them. These streams advance in a planimetric direction that has been shown to be largely determined by the shape of the surrounding groundwater field rather than dynamics within the stream [10,11]. In particular, recent work suggests that this advance maximizes the symmetry of the local groundwater field [12,13]. For a confluence, this predicts an angle of 2*π*/5=72° between confluent streams, and these observations have been observed in the planimetry of a groundwater-fed channel network in Bristol, Florida [10].

However, there is no indication that the orientation of the streams is set once the valley has been initially incised. We therefore consider a third phase of growth: lateral stream rearrangement. We suggest that this rearrangement is also influenced by the surrounding groundwater field. In particular, we hypothesize that lateral migration is driven by variable flux to either side of a stream and the subsequent filling-in of old stream paths through soil movement. Here, we study two morphological features of a groundwater-fed river network to constrain the extent to which this groundwater field can describe network-wide morphology: confluences and planform curvature.

This rearrangement may be of wider interest, both observationally and theoretically. First, recent work has demonstrated that an average branching angle of 72° appears in humid climates [14], suggesting that the influence of groundwater on river planform morphology may be widespread wherever shallow groundwater is present. Second, as we will shortly explain, the height squared of the groundwater table solves Poisson’s equation and, close to the streams, Laplace’s equation. Because these streams can be approximated as effectively one-dimensional, their advance can be understood as an example of thin-finger Laplacian growth [12,15–17]. Our results therefore present a novel mechanism for and example of lateral migration of a thin-finger boundary in response to a Laplacian field.

In what follows, we present evidence of this rearrangement from a network of streams in Bristol, Florida, and discuss how it may arise. We first review a theoretical framework that gives the shape of the groundwater table and the 72° branching angle. We also present an averaging scheme that reveals a striking 120° angle near the junctions of first-order confluences in Bristol, Florida. We then present a model for the migration in response to asymmetry in the groundwater field that predicts the observed deviation. Next, we show how the spectral analysis of streams can differentiate between clockwise and counterclockwise downstream flow, and we use this method to find a length scale of curvature fluctuation that agrees with the length scale of the confluence deviation. Finally, we find evidence for network-wide river migration, suggesting that forcing from groundwater may influence stream planimetry throughout the network.

## 2. Background

### (a) Groundwater flow

Groundwater-fed channels advance as re-emerging groundwater erodes material, causing the intersection point between the groundwater table and the ground to retreat, thereby extending the channel [8]. We begin by reviewing a theoretical framework that can describe these dynamics. Consider a lone channel advancing in a groundwater field. Pressure gradients from the surrounding groundwater table determine where and how the groundwater flows into this channel [18,19]. When these gradients are hydrostatic and flows are approximately horizontal, groundwater flux *q* (per unit length in the horizontal plane) can be simply described through Darcy’s Law in the Dupuit approximation [18,20,21]:

Here, *k* is the hydraulic conductivity and *h* is the height of the groundwater table. Mass conservation between flux out of the groundwater table and precipitation into it in steady state yields Poisson’s equation for *h*^{2} [21]:
*P* is the precipitation rate. Sufficiently close to a stream, precipitation *P* can be neglected, yielding Laplace’s equation [10,11]:

### (b) Growth direction

The solution to (2.3) can be combined with a growth hypothesis to predict how the channels grow forward. Recent results suggest that streams grow in a direction that can be equivalently described by three growth hypotheses: streams maintain symmetry of the groundwater table around the spring as they grow [12,13], streams grow in the direction of the streamline entering the spring [10], and streams grow in the direction of maximum flux entering the tip [13]. When a channel develops a confluence upstream, the streamlines enter the springs without curving when the angle between the confluent streams is 72° [10,11]. If streams advance at an angle narrower than 72°, the streamlines entering the two springs will bend away from each other; if streams advance at an angle wider than 72°, the streamlines entering the springs will bend towards each other [10,11]. The 72° confluence therefore can be considered a stable fixed point for stream advance [13]. An illustration of such a confluence and its streamlines is shown in figure 1.

This predicted shape agrees with the average angle at which confluent streams join in a river network in Bristol, Florida, determined by fitting individual streams to straight lines [10]. However, as we will demonstrate, real confluences deviate from this theory near stream junctions. We now proceed to briefly describe the field site used and methods employed to this end.

## 3. Confluence migration

### (a) Wide angle confluences in Bristol, Florida

Confluences are ubiquitous in river networks, but the often sedimentologically complex dynamics that form them thwart attempts to construct a unifying theory for their development. River confluence geometry can be influenced by a number of parameters: relative flow magnitudes and the depth of the deepest point along a stream bed (thalweg incision) [24], the overall power exerted by the network [25], fluid and sediment transport dynamics [26], scour depth [27], among others. However, these factors may not be the primary factors influencing migration in landscapes formed primarily through erosion by groundwater re-emergence.

We study streams from a network of groundwater-fed streams located on the Florida Panhandle. These rivers are overlain by homogeneous sediment comprising primarily unconsolidated quartz sand [22] which rests atop a layer of impermeable clay [28]. The high hydraulic conductivity *k*∼10^{−4} m s^{−1} of such clean, unconsolidated sand causes rainfall *P*∼5×10^{−8} m s^{−1} to be quickly absorbed, prohibiting surface run-off [22]. This supports a persistent water table in the subsurface. We first consider only first-order streams (those without any upstream confluences). These streams are shallow, with little evidence of scouring or deposition, negligible difference in bed height (bed discordance), minimal difference in discharge and negligible suspended load [28]; these properties make this field site ideal to study the influence of groundwater seepage in isolation from factors that typically influence confluence geometry. An example of such a confluence is shown in figure 2. These streams are roughly approximately 20 cm in width and approximately 1 cm in depth [29].

We extracted the location of streams in Florida by thresholding the curvature of topographic contours, obtained from laser altimetry data with a horizontal resolution of 1.2 m and a vertical resolution of approximately 5 cm [10], as outlined by Devauchelle *et al.* [29]. We ordered our channels according to Horton–Strahler ordering [30,31] and aligned first-order confluence pairs along their bisector, as determined by the sum of the vectors obtained from regression of each of the two confluent streams. To obtain the angle as a function of distance from the junction, the aligned channel coordinates were converted to polar coordinates, and the angle was averaged for each stream as a function of distance from the junction, as illustrated in figure 3.

A network-wide averaging of 1225 first-order confluences from our field site by this method is shown in figure 4*a*, where each point represents an average of points within a radial slice. Figure 4*b* quantifies the angle as a function of distance from the junction by transformation of the shape in figure 4*a* to polar coordinates. The shape we observe exhibits a wide angle deviation (approx. 120°) from the predicted 72° angle near the junction. We note that this averaging scheme is not necessarily representative of the archetypal groundwater-fed stream confluence shape, but a visualization of where channels are located on average.

There are three symmetric geometric adjustments that may give rise to this wide angle deviation: (1) headward growth in a direction that follows the deviant trajectory, (2) movement (advance) of the junction or (3) lateral migration of the channel after growth. We hypothesize that the lateral migration of channels is responsible for the deviation from 72° we observed in figure 4. The average channel location therefore represents a balance between migrational forcing and confinement by valley side walls. We therefore assume that figure 4 represents a statistical steady-state configuration.

To investigate these assumptions, we repeat the averaging scheme of figure 4 for channels sorted by arclength between spring and junction, precluding particularly small (less than 25 m) and large (greater than 120 m) channels from this analysis. Channels averaged in this manner are shown in figure 5, where different average channel lengths correspond to different colours. As channel length increases, the size of the deviation from the 72° lines increases (figure 5*b*), which suggests that the strength of confinement within valleys decreases with channel length. This is consistent with confinement caused by diffusion of the valley sidewalls, which are steepest near the springs. This observation will be discussed further when we talk about the characteristic scales associated with network-wide fluctuations in curvature.

### (b) A mechanism for lateral migration

We hypothesize that forcing from an asymmetric groundwater field downstream is responsible for this lateral stream migration, as streams adjust to the flux difference from their two sides. While self-similar branches with a 72° angle between them have a locally symmetric field around their tips, this field is not symmetric downstream from the tips. An illustration of this asymmetry is shown in figure 6*a*. Consider again the confluence geometry in which the two confluent streams are finite and mirror symmetric about a third, semi-infinite stream downstream. Let *α* be the angle between these two finite branches. Denote the region between them as the ‘inside’ and the regions between each branch and the semi-infinite branch as ‘outside’, as labelled in figure 6*a*. We reference these two regions with the subscripts *i* and *o*, respectively. The flux from the inside is therefore *Q*_{i}(*s*) and, from the outside, *Q*_{o}(*s*), where *s* is the arclength measured from the junction. We define the *flux jump* [*Q*(*s*)] as

Symmetry arguments provide some intuition about the geometry that migration according to [*Q*(*s*)] would produce. When *α*=72°, the field is locally symmetric at the tips [12] and hence [*Q*(*s*)]=0 at the tips. Therefore, no lateral migration would occur at the tips. Close to the junction, the two confluent streams appear infinitely long, becoming locally indistinguishable from the second-order, infinitely long downstream channel. When *α*=120° and all three branches are infinitely long, the field will be symmetric on either side of each channel, producing flux jump [*Q*]=0 everywhere. When the confluent streams are finite, [*Q*(*s*)] increases most slowly with *s* when *α*=120° (see appendix A for proof). Consequently, the flux jump near the junction is minimized when *α*=120°, and the angle *α*=120° can be understood as the most stable angle near the junction. We therefore hypothesize that real channels that migrate according to this flux jump exhibit an approximately 120° angle near the junction and an approximately 72° angle near the tips. These predictions agree with the angles of real confluences near the junction and tips, respectively, shown in figure 4*f*. This agreement suggests that these confluences rearrange primarily in response to the groundwater field rather than flow inside the channel. For example, a model of Howard [25] which optimizes the channel branching angle to minimize stream energy dissipation would predict a confluence angle of 90° for our network, contrary to the 120° angle we observe.

The flux-dependent migration process we have discussed thus far suggests that streams migrate with velocity *v*_{m} perpendicular to themselves at a rate which grows with the flux jump. Assuming a linear response, we then have
*s* with units of m; *Q*_{l} is the flux from the left of the stream when facing downstream and *Q*_{r} is the flux from the right. We use left and right here rather than inside and outside so that ^{2} s^{−1} divided by the precipitation rate. The quantity *v*_{m} is the perpendicular velocity at which the channel is migrating, *ϵ* is a constant with units *s*^{−1}, *ν* is a constant with units m^{2} s^{−1} and *κ* is the local signed curvature [32]. A negative curvature indicates that, when facing downstream, a channel rotates counterclockwise. The *νκ* term represents a stabilizing effect that can be understood as a first-order approximation to the influence of hillslope sediment transport processes as streams migrate against the valley wall.

If this stabilization is balanced by the flux difference across the stream and the streams are in a statistical steady state, we have
^{2} represents the difference in the area drained on either side over the length of the curve. The ratio of *κ* therefore provides a ratio *ν*/*ϵ* with units of m^{2}. We note that this steady-state assumption requires only that *ν* and *ϵ* be of the same sign. However, we expect that the sign of *ϵ* is positive, representing a migration towards the side of the stream that draws greater flux. We will provide evidence of this in §4b.

We calculate *κ* for a number of first-order confluences in our network, shown in figure 7*a*, and for all channels in our network, shown in figure 7*b*. We restrict these measures to a portion of the network for which we have LIDAR data that can clearly resolve curves in these streams. The slope of these lines can be used to obtain λ, which scales with the characteristic difference in distance to the drainage divide on either side of a stream, over the length of a curve in that stream. For all streams in our network, on average, λ∼30 m. As the characteristic distance to the drainage divide is approximately 100 m [11], this suggests that the field asymmetry produces a deviation in stream shape whose size corresponds to roughly 15% of the characteristic distance to the divide. We will now briefly explore the shapes of network-wide curves to explore precisely how this steady-state shape appears on the network scale.

## 4. Network migration

### (a) Network-wide power spectrum

To study how these stream curves appear in the network-wide regime, we measure the periodicity of network-wide fluctuations in curvature. A representative channel coloured by curvature is shown in figure 8.

Determining the periodicity of stream curvature requires a rotation-invariant, differential representation of the geometry. Traditionally, these periodicities are studied through spectral analysis of the angle or curvature [33,34]. However, while these methods are sufficient in isolating the primary frequency of curvature fluctuation, they ignore phase information that might reveal asymmetries in the growth. We therefore offer an alternative method: we rewrite our signal in terms of the 2D Frenet equations. The Frenet equations fully describe the relationship between tangent vectors, normal vectors and curvature on a twice continuously differentiable curve [35]:
** T** and

**are, respectively, the tangent and normal vectors to a curve with arclength**

*N**s*. We calculate the tangent vector over 3 m (as determined by the vector drawn between the two points surrounding the point in question) and curvature over 5 m for each point in all first-order channels, to obtain the d

**/d**

*N**s*. No additional information is carried in d

**/d**

*T**s*, so we can use d

**/d**

*N**s*without loss of generality. We write the tangent vector

**as a single imaginary number**

*T**T*

_{x}+

*iT*

_{y}, where

*T*

_{x}and

*T*

_{y}are, respectively, the

*x*and

*y*components of

**, and**

*T**i*

^{2}=−1. Equation (4.2) becomes

*κ*alone, but the sign of the frequency spectrum also indicates whether the channel is moving counterclockwise (positive) or clockwise (negative) as one travels downstream from the tip. The network-wide effect of this bias is reflected in the slight asymmetry of figure 9

*a*. Positive and negative frequencies are combined in figure 9

*b*and in subsequent figures.

We now calculate the power spectra separately for first-, second- and third-order channels and average them, as shown in figure 9*c*. We only use streams longer than 50 m for this calculation, to diminish the effect of the length distribution on the power spectra. Stream curvatures in our network appear to fluctuate at a dominant frequency that corresponds to a wavelength of 25.6 m for first-order channels, 30.8 m for second-order channels and 30.8 for third-order channels. The higher power spectral density of the second- and third-order channels results from higher values of curvature in second-order channels. However, because a higher value of curvature generates an increased wavelength by simply extending the channel arclength, the correspondence of first-order and second-order channel curvature wavelengths may be even closer than we have found. The agreement of this approximately 30 m wavelength to the length scale over which the wide opening angle of confluences returns to the 72° lines (seen in the intersection between the light blue and red lines in figure 4*f*) suggests that both network-wide fluctuations in curvature and confluence migration may be influenced by a common mechanism.

This spectral analysis resembles analyses conducted on river meanders, which form as flowing water erodes river banks and redeposits eroded sediment downstream [36–38]. However, we contend that the distinct morphology of the streams we study here suggests that they were shaped by different processes than those typically attributed to the formation of river meanders. For meanders confined by valley walls, the meander wavelength is approximately 17 times larger than the width of the meandering river [39]. By contrast, first-order channel widths in the Bristol, Florida, network range from approximately 10 to 20 cm near the springs up to 50 cm near downstream confluences, producing a ratio of wavelength to channel width of approximately 50–500. Streams in our network also have a sinuosity (defined as the ratio of channel length to valley length) less than 1.3, which would generally be classified as ‘straight’, rather than ‘meandering’ [3]. We therefore suggest that the fluctuations reflected in figure 9 are not primarily the result of erosional dynamics typically associated with river meanders, though the increased curvature in higher-order streams may signify an increasing importance of typical meandering dynamics downstream as confinement by valley side walls decreases. This behaviour is consistent with our observation in figure 5 that the size of the deviation of confluences from 72° increases with channel length. Finally, the independence of curvature wavelength with order suggests that this length scale is initially set by factors independent of dynamics within the streams and only later augmented by such dynamics, if at all. While we conjecture that straight streams therefore develop sinuosity as a result of the same processes described in §3, we do not know precisely how this occurs. Future work may shed light on the self-organization of these oscillations.

### (b) Evidence of network-wide migration

We now present field evidence that these channels are indeed migrating in the manner we suggest. We measure slope asymmetries along the valley side walls, which we define as the valley side wall slope on one side subtracted from the other. We hypothesize that if channels are laterally migrating, they should steepen the side wall which they migrate into.

We quantify the slope asymmetry by extracting the channel direction over a window of four metres along all channels in the network. We then obtain values of slope over 10 m on either side of the channel along the axis perpendicular to this direction from LIDAR data. We take the steepest 20% of these slopes and average them, under the assumption that the steepest sections of a valley represent the steepest steady-state configuration sustainable by a balance between migrational forcing and hill slope transport processes. Facing downstream, we call the averaged slope to the left of the stream *S*_{l} and to the right, *S*_{r}. We then calculate the slope asymmetry, which we define as (*S*_{l}−*S*_{r})/(*S*_{l}+*S*_{r}). This quantity is therefore positive when the steeper slope is on the left and negative when on the right.

We compare slope asymmetry, curvature and flux jump as defined in equation (3.3). The signs of curvature and slope asymmetry are illustrated in figure 10. The relationship between slope asymmetry and curvature is shown in figure 11*a*. A negative curvature indicates a negative second derivative: facing downstream, this corresponds to a channel that rotates counterclockwise. Negative values of curvature appear to correspond to positive values of slope asymmetry, suggesting that the valley side into which a channel veers has a higher slope. The data for flux jump against slope asymmetry, shown in figure 11*b*, also scale linearly. This suggests that the steeper side of the valley tends to draw greater groundwater flux. These results suggest that streams migrate, steepening the topography as they do.

## 5. Conclusion

We have presented evidence that stream migration can occur in response to an asymmetric groundwater table. We suggested that this migration is responsible for two planform features: a wide confluence angle and a network wide curvature fluctuation along the reach. First, we presented theoretical results that suggest that confluence migration according to this asymmetry would favour a symmetric configuration with a 120° angle near stream junctions, narrowing to a 72° angle near stream tips. We found evidence of such a feature in real streams. Second, we found that network-wide river curvature fluctuates at a length scale comparable to the length scale of this deviation, and we presented topographic evidence that suggests that this fluctuation is a product of the same migration mechanism. Our results reveal the presence of a novel mechanism in Laplacian growth that may inform understanding of other thin-finger systems.

## Data accessibility

The datasets supporting this article have been uploaded as part of the electronic supplementary material.

## Authors' contributions

R.S.Y., Y.C., O.D., G.G., H.J.S., and D.H.R. designed the research. R.S.Y. derived the theoretical results. R.S.Y., D.H.R., H.J.S., and G.G. analysed data. R.S.Y. and D.H.R. wrote the paper. All the authors gave their final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

This material is based upon work supported by the US Department of Energy, Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division under Award Number FG02-99ER15004.

## Acknowledgements

This work was inspired by valleys of the Apalachicola Bluffs and Ravines Preserve in Bristol, Florida. We thank K. Flournoy, D. Printiss, A. Schmidt, J. Stites and the Nature Conservancy for access to and guidance within the Preserve. We also thank Chris Follett and Eric Stansifer for their insightful comments.

## Appendix A. Optimal migration angle

**(a) Conformal map**

The general form of the conformal map which produces the fishtail shape in figure 1 is
*c* is the angle (*α*) between the two branches as a fraction of 2*π*. We see that *f*(*w*):*w*→*z* maps from the upper half-plane, *w*=*u*+*iv*, to the physical plane *z*=*x*+*iy*. Because the Laplace equation (2.3) is conformally invariant, the imaginary part of *f*^{−1}(*z*)=*g*(*z*):*z*→*w* is harmonic and provides us with the shape of the groundwater field around the bifurcated geometry. This map is illustrated in figure 12.

**(b) Flux**

While we cannot invert *f*(*z*) directly, we can calculate the flux jump implicitly in terms of the mathematical variable *v* by choosing a rotation of the map such that one of the branches lies along an axis. Conveniently, this is the case for (5.1), for which one branch lies along the positive *x* axis, as shown in figure 12. The flux into this branch is therefore d*v*/d*y*, where *v* is the imaginary part of *w*. To this end, the derivative d*z*/d*v* is given by
*y*/∂*v* is the imaginary component of this quantity evaluated at *v*=0:
*x*-axis we are concerned with, the argument within *Re* is always real. The flux *Q* within this domain is given by the reciprocal of this quantity:
*c* is shown in figure 6*c*.

**(c) Flux jump as x→0**

*Q* exhibits a peculiar behaviour: the flux *Q* and, consequently, the flux jump [*Q*] go to 0 regardless of the choice of *c*. To demonstrate this analytically, we can calculate the limit of (5.4) as *u*→0,
*π*). Similarly, for *c*>0 (the bifurcation angle is greater than 0). The flux jump therefore also approaches 0 for all values of *c* as *x*→0.

**(d) Flux jump departs from 0 most slowly for α=2π/3**

We seek to determine the rate at which the flux jump departs 0 from the junction as a function of *c*=*α*/2*π*. This can be found by solving ∂^{2}*v*/∂*x*∂*y*=0 for *c*. To this end, recall that ∂_{z}, the partial derivative with respect to *z*, and *x* and *y* as follows:
^{2}*v*/∂*x*∂*y* can therefore be rewritten in terms of *z* as
*v*/∂*z* in terms of *u* and *v*. We, therefore, require this expression to be rewritten in terms of ∂^{2}*z*/∂^{2}*v*, which we can now calculate directly:
*w*=*u*, for _{x}[*Q*(*x*)]=∂_{x}(*Q*(*u*_{1}(*x*))+*Q*(*u*_{2}(*x*))) as a function of *c* as *x* approaches its roots. Note that the flux jump here is a sum (as opposed to a difference) because the derivatives of *Q* with respect to *x* are signed. To obtain ∂_{x}[*Q*(*x*)], we invert the conformal map (5.1) in the limit as *u* approaches its roots, *v*=0 for equation (5.1). For small *u*, *x* is given by
*u* that correspond to a single value of *x*:
_{x}[*Q*(*x*)] explicitly for small *x*:
_{x}[*Q*(*x*)] for small *x*:

To understand the behaviour of this function, note that, within the range *α* and *β* are positive. Since we are concerned with the behaviour of (A 25) in the vicinity of *x* near the junction is at a minimum. This can be seen by considering the quantity
*x* behaves as
*ρ*(*c*) is shown in figure 13, and as expected by our reasoning above, *ρ* has a minimum at

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3918115.

- Received August 3, 2017.
- Accepted October 12, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.