## 1 Introduction

When rain falls on a sandy landscape, it flows into the aquifer, eventually re-emerging into streams. As groundwater flow accumulates in streams, sediment is removed and the landscape is eroded. This interplay between subsurface flow, erosion of the surrounding landscape, and the removal of sediment through streams causes existing stream heads to grow forward and new streams to form. In this way, groundwater-fed streams grow into ramified networks (Dunne, 1980).

Proposed examples of these so-called seepage channels can be found on both Earth (Orange et al., 1994; Schumm et al., 1995; Wentworth, 1928) and Mars (Higgins, 1982). Of these systems, a seepage network on the Florida Panhandle stands out as a superb example of the type. A topographic map of the network is shown in Fig. 1 (Abrams et al., 2009).

This kilometer-scale network is incised into a 65-m deep bed of sand overlying an impermeable layer of muddy marine carbonates and sands (Schumm et al., 1995). The sand layer corresponds to the ancient delta of the Appalachicola river, which was deposited from Late Pliocene to Pleistocene (Rupert, 1991; Schmidt, 1985). The incision of ravines into the sand plateau probably started about 1 My ago (Abrams et al., 2009). This origin explains the remarkable homogeneity of the lithology, only perturbed by localized clay lenses. The vertical localization of the streams seems not to depend on stratigraphy (Schumm et al., 1995). Groundwater flows through the sand above the impermeable substratum, into the network of streams, and drains into a nearby river (Abrams et al., 2009; Schumm et al., 1995).

The dynamics shaping groundwater-fed networks are simpler than the more common type of drainage networks that are fed by overland flow. In particular, the flow of the groundwater into streams is determined by the shape of the water table, which is a solution to a Poisson equation (Bear, 1979). Because the growth of streams fed by groundwater is an example of network growth in a Poisson field, the analysis of landscapes shaped by this process benefits from a substantial literature on interface growth in a harmonic fields (Carleson and Makarov, 2002; Hastings and Levitov, 1998; Saffman and Taylor, 1958). Past work (Petroff et al., 2011) has shown that the flow of groundwater into the Florida network is accurately described by the Poisson equation.

Here we use the specific example of the Florida seepage network to explore the general connection between the geometry and growth of drainage networks. This exploration takes the form of four thematically-related, but independent, exercises in which the dynamic simplicity of seepage erosion is exploited to develop a quantitative understanding of field observations. The focus of this collection is the realization that groundwater flow and landscape erosion are coupled to one another by the stream network, and the empirical validation that, close to equilibrium, this coupling takes a simple form. Related oddities and novelties are also briefly discussed.

## 2 Estimating groundwater flow from network geometry

We first investigate the relationship between a heuristic geometrical model for groundwater flow and the Dupuit approximation (Bear, 1979), which is based on well-defined physical assumptions. We begin by considering a spring, the head of a stream where groundwater first reemerges to the surface. In particular, we compare two models that relate the position of a spring in the network to the flow of water into it. In the field site, we consider here, the ground surface is permeable enough to instantaneously absorb most of the rainfall. As a consequence, the water is transported to the stream almost exclusively as groundwater (Abrams et al., 2009; Schumm et al., 1995).

Our first model is the continuum model alluded to in the introduction. According to the Dupuit approximation of Darcy's Law (Bear, 1979), the water flux **q** is related to the height h of the water table above the impermeable layer as

$$\text{q}=-Kh\nabla h\text{,}$$ | (1) |

$$\frac{K}{2}{\nabla}^{2}{h}^{2}=-P\text{.}$$ | (2) |

Because the change in elevation along the Florida network is small (the median slope S ∼ 10^{−2}), we approximate the height of the streams above the impermeable layer as constant. Because this equation is linear in the variable ϕ = Kh^{2}/2P, this boundary-value problem around the streams of the Florida network can be solved using the finite element method (Hecht et al., 2005; Petroff et al., 2011).

The second model of subsurface flow is a heuristic geometrical model (Abrams et al., 2009). According to this approximation, all groundwater flows towards the nearest point on the stream network. Thus, to find the flux of water into a section of the network, one first identifies the area a around the streams that is nearer to the selected section than to any other point on the streams. This drainage area is not related to the topographic drainage area, which is commonly used to model surface runoff. From conservation of mass, the discharge of groundwater into this section is Q_{g} = Pa.

