Controls on Zero‐Order Basin Morphology

Zero‐order basins are common features of soil‐mantled landscapes, defined as unchanneled basins at the head of a drainage network. Their geometry and volume control how quickly sediment may reaccumulate after landslide evacuation and, more broadly, zero order basins govern the movement of water and sediment from hillslopes to the fluvial network. They also deliver water and sediment to the uppermost portions of the fluvial network. Despite this role as the moderator between hillslope and fluvial processes, little analysis on their morphology has been conducted at the landscape scale. We present a method to identify zero‐order basins in landscapes and subsequently quantify their geometric properties using elliptical Fourier analysis. We deploy this method across the Coweeta Hydrologic Laboratory, USA. Properties such as length, relief, width, and concavity follow distinct probability distributions, which may serve as a basis for testing predictions of future landscape evolution models. Surprisingly, in a landscape with an orographic precipitation gradient and large hillslope to channel relief, we observe no correlation between elevation or spatial location and basin geometry. However, we find that two physiographic units in Coweeta have distinct zero‐order basin morphologies. These are the steep, thin soiled, high‐elevation Nantahala Escarpment and the lower‐gradient, lower‐elevation, thick soiled remainder of the basin. Our results indicate that basin slope and area negatively covary, producing the distinct forms observed between the two physiographic units, which we suggest arise through competition between spatially variable soil creep and stochastic landsliding.