Using these two models, we proceed to compare the estimated groundwater flux into 82 springs from the Florida network. To do so, it is useful to define the Poisson flux ${q}_{p}=\u2225\nabla \varphi \u2225$. Physically, q_{p} is the area draining into a section of stream per unit length. By analogy, the geometric flux into a section of network of length δs and drainage area a is q_{g} = a/δs. Because the area-driven model allows a finite area to drain into a point, this definition requires a finite value of δs; we take δs = 11 m. As shown in Fig. 2, there is a power law relationship (R^{2} = 0.94, p < 10^{−35}) between these quantities, ${q}_{g}=B{q}_{g}^{\alpha}$. The scaling exponent α = 0.69 ± 0.06. The proportionality constant, B = 4.0 ± 1.6 m^{1−α}, likely depends on the choice of δs.

Although there is no general quantitative relationship between these measures, there is a formal relationship. Let us imagine that each rain drop falls randomly, with a uniform distribution over the domain, and then undergoes a random walk until it reaches a boundary. The probability distribution of the random walkers’ positions satisfies the Poisson equation, which can be interpreted as a diffusion equation with a uniform source term (Doyle and Snell, 1984). As a consequence of its random walk through the domain, this imaginary rain drop has a continuous and non-vanishing probability distribution p of exit locations along the boundary, with a maximum at the closest point on the boundary. The mean flux of walkers through the section of channels gives the Poisson flux. To find the geometric drainage area, one sends each random walker to the mode of p. Thus, the precise relationship between q_{g} and q_{p} depends on the geometry of the system and, consequently, the scaling exponent α relating these quantities is likely not universal.

We conclude that the area-driven model can be safely used as a conceptual tool with which to gain intuition. Nevertheless, the Poisson equation should be used whenever a quantitative prediction of the seepage intensity or the relative fluxes into different parts of the network is required.

## 3 The three dimensional structure of the network

Having discussed the flow of water into a stream, we now consider the response of sediment within the stream to the flowing water. Flowing water erodes the bottom of the sandy stream when the shear force exerted on a sand grain is sufficient to overcome friction. Thus, there is a threshold that the shear must exceed for any sediment to be transported. This threshold is often expressed as a function of the Shields parameter (which compares the grain weight to the flow-induced shear) in bedload transport laws (Sheilds, 1936; Vanoni and Keck, 1964). Many streams, and the Florida network in particular, are thought to adjust to a state in which every grain is on the threshold of motion (Devauchelle et al., 2011; Henderson, 1961; Lane, 1955). If we further approximate the stream with a straight and shallow channel formed in uniform sediment, the above constraint on the shear can be re-expressed as a relationship between the slope of the stream bed S and the discharge Q (units of volume/time) of water in the stream as

$$Q{S}^{2}={Q}_{0}\text{,}$$ | (3) |

_{0}is a constant with units of discharge the value of which depends on parameters such as the grain size (Devauchelle et al., 2011; Henderson, 1961; Lane, 1955). Thus the threshold theory provides a new interpretation to the typical relationship between slope and area (Montgomery, 2003). Using the conservation of mass and momentum, this equation becomes equivalent to Lacey's relation, which states that the width of a stream scales like the square root of its discharge (Lacey, 1929).

Because both the discharge of a stream and the slope affect the flow of groundwater into it, Eq. (3) imposes a boundary condition for the flow of groundwater into the streams, from which one can derive the profile of an isolated stream (Devauchelle et al., 2011). Here we show that this boundary condition is satisfied throughout the network.

To test the applicability of Eq. (3), we compare the measured slope of streams to the stream discharge. Estimates of both the slope and the discharge require a three-dimensional description structure of the drainage network, which is extracted from the high-resolution topographic map of the network shown in Fig. 1. A depiction of part of the network is shown in Fig. 3. Given the height H_{s} of each point of the network, the slope is measured from the change in height along the stream. To estimate Q, we first solve Eq. (2) subject to the boundary condition on the streams

$$h={H}_{s}\text{.}$$ | (4) |

Because Eq. (2) is a non-linear function of h, its solution depends on the height of the streams above the impermeable mudstone layer. The elevation of the impermeable layer above sea-level is h_{0} = 0. Fig. 4 shows the height of the water table around the network taking P = 5 × 10^{−8} m/s, K = 1.6 × 10^{−5} m/s. The value of the precipitation rate P corresponds to the yearly average rainfall in the concerned area. The value of the conductivity K was adjusted independently to match the distribution of discharges measured in the field (Petroff et al., 2011). Given this solution, Q is estimated throughout the network as

$$Q({x}_{0},{y}_{0})=-K\underset{N({x}_{0},{y}_{0})}{\oint}h[\frac{\partial h}{\partial n}]\text{d}s\text{,}$$ | (5) |

_{0}, y

_{0}) is the network upstream of the point (x

_{0}, y

_{0}), s is a curvilinear coordinate along the streams, n is the local normal, and [·] represents the “jump”, which accounts for the flow of groundwater into either side of the one-dimensional stream.

Fig. 5 shows the comparison between the measured slope-discharge relation (blue points) and the quadratic relation predicted by Eq. (3). The value of Q_{0} is fit to the data and corresponds to a grain size of d_{s} = 1.7 mm, consistent with observation. The slight disagreement between theory and observation at high and low values of discharge is likely the result of errors in the extraction of the network from the elevation map. The measurement of the discharge close to a spring is sensitive to the position and elevation of the spring in the network. The measured elevations of a substantial fraction of the springs are above the solution of the water table elevation, leading to a negative value of Q very close to the spring. The measurement of very small values of the slope is also affected by error in the original topographic map. Finally, the dependence of the friction coefficient on flow depth or the downstream fining of sediments can explain the slight overestimation of the exponent by the simplest version of the threshold theory.

The approximate agreement between theory and observation leads us to conclude that the streams throughout the network are poised at the threshold of sediment motion. The solution of Eq. (2) subject to boundary condition (3) determines the height of the streams throughout the network (Devauchelle et al., 2011).

## 4 Response of the surrounding topography

In the Florida network, streams are incised into a bed of unconsolidated sand. Since there is no surface runoff except for the streams themselves, the topography around them evolves as a function of slope only. The dependence of the sediment flux on slope is likely to be non-linear where the slope is close to critical. Far from the channels, however, we can represent the relaxation of the landscape by linear diffusion (Culling, 1960). We have not evaluated this hypothesis, but this assumption allows us to illustrate how the effect of sunlight on diffusion can be quantified. According to this hypothesis, the sediment flux **j** is related to the height of the topography H as

$$\text{j}=-D\nabla H\text{,}$$ | (6) |

$$\frac{\partial H}{\partial t}=D{\nabla}^{2}H\text{,}$$ | (7) |

In general, D is determined by processes including rainfall, animal activity and vegetation. In this section, we consider variations in the diffusion coefficient (Bass, 1929; Burnett et al., 2008; Emery, 1947; Istanbulluoglu et al., 2008).

As shown in Fig. 6, the valleys cut by the drainage network are asymmetric; the southern valley wall is steeper than the northern wall. This asymmetry can be interpreted in two ways. Either there is a substantial north-south variation in the rate the landscape is diffusing or the the erosion rate is symmetric but the value of D is different. Because this network is in the Northern Hemisphere, the average position of the sun is to the south of the channels. Consequently, the southern walls of the valley are more shaded than the northern walls. Past studies have shown that there is a difference in vegetation between the northern and southern sides of the valley (Kwit et al., 1998).

We characterize the influence of the sun with the projection of sun light onto the landscape. If $\stackrel{\u02c6}{\sigma}$ is a unit vector pointing towards the sun and $\stackrel{\u02c6}{N}(x,y)$ is the unit vector that is normal to the landscape at the point (x, y), the dimensionless light intensity is

$$I(x,y)=\stackrel{\u02c6}{N}(x,y)\cdot \stackrel{\u02c6}{\sigma}\text{.}$$ | (8) |