Several models of hollow and zero-order basin evolution have been proposed, considering the buildup and evacuation of material, the potential for basins to laterally migrate, or to appear as static or transient landscape features. Dietrich and Dorn (1984) review many of these theories, which argue that hollows form as a consequence of gullying or deep-seated landsliding, which both create topographic convergence. Bryan (1940) suggests that following this initial formation the base of the gully is armored, resulting in future gullying, which erodes hollow noses and produces a steady rate of hillslope retreat in a process termed gully gravure. From a similar observation of the potential for bedrock armoring, Mills (1981) argues that the oversteepening  Hack and Goodlett (1960), which are used to define a zero-order basin (composed of a hollow, planar side slopes, and divergent noses) in a landscape. Profile follows the path of the red dashed line (A-A ′ ) in (a). of the evacuated hollow coupled with the armoring of the base of the hollow will lead to the lateral migration of zero-order basins. Work by Stock and Dietrich (2006) in the western United States demonstrated that zero-order basins can evolve due to the exposure of bedrock by debris flows and the subsequent preferential weathering of these surfaces.
Zero-order drainage basins have long been recognized as important for understanding the patterns of discharge within mountain catchments (Hewlett & Hibbert, 1967) and potentially across large drainage basins (MacDonald & Coe, 2007). Hydrologists have suggested that the storage of water within zero-order basins controls the flashiness of discharge within mountain catchments (Tsukamoto, 1963) and that the size and shape of zero-order basins can strongly affect flow regimes (e.g., the variable source concept; Hewlett & Hibbert, 1967). One important component of the hydrological regime that is controlled by the shape of zero-order drainage basins is the position of the channel head, which represents a point in space where geomorphic and hydrologic processes transition from hillslope (and throughflow) to fluvial (and channel flow) processes (Clubb et al., 2014;Jefferson & McGee, 2013;Julian et al., 2012;Montgomery & Dietrich, 1988, 1989Tarolli & Dalla Fontana, 2009). Defining the location of the channel head (and therefore the bottom of the zero-order drainage basin) has been challenging due to the complex interplay of hillslope processes, overland flow, and fluvial processes (Dunne et al., 1991).
Zero-order drainage basins are the setting for a range of geomorphic processes occurring within several linked topographic domains (hollows, side slopes, and noses; ; Hack & Goodlett, 1960). A zero-order basin contains a colluvial hollow, a zone of topographic convergence, bounded by zones of topographic divergence termed hillslope noses, which pass into planar side slopes that drain to the hollow axis ( Figure 1). Noses and side slopes are dominated by creep-like processes, such as tree throw and animal burrowing, whose sediment flux to the hollow axis is either a linear or nonlinear function of slope and is mediated by sediment thickness (Heimsath et al., 2005). Sediment accumulates within the hollow axis for tens to thousands of years, before being removed by one of a number of processes that include overland flow erosion , shallow landsliding (Alger & Ellen, 1987;Benda & Dunne, 1997;Dietrich & Dunne, 1978;D'Odorico & Fagherazzi, 2003;Hack & Goodlett, 1960;Parker et al., 2016;Reneau et al., 1989), and the upslope migration of transient channel heads into the zero-order basin (Dietrich & Dunne, 1993). Below the channel head, erosion and deposition occur by a mixture of debris flows and fluvial processes . The complex interaction of geomorphic processes occurring between zero-and first-order drainage basins has limited the development of simple erosion laws ) that can explain the shape of these features.
Given the importance of the size, area, and topography of zero-order drainage basins for the hydrology and erosion of catchments, it is surprising how little we understand about the shape of these features. When most authors conceptualize a zero-order basin, it has an approximately teardrop-shaped planform (e.g., Dietrich & Dorn, 1984;Reneau et al., 1986Reneau et al., , 1989Thorne et al., 1987), yet this conceptualization has never been systematically tested. The small number of field observations that have been collected suggest that zero-order basins can have complex shapes including branching forms . This paper seeks to systematically understand how we might measure the shape of zero-order drainage basins and apply our measurements to an Appalachian Mountain catchment, chosen because the Appalachians feature well-defined ridge and valley topography (Hack & Goodlett, 1960) and in addition our specific study catchment, at Coweeta North Carolina, has been the site of many decades of studies into both hydrolgic connectivity (e.g., Band et al., 2012;Hewlett & Hibbert, 1967) and sediment transport (e.g., Hursh, 1941;Wooten et al., 2008Wooten et al., , 2007. We define and model the distribution of zero-order drainage basins using the following criteria: (i) they are areas of topographic convergence high on hillslopes, which concentrate flows of water and sediment into their apex (Dietrich et al., 1982;Hack & Goodlett, 1960;Reneau et al., 1986;Sólyom, 2011), with a defined point of maximum convergence within a hollow (Parker et al., 2016). (ii) There is no stream channel present, with the channel head located at or below the base of the zero-order basin. (iii) The feature must contain convergent topography separated by areas of divergent topography. Our goal is to produce reproducible information about zero-order basins that can be compared to model predictions in order to gain insight into geomorphic processes.

Coweeta, North Carolina
The Coweeta Hydrologic Laboratory in North Carolina ( Figure 2) has been managed as an experimental forest since 1935 and as part of the Long-Term Ecological Research network since 1980 (Douglass & Hoover, 1988). The catchment covers 50 km 2 of the steeper, eastern side of an asymmetrical escarpment (locally called the Nantahala Mountains Escarpment; Wooten et al., 2008) with considerable relief (∼900 m over the 5-km length of the basin) and a strong orographic precipitation gradient (a mean annual precipitation range of 1,800 to 2,300 mm; Swift et al., 1988). The Coweeta basin has a physiography that is typical of basins which drain the Blue Ridge Escarpment (Prince et al., 2010), with steep upper slopes along a linear escarpment, and lower slopes demonstrating a more typical ridge and valley topography. In Coweeta these two broad physiographies are represented by (i) a steep, thinly soil mantled, northern hardwood-dominated escarpment, and (ii) a lower-gradient, lower elevation catchment with thick soils and saprolites, dominated by oak-pine and cove hardwood forest. Consequently, these two physiographic regions can be used as natural laboratories to explore the factors which control zero-order basin formation, evolution, and evacuation.
Coweeta is underlain by high-grade metamorphic rocks, primarily biotite and quartz diorite gneisses, with minor schists and metasedimentary rocks that have been strongly folded and faulted throughout the basin (Hatcher, 1979). Authors have suggested that this area may diverge from steady state due to rejuvenation of uplift in the Miocene (Gallen et al., 2013) or drainage reorganization associated with escarpment retreat (Prince & Spotila, 2013). However, erosion rates measured in the basin are consistent with rates measured in the nearby Great Smoky Mountains, suggestive of a mountain range approaching steady state. The basin topography is characterized by the nose and hollow topography typical of Appalachian uplands (Hack & Goodlett, 1960; see Figure 1 for a schematic of a typical zero-order basin cross section). Field observations of tree throw and animal burrowing in Coweeta suggest that sediment transport into hollows is dominated by creep processes, whereas evacuation of hollows is predominantly by landsliding, consistent with historical observations (Hursh, 1941;Wooten et al., 2008Wooten et al., , 2007. Coweeta is likely to have maintained a soil mantle through much of the Quaternary, with three main forest assemblages being present since the last glacial maximum (LGM): boreal forests on upper slopes during the LGM, northern hardwood forests (dominated by sugar maple and oaks) at the highest elevations today, and cove hardwoods (dominated by chestnut, hemlock, and tulip poplar) in lower elevations today (Delcourt et al., 1982). During the LGM, there may have been some periglacial activity on the upper slopes (Braun, 1989;Clark & Ciolkosz, 1988). The thick forest mantle remained intact through to the late nineteenth century, when a short period of clearfelling and selective logging occurred on Coweeta's hillslopes before establishment of the Hydrologic Laboratory in 1934 (Douglass & Hoover, 1988). An increase in sedimentation rates have been recorded in alluvium in the area (Wang & Leigh, 2015), but without producing an obvious increase in landsliding rates (Eschner & Patric, 1982;Parker et al., 2016). The short period of increased sediment flux is unlikely to have significantly affected hillslope morphology. Observed differences in species composition through the catchment have been demonstrated to have little effect on the distribution of root biomass (Hwang et al., 2015) or root cohesive strength (Hales & Miniat, 2017;Hales et al., 2009).

Methods
Here we present a new method for extracting zero-order basin morphologies from high-resolution topographic data through the identification of the maximum upslope extent of the channel network and the extraction of hydrologically connected patches of hillslope draining into these points. These patches are then converted to vector outlines for further analysis and processing.

Point Cloud Processing
Topographic data were generated for our study site directly from point cloud data, provided by OpenTopography, which has been manually postprocessed to distinguish between vegetation and ground returns within the point cloud. This is a vital step in the generation of digital elevation models (DEMs). The accurate separation of the point cloud into ground and nonground returns is particularly important in our heavily vegetated study sites. Automated routines have been developed to process unfiltered point clouds (e.g., Evans & Hudak, 2007); however, it has been demonstrated that as the vegetation cover increases, more postprocessing and interpolation is needed to generate a final surface (e.g., Liu, 2008;Meng et al., 2010;Tinkham et al., 2012). Therefore, we use point clouds that have undergone a supervised classification in conjunction with automated filtering, both performed by the National Center for Airborne Laser Mapping (http://www.ncalm.org). The accuracy of such filtering has been evaluated in Santa Cruz Island, CA, where Perroy et al. (2010) demonstrated a vertical accuracy of 0.067 m at a 95% confidence interval across the Island, suggesting that with this type of processing the surface morphology can be accurately represented at the scale with which we intend to measure features, even considering the potential for increased vertical inaccuracy due to Coweeta's dense forest cover.
The Coweeta point cloud has a point density of 8.91 points per m 2 , a reported vertical accuracy of 0.13 m and a horizontal accuracy of 0.11 m. It has been shown in previous studies that the Coweeta point cloud is able to support 1-m resolution data (Grieve, Mudd, Hurst, & Milodowski, 2016;, and as such, the data set was gridded to a 1-m cell size. The classified point cloud is gridded using a local binning algorithm, which searches for points within a circular window with a radius defined by Kim et al. (2006) as where R is the desired output resolution. Within each circular window an inverse distance weighting is employed to all of the points found, calculating the elevation value for each grid cell. Different methods of scaling the search radius have been proposed, but Kim et al. (2006) suggest that equation (1) is the most parsimonious solution, which will yield a continuous surface in most cases. This gridding method has been used successfully in previous studies Grieve, Mudd, Hurst, & Milodowski, 2016; and produces a surface with few data gaps and thus little need for interpolation, which may impact the quantification of zero-order basins properties. When data gaps are present, they are filled using an inverse distance weighted focal mean of the surrounding cell elevations.

Hydrological Correction
To extract zero-order basins from topography the channel network must first be defined, as models of ridge-hollow topography typically identify the colluvial-fluvial transition as the lowest point of a hollow (Reneau et al., 1986). To perform this, we use a modified implementation of the DrEICH algorithm (Clubb et al., 2014) parameterized to extract the bases of zero-order basins rather than the channel network.
In order to perform any hydrological analysis on a DEM, it must first be hydrologically corrected to ensure that cells with no downslope neighbors, known as pits or sinks, are removed (Mark, 1984). In early topographic analysis such pits were filled using multiple passes of a smoothing window (O'Callaghan & Mark, 1984) until the pit was removed, but this method fundamentally changes the morphology of the surface, particularly in areas of high topographic complexity (Band, 1986), which are particularly evident in high-resolution topographic data (Purinton & Bookhagen, 2017).
In higher-resolution data, pits can be removed through the use of a constructive algorithm, which increases the elevation of each pit until flow can pass across it unobstructed (e.g., Jenson & Domingue, 1988;Tarboton, 1997); a destructive algorithm, which reduces the elevation of points surrounding a pit until the pit is removed (e.g., Martz & Garbrecht, 1998); or by a combination of destructive and constructive methods, designed to limit the amount of modification of the topographic surface (Lindsay & Creed, 2005;Soille, 2004).
Here we implement the optimized constructive algorithm of Wang and Liu (2006), which utilizes a priority queue to incrementally identify the outlets of depressions from the DEM edge, upslope toward the center cells. Pits are filled until they reach a threshold gradient in the downslope direction of 0.0001, selected as it ensures realistic surface flow across filled pits with minimal change to the topographic surface. Minimal alteration of the topographic surface is produced using this method, and it has been used successfully for several geomorphic applications in a diverse range of settings (Clubb et al., 2014(Clubb et al., , 2017Mudd et al., 2014;Grieve, Mudd, Hurst, & Milodowski, 2016;Milodowski et al., 2015).

Surface Fitting
To extract the channel network, the tangential curvature of the surface must be extracted in order to identify potentially channelized portions of the landscape (e.g., Pelletier, 2013). Hurst et al. (2012) demonstrated the suitability of extracting surface derivatives including curvature from high-resolution topographic data by fitting a quadratic function of the form to elevation values within a moving circular window, using least squares regression (Evans, 1980). Where x and y are horizontal coordinates, is the elevation and a, b, c, d, e, and f are fitting coefficients. The scale of this window is selected by identifying scaling breaks in the standard deviation and interquartile range of curvature as the window size is increased (Hurst et al., 2012;Lashermes et al., 2007;Roering et al., 2010). Using such a window size ensures that measurements of tangential curvature from the polynomial surface are not influenced by microtopographic variations that may be generated by a combination of natural processes roughening the landscape such as animal burrowing or tree throw, or from measurement noise generated during lidar data capture (Hurst et al., 2012;Roering et al., 2010), but rather represent the hillslope scale morphology with which we are concerned. The fitted coefficients of equation (2) can be employed to calculate the tangential curvature, given by Mitasova and Hofierka (1993) as In addition to calculating the tangential curvature, total hillslope curvature and gradient can also be calculated using the fitted coefficients and each of these measurements are used to characterize zero-order basins throughout this study. Total curvature, C Total , is calculated as 10.1029/2017JF004453 and topographic gradient, S, is calculated as,

Zero-Order Basin Extraction
To extract zero-order basins, we must find the transition between hillslopes and channels. It is possible to use a drainage area threshold for this step, but this method has been shown to be particularly problematic in high-resolution topography when contrasted with field-mapped channel heads (Clubb et al., 2014;Orlandini et al., 2011;Passalacqua et al., 2010;Pelletier, 2013;Sofia et al., 2011). Instead, we begin by computing a provisional network based on using planform curvature to extract valley heads, which we have previously shown to be robust at a range of data resolutions . We first employ a Wiener filter (Wiener, 1949) to remove noise from the hydrologically corrected DEM (Pelletier, 2013). A tangential curvature threshold is then identified through analysis of its quantile-quantile plot to identify the point at which curvature values deviate from a normal distribution (e.g., Lashermes et al., 2007;Passalacqua et al., 2010). This curvature threshold is then employed to identify discrete channelized portions of the landscape.
These channelized patches of the landscape are then merged into a contiguous channel network by employing a connected components algorithm (He et al., 2008), which merges these discreet patches into a contiguous channel network. A threshold of 10 pixels is applied to ensure that small patches of positive tangential curvature caused by surface roughness are excluded. This connected components network is then thinned to a single pixel wide network using the algorithm of Zhang and Suen (1984) and the upstream limits of this network are identified as the input channel heads for the DrEICH algorithm (Clubb et al., 2014).
The final processing step for these input channel heads is to perform a steepest descent flow routing algorithm (O'Callaghan & Mark, 1984) to generate a channel network, and identifying any channels which are composed of a single pixel. Such single-pixel channels are considered to be a product of a combination of artificial and natural topographic noise common in high-resolution topography (e.g., Roering et al., 2010;Sofia et al., 2013) and are subsequently removed from the input channel heads prior to their use in the DrEICH algorithm.
The DrEICH algorithm has been evaluated against field-mapped channel heads in several locations of varying geomorphic character (Clubb et al., 2014). This algorithm identifies channel heads as the point at which the topographic signal transitions from fluvial to hillslope-dominated processes. This is performed through the transformation of traditional river profiles by integrating over drainage area (Perron & Royden, 2013;Royden et al., 2000). Such transformed profiles, termed plots, produce linear profiles when is plotted against elevation. Mudd et al. (2014) developed a statistical technique, which identifies best fit linear segments within these plots to facilitate the identification of landscape transience. When the technique is applied to nonfluvial topography, the segments become increasingly nonlinear. Therefore, the algorithm of Clubb et al. (2014) is designed to identify the transition between linear and nonlinear -elevation plots, with the spatial location of this transition identified as the channel head.
We modify the input parameters for the DrEICH algorithm to only consider first-order channels and to identify the uppermost signal of fluvial incision upon a hillslope, thereby extracting the bases of zero-order basins from high-resolution topography rather than a network of channel heads initiating at the point where fluvial incision dominates over hillslope processes. This results in a series of points on the landscape identified as the transition between hillslope and fluvial processes which we can define as the base of our zero-order basins.
The location of channel heads in the landscape is challenging to identify, particularly as such features may be dynamic or transient (e.g., Dietrich & Dunne, 1993). In some landscapes it is also challenging to distinguish between channeled and unchanneled valleys, carved by debris flow action, occurring far above the fluvial network (e.g., Dunne et al., 1991;Penserini et al., 2017). However, the relatively small number of recent debris flows in the area, and the challenge of identifying debris flows even after recent debris flow events (Band et al., 2012), suggests that channels heads represent the lower limit of hillslope processes. Clubb et al. (2014) found good agreement between field-mapped channel heads and the DrEICH algorithm's results, although there is still potential uncertainty within the extracted zero-order basin bases. To attempt to constrain these potential location errors, two further sets of bases are produced. The initial set of bases produced using the DrEICH algorithm represent the upper limit of convergence in a landscape. Attempts to extract zero-order basins from above these points typically fail as the amount of hillslope scale convergence is decreased relative to smaller scale topographic disturbances. As a consequence of this the upslope areas extracted from points above this base are not representative of the ridge-hollow geometry of the landscape in question. Therefore, a steepest descent algorithm is used to move the base downslope in 5-m increments, to produce two larger zero-order basins, which can be used to constrain the uncertainties surrounding zero-order basin extraction from digital topography. We demonstrate below that zero-order basin size is relatively insensitive to changes of channel head location of the order of 10 m, which is the typical range of uncertainty in the field verification of channel heads identified by the DrEICH algorithm (Clubb et al., 2014), and in the literature more broadly (Julian et al., 2012).
In order to delineate the zero-order basin morphology from the base point defined by the DrEICH algorithm, we use an upslope flow accumulation method to identify all the cells of the DEM, which flow into the zero-order basin base ( Figure 1). Tests using flow routing algorithms such as multidirection dispersive methods (Freeman, 1991;Quinn et al., 1991) or the D∞ method (Tarboton, 1997) produced large zero-order basins, which often crossed between ridges, joining two patches, which would be field mapped as discrete basins as a single larger basin. Consequently, we employed a steepest descent upslope contributing area technique (O'Callaghan & Mark, 1984), which extracts topographically connected patches of the landscape. This algorithm is employed on each of the three confidence interval bases defined from the topography, which allows us to understand the influence that the initial starting point of our zero-order basins will have on the extracted morphologies and properties. Figure 3 demonstrates the spatial distribution of the middle bound zero-order basins in both of the physiographic units identified in Coweeta.
Once the zero-order basins are extracted for the three confidence intervals, their raster outlines are converted to vectors using GDAL (GDAL Development Team, 2013). Using these zero-order basin outlines, a number of topographic parameters can be extracted for each triplet of zero-order basins, allowing the analysis of landscape wide trends in addition to changes within triplets, which would indicate a strong sensitivity to channel head location in basin morphology.

Elliptical Fourier Analysis
To understand the 2-D geometry of zero-order basins, it is instructive to compute the average shape of a basin and to compare it to models of basin formation, which mainly focus on the evacuation of material and consequently assume a constant geometry (D'Odorico & Fagherazzi, 2003;D'Odorico et al., 2005;Parker et al., 2016) and theories of zero-order basin evolution, which argue for either basins that are stable or dynamic in space (Dietrich & Dorn, 1984). This can be achieved through elliptical Fourier analysis (EFA), which allows the analysis and comparison of 2-D shapes regardless of orientation, scale, coordinate density, spacing, or origin. This technique, initially presented by Kuhl and Giardina (1982), decomposes a closed contour as a series of ellipses, referred to as harmonics. The sum of these harmonics can be used to reconstruct the original closed contour, with increasing numbers of harmonics capturing increasingly fine variations in the shape of the original contour (Carlo et al., 2011).
This method has been applied in many disciplines where the 2-D morphology of objects can be used as a diagnostic tool, or to classify objects into predefined categories, including material science (Raj & Cannon, 1999), agricultural science (Costa et al., 2009(Costa et al., , 2011, biology (Carlo et al., 2011;Yoshioka et al., 2004), and paleontology (Crampton, 1995). Recently, this method has been demonstrated to have specific utility in the analysis of geomorphic objects such as watershed boundaries (Bonhomme et al., 2013), and a new software package has been developed with the aim of better applying EFA in the geosciences, which has been used in this analysis (Grieve, 2017).
Another advantage of EFA, compared to other shape description methods, such as the comparison of major and minor axes, is that the input data do not have to be evenly spaced, allowing an increase in coordinate density along complex sections of a shape and sparser sets of coordinates in simpler sections (Crampton, 1995). The use of EFA also allows the description of shapes at increasing levels of complexity. In some shapes high-frequency variations in the outline are important, whereas in others such variations are a product of the uncertainty of the digitization process. Consequently, EFA can be applied to both simple and complex shapes with users able to separate valuable high-frequency information from high-frequency noise (Crampton, 1995). This method is size invariant and is important for the analysis of zero-order basins as changes in area will constitute most of the variation in harmonic amplitudes and would consequently make up the majority of the statistical variance between basins (Crampton, 1995). For example, a perfectly circular basin with radius 5 m is identical to another circular basin with a radius of 100 m. This sets the EFA method apart from other shape description techniques as it facilitates the independent analysis of planform morphology and area.
The vector outlines generated in section 3.3.1 can be used for this analysis. The polyline of the perimeter of each zero-order basin was translated into an ordered series of Universal Transverse Mercator (UTM) x, y coordinates, representing each node on the perimeter polyline. The centroid of each basin is calculated as where and x and y are lists of UTM coordinates where (x 0 , y 0 ) = (x k , y k ) and k is the number of coordinate pairs in each zero-order basin outline. These coordinates are then normalized about this centroid to ensure no area bias in the analysis of the planform geometries is introduced, where x min and y min are the minimum x and y coordinates and x len and y len are the dimensions of the normalized bounding box in the x and y directions (Figure 4). This normalization has the additional benefit of lowering the computational burden of this technique when working with large data sets. These normalized coordinates are then rotated by an angle, , so that the outlet of each basin flows south, ensuring consistency between basins extracted from hillslopes of differing aspects, As demonstrated by Kuhl and Giardina (1982), for each harmonic of a contour, two Fourier coefficients can be generated for the x and y coordinates, resulting in four coefficients for each harmonic. The four coefficients, A n and C n , which describe the symmetry of a shape, and B n and D n , which describe the asymmetry of a shape, are calculated as where T is the period, equivalent to the length of the perimeter of the contour, t p is the total distance along the contour at p, x p is the distance along the x axis at point p, y p is the equivalent for the y axis. The complete derivation of these coefficients can be found in Kuhl and Giardina (1982) and many subsequent works, for example, Raj and Cannon (1999), and for the sake of brevity is not reproduced here.
To ensure that the set of coefficients for each zero-order basin can be compared, three separate normalizations must be applied (Kuhl & Giardina, 1982). The first removes the influence of the location of the origin of each contour from the coefficients The second stage requires normalization with respect to the rotation of individual contours to their location in space, ensuring that the major axes of all zero-order basins are aligned along the same plane.
As the zero-order basins have already been oriented to flow south, this will only make fine-scale adjustments to the coefficients, but both rotations must be performed to ensure all basins are aligned to the same major axis and flow to the south. The final operation normalizes the coefficients with regard to the absolute value of the coefficient A 1 Note that throughout the remainder of this paper the notation A * * * n will not be used, with the normalized coefficients being represented as A n for clarity.
In theory the total number of harmonics, which can be computed for a given shape, is equal to the Nyquist frequency (k∕2); however, in practice this will often yield an overfitted result with high-frequency noise accounting for much of the higher harmonic coefficient values. The Fourier power of a harmonic can be considered proportional to a measure of the amount of shape information provided by a given harmonic (Crampton, 1995). The Fourier power of a harmonic (P n ) can be calculated as P n = A 2 n + B 2 n + C 2 n + D 2 n 2 (17) and the total power (P T ) is given as Using the value of the total power, the cumulative power of increasing harmonics can be calculated until it reaches a desired fraction of the total power, at which point the series can be truncated and a limit on the number of harmonics has been identified. In this study we follow Crampton (1995) and select 0.9999 as the threshold beyond which further harmonics are not required, which for Coweeta corresponds to only working with the first 17 harmonics in our analysis.
To reconstruct a representation of a zero-order basin from a series of coefficients, an inverse Fourier transform can be applied (Kuhl & Giardina, 1982) x which yields a series of x and y coordinates representing the reconstructed contour generated using n harmonics.
To perform analysis of zero-order basins using this technique, it is useful to be able to average populations of outlines, to develop information about the typical geometry produced by a given landscape property. The average basin shape can be generated by averaging the normalized Fourier coefficients and solving equation (19) using the averaged values (Raj & Cannon, 1999). The average of a coefficient, C, is given as and the standard deviation, C , of an averaged zero-order basins is where N is the total number of basins being studied. This averaging process produces the average and standard deviation of a collection of coefficients, which can be transformed into average coordinates, to graphically represent the geometry of the average zero-order basin which represents a given landscape parameter or spatial location.

Zero-Order Basin Apex and Width Extraction
The apex of each extracted basin is identified by using a least cost algorithm to route a path between the highest elevation and lowest elevation point in each basin. The cost surface used for this routing is generated as |ln A| where A is the upslope contributing area, calculated using the D∞ algorithm (Tarboton, 1997). The D∞ algorithm is applied here rather than the steepest descent algorithm applied in section 3.3.1 as it provides the most hydrologically significant representation of surface flow accumulation on hillslopes (Shelef & Hilley, 2013). The aim is to identify the path of maximum overland flow accumulation through each basin, rather than simply the shortest flow path from top to bottom. Following the extraction of apex lines, a test is performed to ensure that the identified apex line falls within the original basin bounds. In some cases the least cost path will fail to generate a valid route through the basin, while in other cases the trace will exit the basin, rendering the measurement meaningless. Such cases are identified by comparing the length of the apex trace to the length of the trace, which falls within the basin outline, if these values are not similar, the basin is excluded from further analysis. An exact match is not required; rather, a threshold of a 90% match is used, in order to account for edge effects resulting from the conversion of raster outlines to vector data.
Zero-order basin width is extracted from the rotated vector outlines of basins produced in section 3.4 by projecting a line perpendicular to the hollow apex through the centroid of the basin, identified using equation (6), until it intersects the basin outline on both sides of the centroid. The distance between these two intersections is computed and given as the width of the basin. This method is preferred over dividing the basin area by the basin length to give an average width as it is analogous to the process of measuring zero-order basin width in the field.
We validated our estimates of zero-order basin shape against a data set of zero-order basin widths measured in the field as described by Parker et al. (2016). For each basin, we identified the hillslope noses as areas of convex upward topography. We then stretched a tape measure between the two noses, perpendicular to the hollow axis. This field-mapping campaign was supported by a curvature map derived from a 6-m resolution DEM. Candidate basins were identified on the basis of their concavity, with a minimum width of four concave pixels (or 24 m). Consequently, there are no basins below this threshold in this data set. Figure 5 presents the distribution of zero-order basin properties across Coweeta. In each of the subplots there is little variation between the three zero-order basin bounds, with each bound showing the same broad patterns with only limited small-scale variations. The basin area data exhibit an initial exponential decay, with the majority of basins only covering a small spatial area and a small number of outlying larger basins. Basin average topographic gradient exhibits a Gaussian-like distribution, with a median gradient of 26 ∘ and a maximum average gradient of 39 ∘ . The percentage of concave pixels in each basin shows a broadly normal distribution, with the majority of basins having approximately 50% concave pixels. Such a pattern conforms to the conceptual model of zero-order basins presented by Hack and Goodlett (1960), where colluvial hollows are bounded by planar side slopes and divergent noses (Figure 1), limiting the total percentage of concave pixels in each zero-order basin. The basin width data show a large proportion of narrow basins with a width of 10 m or less, and very few basins exceeding 100 m wide, which corresponds well with field observations (Parker et al., 2016). The length data show a high incidence of short basins, coinciding with the large number of basins with small areas but then decreases between approximately 40 and 90 m before increasing at 100 m and then declining in a similar manner to the area data. Even the smallest basins have some relief, so the probability of basins with little relief is low, rising to a maximum at intermediate relief values with a plateau of maximum probability between 10 and 70 m. There is then a long tail of basins with high relief, mirroring the presence of long zero-order basins.

Zero-Order Basin Properties
There is no clear difference between the three basin bounds we have extracted ( Figure 5), showing that the trends observed in morphology are not driven by the meter-scale placement of the base of a zero-order basin. Consequently, the data presented in the rest of the paper use the middle bound basin measurements.

Contrasts Between Zero-Order Basins in Different Physiographic Units
The two physiographic units present in Coweeta contain zero-order basins with distinct properties (Figure 6). The population of basins located on the escarpment have lower areas, and correspondingly lower widths and lengths, but span a wider range of basin average gradients and a larger number of concave pixels per basin than the population of basins located in the remainder of the Coweeta basin. The only property that does not vary significantly between the two populations is relief, where both data sets show similar ranges and median values.

Average Planform Shape of A Zero-Order Basin
The normalized average basin shape presented in Figure 7a demonstrates that the majority of zero-order basins are elliptical in nature, with little variation in shape ±1 standard deviation and the long axis aligned with downslope direction. Figures 7b and 7c show the average shape of basins on the Nantahala Escarpment and in the rest of Coweeta, and these demonstrate a fundamental difference in planform morphology between these two physiographic units, with the escarpment dominated by narrow basins and the remainder of Coweeta dominated by wider, more elliptical basins.
Segmenting the full population of zero-order basins by landscape properties shows similar patterns of differing morphology, which, when combined with the data presented in Figure 6 suggests that these patterns are merely functions of the broader physiographic patterns within the landscape. Figures 8a and 8b demonstrate that the majority of narrow basins have areas below the median basin area for the landscape (1,548 m 2 ), resulting in a much more uniform elliptical form for larger basins. Parker et al. (2016) demonstrated for a larger area of the Southern Appalachians, including Coweeta, the median colluvial hollow apex gradient is 28 ∘ . We use this value to segment the population of zero-order basins into shallow (< 28 ∘ ) and steep (≥ 28 ∘ ) categories, to explore the variation in morphology with gradient, which shows that steeper gradient basins are more narrow than shallower gradient basins (Figures 8c and 8d). There are many more east facing than west facing basins in Coweeta, and this aspect pattern results in distinct average shapes, with the majority of elliptical basins having an eastern aspect and narrower basins having a western aspect ( Figure 9); however, this pattern is impacted by the presence of the Nantahala Escarpment, which controls the aspect of soil-mantled hillslopes in much of this region.  Figure 10 gives a simple representation of the complexity of basin outlines for the whole landscape, with the majority of basins being represented to within 99.99% of their actual shape by between 14 and 21 harmonics. There are also a large number of basins, which can be adequately represented using fewer than three harmonics, suggesting that these basins are very elliptical in shape, with little small-scale variation in their outlines. If the basins are segmented by landscape properties such as aspect, gradient, curvature, or area, there is no significant change in the number of harmonics required to describe their outlines, suggesting that the complexity of a zero-order basin outline is not controlled by landscape property and, at least within the limitations of the data resolution and extraction method, basin outline complexity is broadly constant across Coweeta.

Evaluation of Extracted Zero-Order Basins
The widths of zero-order basins extracted from the DEM using the method described above can be evaluated against measurements taken by Parker et al. (2016) of 55 basins across the southern Appalachians. As only six of these mapped basins fall within the Coweeta basin, it is not feasible to perform a comparative spatial analysis on the locations of basins; the objective of the field-mapping campaign undertaken by Parker et al. (2016) was not to map every zero-order basin in the landscape, but rather to sample basins that could potentially generate shallow landslides. However, it is possible to perform a Monte Carlo analysis on the distribution of field measured and automatically extracted widths, to assess the similarity of the geometries extracted from high-resolution topographic data, to the extant field measurements.
To undertake this analysis, a subset of the full Coweeta zero-order basin data set was generated. Basins retained in the subset had a minimum width of 24 m, which was the minimum width of basins that were Figure 11. Probability density of the mean absolute error between field-mapped zero-order basin widths and those extracted from high-resolution topographic data. The probability density is generated using the Monte Carlo method to compare field data from Parker et al. (2016) with random selections of basins extracted from across Coweeta. The red dashed line indicates the median absolute error between these two measurement types.
field mapped. No other basins were excluded. From this subset, one million random subsamples of this processed data set were then generated, sorted into rank order and the absolute deviation between each random basin and a field measured basin was calculated. The results of this analysis are presented in Figure 11. In the vast majority of cases the absolute error between the automated and field measurement techniques fall below 5 m, comparable to the typical accuracy of handheld GPS devices (Julian et al., 2012). This suggests that the widths of the features extracted by the automated algorithm are similar to those that geomorphologists would identify in the field.

What Shape is A Zero-Order Basin?
The distributions of zero-order basin properties reported in Figure 5 demonstrate the range of topographic and morphometric properties observed in a population of basin in a steady state landscape. Zero-order basin width, length, and area all exhibit an initial exponential decay, with large numbers of small and narrow basins, and a smaller number of outlying larger, wider basins. The percentage of concave area in each basin and the topographic gradient are both gaussian-like and the basin relief is similar to a Wiebull distribution. These distributions provide a set of diagnostic zero-order basin parameters which can be used to test predictions of existing landscape evolution models and to develop tests for future models to better capture the spatial and geometric va riability within and between basins in a steady state landscape.
Observations of uniform valley spacing and the uniformity of first-order valleys (Perron et al., 2008(Perron et al., , 2009 have in the past lent support to the idea that a characteristic zero-order basin area or width should be observable in many steady state soil-mantled landscapes. However, the distribution of zero-order basin areas and widths presented in Figure 5 does not suggest that a characteristic basin geometry is extant in Coweeta. Perron et al. (2008) suggest that uniform valley spacing is driven by the interplay between advective and diffusive processes. The lack of similarity between zero-order basin areas across Coweeta suggests that in this landscape, differing rates of processes may be operating on the two physiographic units which can enhance the diversity of basin properties and morphologies observed at the landscape scale.
As highlighted in Figure 3, the spatial distribution of zero-order basins does not appear to form any kind of pattern, with basin morphology and geometry varying apparently randomly in space. This suggests that across the landscape the processes that form, fill, and empty zero-order basins are set by stochastic processes, independent of the properties measured in this study, and are modulated by the interplay of landscape, biological and climatic processes at the scale of individual zero-order basins (tens to hundreds of meters). These interpretations are supported by field observations, where pits dug into hollow apexes show a range of subsurface hydrologies, with some having significant amounts of flowing water while others remain dry to the bedrock interface. Similarly, the complex tectonics of the area lead to a range of fine scale geological structures that appear to affect individual basins.

Physiographic Controls on Zero-Order Morphology
The two physiographic units differ in that the escarpment has thinner soils and is at a higher elevation than the thicker soil-mantled remainder of the Coweeta basin. The forest type also differs between the two units, although the root cohesion of both units has been shown to be similar (Hales & Miniat, 2017). When the zero-order basin populations are segmented into these units, clear patterns emerge between the two sets of basins, as demonstrated by Figures 6 and 7: escarpment basins are characteristically narrower, more concave, and cover a smaller area than typical basins from the remainder of the Coweeta basin.
Landscape evolution has long been modeled, both numerically and conceptually, as the competition between advective and diffusive processes, which disturb and smooth the topographic surface, respectively . Perron et al. (2008) invoke a balance between advective fluvial incision and diffusive transport of hillslope material to account for the uniformity of first-order drainages in their study sites. However, given the 10.1029/2017JF004453 stochastic nature of zero-order basin erosion processes (Benda & Dunne, 1997;Wooten et al., 2008), it seems possible that areas with similar rates of diffusive sediment transport could produce a wide range of zero-order basin morphologies.
Diffusive sediment transport rates are fairly consistent within the Coweeta catchment ), yet we can demonstrate a diverse range of morphologies within zero-order basins. Potential advective processes that could drive the evolution of Coweeta zero-order basins include shallow landsliding and debris flow activity (Wooten et al., 2008) and possibly overland flow, which the authors have observed during intense convective storms. There is also a nonuniform cover of saprolite in zero-order basins, suggesting that different rates of chemical weathering may occur. There is no obvious spatial control on these processes. Similarly, there are no systematic differences in lithology (Hatcher, 1979) and root cohesion (Hales et al., 2009) across the basin.
Nevertheless, there are systematic differences in the shapes of zero-order basins across our two physiographic units highlighted by the variation in curvature distribution between hollows shown in Figure 12. Field observations suggest that there are corresponding differences in rates of sediment transport between the units, exemplified by thinner soils and more bedrock outcrops on the Nantahala Escarpment. However, the relationship between morphology and process is not simple. The distributions of gradients for the two populations of zero-order basins share similar mean values ( Figure 6), suggesting that these differences cannot be explained through slope-driven differences in the occurrence of linear versus nonlinear sediment transport (e.g., Roering et al., 1999). However, rates of diffusive transport may scale with the depth-slope product (Heimsath et al., 2005), which, for similar slopes (Figure 6b), suggests that differences in soil depth between the physiographic units could drive differences in the hillslope sediment transport rate, which may have implications for basin shape.
One possible explanation for the physiographical differences in basin morphology could be through the relationship between a depth-dependent diffusive transport and a debris flow-driven advective transport. In locations with thicker soils, such as the non-escarpment sections of the Coweeta basin, depth-integrated sediment flux is higher and basins accumulate material more rapidly, having the effect of stabilizing the zero-order basin by reducing concavity (Dietrich & Dunne, 1978;Sidle, 1984). Such features may reach a threshold soil thickness which precludes further evacuation, creating a population of basins which have reduced concavity and larger widths, as identified in Figure 6. In contrast, on the escarpment, where hillslope sediment flux may be lower due to a thinner soil mantle, basin concavities will remain high, potentially causing more frequent hollow evacuation than in the remainder of the basin. Because soil thickness and concavity modulate the propensity for debris flow occurrence (Dietrich et al., 1995), and corresponding basin erosion (Stock & Dietrich, 2006), it is possible that the narrowing of escarpment basins is driven by frequent debris flows incising preferentially in the hollow apex.

Contrasts Between Models of Zero-Order Basins and Topographic Measurements
The wide distribution of zero-order basin areas and morphologies demonstrated here highlight the challenges in attempting to model basin evolution and associated landslide hazard using a simple geometry. The filling and emptying trough model employed in many studies (e.g., D'Odorico & Fagherazzi, 2003;D'Odorico et al., 2005;Parker et al., 2016) is valuable for modeling landslide hazard but does not incorporate the wide spatial variability in zero-order basin geometries observed in our data. This indicates that it may be necessary to consider basin geometry in conjunction with models of sediment accumulation in order to better capture the signal of shallow landsliding in soil-mantled landscapes such as Coweeta.
Uniform valley spacing and a uniformity of drainage area can be observed in many landscape evolution models (e.g., Hobley et al., 2017;Tucker et al., 2001) and are predicted in many theories of landscape evolution (Perron et al., 2008), and drainage basins extracted from such models are typically elliptical in morphology.
When contrasted with measurements of zero-order basins extracted from high-resolution topography, the disparity between the two is apparent, with extracted basins exhibiting much more variability than their modeled counterparts. An explanation for this variability is the wide range of axis gradients which a population of zero-order basins will have, as demonstrated by Parker et al. (2016). Such a range is driven by a wide range of material properties such as soil and root cohesion and friction angle extant across a population of basins. We postulate that the high variability of zero-order basins observed at a fixed point in time from the high-resolution topographic data reflect the temporal stochasticity of landsliding across millennial timescales and the spatial heterogeneity of the processes which control the evolution of zero-order basins and their adjacent hillslopes.
However, this is not necessarily a criticism of the implementation of landscape evolution models and theories of landscape evolution in general but rather highlights the difference inherent in observing a complex natural system at a fixed point in time rather than a numerical model, which has evolved to steady state. The spatially averaged zero-order basin presented in Figure 7a demonstrates this, whereby averaging all of the basin morphologies from a fixed point in time presents an elliptical form which is more aligned with the uniform features identified in landscape evolution models. This suggests that models do indeed capture a significant proportion of the variability of natural systems, and although there is still much complexity to understand, these models are a valuable tool to explore the evolution of landscapes. A future development of such models will be to incorporate a stochastic advective process to evacuate modeled zero-order basins in competition with a diffusive filling process, with their respective rates set by a combination of topographic, biotic and climatic parameters unique to each basin.

Morphological Consistency of Extracted Zero-Order Basins
The technique employed here to extract zero-order basins from high-resolution topography incorporates three sets of channel heads to capture potential uncertainty in the precise location of channel heads across the landscape. Such features are well understood to be transient in nature with seasonal variability in their location on the meter scale observed in both field and experimental studies (Dietrich & Dunne, 1993). The comparison of zero-order basin morphologies extracted for the middle, upper, and lower bound channel heads presented in Figure 5 demonstrate the stability of the such features within the studied landscape, with no significant morphological changes identified based on movement of the channel head ±5 m. This value is selected as it corresponds to the maximum horizontal GPS error reported in the field verification of the DrEICH algorithm performed by Clubb et al. (2014).
The analysis of the number of harmonics required to represent 99.99% of the variation in a zero-order basin outline presented in Figure 10 demonstrates that the extraction of zero-order basin outlines from high-resolution topography is consistent across the landscape. This consistency in outline complexity reflects the ability of the zero-order basin extraction method to identify the signal of a zero-order basin from complex topographic data, without introducing meaningless noise to the perimeter measurements. These analyses provide confidence that the spatial and geometric properties of basins being studied here are not significantly impacted by the uncertainties inherent in the extraction methodology employed. Such stability of the features in Coweeta also suggests that it will be possible to perform similar analyses on other landscapes dominated by ridge-hollow terrain.

Conclusions
We present a technique to extract the outlines of zero-order basins from high-resolution topography, by identifying the upslope extent of the channel network and extracting the zero-order drainage above that point. Using this technique, a data set of over 1,000 zero-order basins from across Coweeta was created and used to understand the variations in spatial location and landscape properties, which exist across this data set. Diagnostic distributions of basin parameters were identified, which indicate the spatial variability in zero-order basin geometry and call into question the value of the use of a single characteristic length scale as a descriptive metric for zero-order basins. These features have complex forms, and only by studying them in their entirety will we obtain an understanding of the processes which govern their topographic development.
Two physiographic units were identified within the Coweeta basin, with distinct soil thicknesses and average gradients and the extracted basins were divided into two populations from these two units to further explore the controls on zero-order basin distribution and morphology.
The application of EFA allowed the quantitative analysis of zero-order basin planform morphology, which provided an insight into the elliptical nature of the landscape average basin, which corresponds well to predictions made by conceptual and numerical models of landscape evolution. Zero-order basins located on the Nantahala Escarpment were shown to be more concave and narrower than those located in the remainder of the basin, with characteristic planform morphologies for each physiographic unit identified.
Finally, a conceptual framework was presented, highlighting the competition between advective landsliding and diffusive sediment transport as the mechanism driving variability in zero-order basin morphology both