It is straightforward to measure $\stackrel{\u02c6}{N}$ from the topographic map shown in Fig. 1. Given the latitude and longitude of the Florida network, the annual average solar zenith and azimuthal angles are 26.6° and 180° respectively (Reda and Andreas, 2004). Combining these values with Eq. (8) provides an estimate of I at each point in the field site.

We hypothesize that differences in average light intensity give rise to differences in D, resulting in asymmetric valley walls. Expanding D to first order in I gives

$$D\approx {D}_{0}{I}_{0}(1+{D}_{1}{I}_{1})\text{,}$$ | (9) |

_{0}is the mean diffusion coefficient, I

_{0}is the mean light intensity, I

_{1}= (I−I

_{0})/I

_{0}is the fluctuation in the light intensity, and D

_{1}is the correction due to light.

Here we estimate D_{1} from the shape of the topography assuming that there is no north-south asymmetry in the erosion rate. It is useful to express $j=\u2225\text{j}\u2225$ and I_{1} in terms of the orientation ω. Here, ω is defined as

$$\mathrm{cos}\text{\hspace{0.28em}}(\omega )=\frac{\nabla H}{s}\cdot \stackrel{\u02c6}{x}\text{,}$$ | (10) |

_{1}, estimated throughout the field site, with orientation is shown in Fig. 7. The distribution of slopes with respect to orientation shows two modes. The mode with period 2π is in phase with sun light intensity. The mode with period π/2 is caused by favored orientation of channels in the network, visible on Fig. 4. Indeed, a slope facing one of the cardinal directions is more likely to belong to a channel wall than if it were facing another direction. Now, channel walls have a much steeper slope than the rest of the topography, hence the higher average slope associated to cardinal directions on Fig. 7.

If there is no net asymmetry in the erosion rate, then

$$\underset{0}{\overset{2\pi}{\int}}j\left(\omega \right)\mathrm{sin}\left(\omega \right)d\omega =0\text{,}$$ | (11) |

**j**or in relation to the first Fourier coefficient. Recognizing that j = D

_{0}I

_{0}(1 + D

_{1}I

_{1})s, it follows that

$${D}_{1}=-\frac{\underset{0}{\overset{2\pi}{\int}}{I}_{1}(\omega )s(\omega )\mathrm{sin}(\omega )d\omega}{\underset{0}{\overset{2\pi}{\int}}s(\omega )\mathrm{sin}(\omega )d\omega}\text{.}$$ | (12) |

Combining the measured values of s and I from the Florida network with Eq. (12), yields D_{1} = 0.53 ±0.02. Thus, the fractional variation in the diffusion coefficient needed to account for the slope asymmetry is ${D}_{1}{\langle {I}_{1}^{2}\rangle}^{1/2}=0.048\pm 0.001$.

## 5 Growth of a stream

Having discussed the flow of groundwater into streams and the resulting adjustment of stream shape, and erosion of the landscape, we now turn our attention to the growth of a spring in response to the flow of groundwater. Because springs in the Florida network grow slowly, at a speed of ∼ 1 mm/year (Abrams et al., 2009; Schumm et al., 1995), this section relies on experimental observations.

Seepage channels are grown in a previously described experimental apparatus (Lobkovsky et al., 2007). In this experiment, a hydraulic head of 19.6 cm is used to push water through a bed of sand (grain diameter of 0.5 mm), having an initial slope of 7.8°. As water flows out of the sand bed, it entrains sand grains, thus forming a seepage channel. The channel is initialized with a small indentation in the otherwise flat bed of sand.

We characterize the growth of a channel by the shape of an elevation contour. A map of the experiment showing the shape of an elevation contour and its position in the experiment is shown in Fig. 8. Given this characterization, we ask how the growth of the contour outwards depends on the flux of water into it. To this end, we compare velocity of the contour with the solution of Eq. (2) around the channel.

The normal velocity of the contour is measured by comparing the shape of the contour to its shape a moment later. Given an elevation, we determine the normal vector at each point along the associated contour. We then compare the shape of a contour to the shape of the same contour three minutes later to determine the amount that each point on the contour grew in the normal direction. Two elevation contours at three minute intervals are shown in Fig. 9.

To determine the flux of water entering each part of the channel, we solve for the shape of the water table. In these experiments, there is no rainfall, thus Eq. (2) simplifies to the Laplace equation

$${\nabla}^{2}\psi =0\text{,}$$ | (13) |

_{c}) relative to the size of the hydraulic head Δh = 19.6 cm. By construction, ψ = 0 on the channel and ψ = 1 at the back wall. Because Eq. (13) is linear, this boundary-value problem can be solved using the finite element method (Hecht et al., 2005). We solve for the shape of the water table around a growing elevation contour at one minute intervals. Fig. 9 shows the dimensionless flux $q=\ell \u2225\nabla \psi \u2225$, where ℓ = 90 cm is the length of the experiment, into different parts of the channel.

As shown in Fig. 10, the velocity at which the channel grows outward is linearly related to the flux of water into it. This is unexpected, as typical transport relationships relate sediment transport to shear nonlinearly. Moreover, according to the least-squares fit of the data, the intercept of this relationship is negative, meaning that a finite flux of water is required for the channel wall to grow forward. This result is superficially similar to the finite shear required to transport sediment in a stream (Sheilds, 1936; Vanoni and Keck, 1964). Because the growth of the channel requires that sediment be transported through the channel, this relationship between flux and growth does not reflect a simple force balance on a sand grain. Rather, this relationship arises from the coupling between sediment transport within the channel and subsurface flow around the channel. It is therefore remarkable that this relationship takes a form in which the growth at a point depends only on the flux at that point. As the channel wall height remains roughly constant during the experiment, this result implies that the increase in the network volume is roughly proportional to the total delivery of water (neglecting the channel sides, where the water flux is small).

## 6 Discussion

This collection of remarks has been arranged so as to connect – the way a skipping stone connects the banks of a wide river – the flow of groundwater through the aquifer to the flow of sediment over the landscape. We began, in section 2, by relating a heuristic model of subsurface flow to Darcy's law. In section 3, we showed that the boundary condition for the flow of groundwater is determined by the condition that sand grains in the streams are at the threshold of motion; this boundary condition determines the height of streams (Devauchelle et al., 2011). In section 4, we proposed a method to measure the influence of sun exposure on landscape diffusion from topographic data. Finally, in section 5, we showed that the local erosion rate of experimental seepage channels is linearly related to the water flux, beyond a threshold value.

We now pause to consider the metaphoric “wide river”. The first three results represent three aspects of a single boundary value problem specified by the principal equations of these sections,

$$\frac{K}{2}{\nabla}^{2}{h}^{2}=-P\text{,}$$ | (14) |

$$\frac{\partial H}{\partial t}=D{\nabla}^{2}H\text{,}$$ | (15) |

$$Q{S}^{2}={Q}_{0}$$ | (16) |

$$H=h\text{.}$$ | (17) |

Given the current shape of the landscape, these equations describe the evolution. It is an open question if this set of equations describes the advance of the spring. In particular, we ask if the observed relationship between channel growth and groundwater flux in section 5 represents a fourth aspect of this boundary value problem.

The answer to this question deeply influences our conceptualization of how drainage networks grow. If these equations are complete, network growth represents a balance between the accumulation of water in the streams and the erosion of the surroundings. In this case, sediment transport is simply the mechanism that maintains this balance. If these equations are not complete, the way a network grows and the surrounding landscape changes depend on the details of how sediment moves through the network. Once this question is answered, the details of this model can be adapted to suit the details of a more complex world.

## Acknowledgments

We would like to thank The Nature Conservancy for access to the Apalachicola Bluffs and Ravines Preserve, and K. Flournoy, B. Kreiter, S. Herrington and D. Printiss for guidance on the Preserve. Thanks to D.M. Abrams and H. Seybold for fruitful discussions. We thank B. Smith for her experimental work. Finally, we thank the referees C. Paola, P. Davy, S. Douady and Y. Couder for their careful reviews, which helped us improve and clarify the manuscript.

This work was supported by Department of Energy grant FG02-99ER15004. O. D. was additionally supported by the French Academy of Sciences. A.K. was additionally supported by Department of Energy grant FG02-02ER15367.