Pattern and structure of microtopography implies autogenic origins in forested wetlands

Pattern and structure of microtopography implies autogenic origins in forested wetlands Jacob S. Diamond, Daniel L. McLaughlin, Robert A. Slesak, and Atticus Stovall Quantitative Ecohydrology Laboratory, RiverLy, Irstea, Lyon, 69100, France Continental Hydrogeosytems Laboratory, University of Tours, Tours 37200, France 5 School of Forest Resources and Environmental Conservation, Virginia Tech, Blacksburg, 24060, USA Minnesota Forest Resources Council, St. Paul, 55108, USA NASA Goddard Space Flight Center, Greenbelt, 20771, USA


Introduction
Microtopography, or the small-scale structured variation (10 -1 -10 0 m) in ground surface height, is common to many ecosystems. Wetland microtopography is particularly well studied, and is found in freshwater marshes (Van de Koppel et al., 2006), fens (Sullivan et al., 2008), peat bogs (Nungesser, 2003), forested swamps (Bledsoe and Shear, 2000), tidal freshwater swamps (Duberstein et al., 2013), and coastal marshes (Stribling et al., 2007). Wetland 35 microtopography is common enough that researchers in disparate systems collectively refer to local high points as In this study, we evaluated wetland soil elevations, hummock spacing, and hummock sizes and their associations with hydrologic regimes in black ash (Fraxinus nigra Marshall) forested wetlands in northern Minnesota, U.S.A. To do so, we characterized microtopography with a 1-cm spatial resolution dataset from a terrestrial laser scanning (TLS) 95 campaign. We also evaluated subsurface mineral layer topography and daily water levels to determine the extent that that these variables influenced observed surface microtopography. Specifically, we tested the following predictions: 1. Elevation will exhibit a bimodal distribution, but the degree of bimodality and the overall variability in elevation will be greater in wetter sites than drier sites.
2. Surface topography will not reflect subsurface mineral topography, but will instead be representative of self-100 organizing processes at the soil surface.
3. Hummock heights will be positively correlated with water levels at site and within-site scales. 4. Hummock patches will exhibit spatial overdispersion, which will be more evident at wetter sites.
5. Cumulative distributions of hummock areas (and perimeters and volumes) will correspond to a family of truncated distributions (e.g., exponential or lognormal), indicating a characteristic patch size, with wetter 105 sites exhibiting more large area hummocks than drier sites.

Site descriptions
To test our hypotheses, we investigated ten black ash wetlands of varying size and hydrogeomorphic landscape position in northern Minnesota, U.S.A. (Figure 2; Table 1). Thousands of meters of sedimentary rocks overlay an 110 Archean granite bedrock geology in this region. Study sites are located on a glacial moraine landscape (400-430 m ASL) that is flat to gently rolling, with the black ash wetlands found in lower landscape positions that commonly grade into aspen or pine-dominated upland forests. The climate is continental, with mean annual precipitation of 700 mm and a mean growing season (May-October) temperature of 14.3°C (mean annual temperature = -1.1°C -4.8°C; WRCC 2019). Annual precipitation is approximately two-thirds rain and one-third snowfall. Potential 115 evapotranspiration (PET) is approximately 600-650 mm per year (Sebestyen et al., 2011). Detailed site histories were unavailable for the ten study wetlands, but silvicultural practices in black ash wetlands have been historically limited in extent (D'Amato et al., 2018). Based on the available information (e.g., Erdmann et al., 1987;Kurmis and Kim, 1989), we surmise that our sites are late successional or climax communities and have not been harvested for at least a century. 120 As part of a larger effort to understand and characterize black ash wetlands (D'Amato et al., 2018), we categorized and grouped each wetland by its hydrogeomorphic characteristics as follows: 1) depression sites ("D", n = 4) characterized by a convex, pool-type geometry with geographical isolation from other surface water bodies and surrounded by uplands, 2) lowland sites ("L", n = 3) characterized by extensive wetland complexes on flat, gently sloping topography, and 3) transition sites ("T", n = 3) characterized as flat, linear boundaries between uplands and 125 black spruce (Picea mariana Mill. Britton) bogs ( Figure 3). The three lowland sites were control plots from a longterm experimental randomized block design on black ash wetlands (blocks 1, 3, and 6;Slesak et al., 2014;Diamond et al., 2018). We considered hydrogeomorphic variability among sites an important criterion, as it allowed us to capture expected differences in hydrologic regime and thus differences in the strength of our predicted control on microtopographic generation (Figure 1). Ground slopes across sites ranged from 0-1%. Black ash wetlands are 130 typically hydrologically disconnected from regional groundwater and other surface water bodies, resulting in precipitation and evapotranspiration (ET) as dominant components of the water budget, with no indication of extreme surface flows (Slesak et al., 2014). Water levels follow a common annual trajectory of late-spring/early-summer inundation (10-50 cm) followed by ET-induced summer drawdown and belowground water levels (Slesak et al., 2014;Diamond et al., 2018). However, the degree of drawdown depends on local hydrogeomorphic setting; we observed 135 considerably wetter conditions at depression and transition sites than at lowland sites.

Vegetation
Overstory vegetation at the ten sites is dominated by black ash, with tree densities ranging from 650 stems ha -1 (basal area = 195 m 2 ha -1 ) at the driest lowland site to 1600 stems ha -1 (basal area = 40 m 2 ha -1 ) at a much wetter depression site (across-site mean = 942 stems ha -1 ; Diamond et al. 2019). At the lowland sites, other overstory species were 140 negligible, but at the depression and transition sites there were minor cohorts of northern white-cedar (Thuja occidentalis L.), green ash (Fraxinus pennsylvanica Marshall), red maple (Acer rubrum L.), yellow birch (Betula alleghaniensis Britt.), balsam poplar (Populus balsamifera L.), and black spruce (Picea mariana Mill. Britton). Except at one transition site (T1), where northern white cedar represented a significant overstory component, black ash represented over 75% of overstory cover across all sites. Black ash also made up the dominant midstory component 145 in each site, but was regularly found with balsam fir (Abies balsamea L. Mill.) and speckled alder (Alnus incana L. Moench) in minor components, and greater abundances of American elm (Ulmus Americana L.) at lowland sites.
Black ash stands are commonly highly uneven-aged (Erdmann et al., 1987), with canopy tree ages ranging from 130-232 years, and stand development under a gap-scale disturbance regime (D'Amato et al., 2018). Black ash are also typically slow-growing, achieving heights of only 10-15 m and diameters at breast height of only 25-30 cm after 100 150 years (Erdmann et al., 1987). The relatively open canopies of black ash wetlands (leaf area index < 2.5; Telander et al., 2015) allow for a variety of graminoids, shrubs, and mosses to grow in the understory. However, the majority of understory diversity and biomass tends to occur on hummocks that are occupied by black ash trees (Diamond et al., 2019). Hollows exhibit relatively little plant cover and are typically bare soil areas, but may be covered at times of the year by sedges (Carex spp.) or layers of duckweed (Lemna minor L.), especially after recent inundation. 155

Soils
Soils in black ash wetlands in this region tend to be Histosols characterized by deep mucky peats underlain by silty clay mineral horizons, although there were clear differences among site groups (NRCS 2019). Depression sites were commonly associated with Terric haplosaprists of the poorly drained Cathro or Rifle series with O horizons approximately 30-150 cm deep (Table 1). Lowland sites were associated with lowland Histic inceptisols of the 160 Wildwood series, which consist of deep, poorly drained mineral soils with a thin O horizon (< 10 cm) underlain by clayey till or glacial lacustrine sediments. Transition sites typically had the deepest O horizons (> 100 cm), and were associated with typic haplosaprists of the Seelyeville series and Typic haplohemists (NRCS 2019). Both depression and transition sites had much deeper O horizons than lowland sites, but depression site organic soils were typically muckier and more decomposed than more peat-like transition site soils. 165

Data collection
To characterize the microtopography of our sites, we conducted a terrestrial laser scanning (TLS) campaign from October 20-24, 2017. We chose this period to ensure high-quality TLS acquisitions, as it coincided with the time of least vegetative cover and the least likelihood for inundated conditions. During scanning, leaves from all deciduous 170 canopy trees had fallen and grasses had largely senesced. Standing water was present at portions of three of the sites and was typically dispersed across the site in small pools (ca. 0.5-2 m 2 ) less than 10 cm deep. We used a Faro Focus 8 Stovall et al., 2019 for exact methodological details). For each site, we merged our plot-level TLS data to a single ~900 m 2 site-level point-cloud using 30 strategically placed and scanned 7.62 cm radius polystyrene registration 175 spheres set atop 1.2 m stakes. We referenced each site to a datum located at each site's base well elevation (see section 2.3.1).
To validate the TLS surface model products, we installed sixty 2.54 cm radius spheres on fiberglass stakes exactly 1.2 m above ground surface at each site. With the validation locations we could easily calculate the exact surface elevation (i.e., 1.2 m below a scanned sphere) of 60 points in space. We installed 39 (13 at each plot) validation spheres at points 180 according to a random walk sampling design, and placed 21 (7 at each plot) validation spheres on distinctive hummock-hollow transitions. We placed the 1.2 m tall validation spheres approximately plumb to reduce errors due to horizontal misalignment.
We processed the point clouds generated from the TLS sampling campaign to generate two products: 1) site-level 1cm resolution ground surface models, and 2) site-level delineations of hummocks and hollows. The details and 185 validation of this method are described completely in Stovall et al. (2019), but a brief summary is provided here.

Surface model processing and validation
For each site, we first filtered the site-level point-clouds in the CloudCompare software (Othmani et al., 2011) and created an initial surface model with the absolute minima in a moving 0.5 cm grid. We removed tree trunks from this initial surface model using a slope analysis and implemented a final outlier removal filter to ensure all points above 190 ground level were excluded. Our final site-level surface models meshed the remaining slope-filtered point cloud using a local minima approach at 1 cm resolution. We validated this final 1-cm surface model using the 60 validation spheres per site.
Before we analyzed surface models from each site, we first detrended sites that exhibited site-scale elevation gradients (e.g., 0.02 cm m -1 ). These gradients may obscure analysis of site-level relative elevation distributions (Planchon et al., 195 2002), and our hypothesis relates to relative elevations of hummocks and hollows and not their absolute elevations.
We chose the best-detrended surface model based on adjusted R 2 values and observation of resultant residuals and elevation distributions from three options: no detrend, linear detrend, and quadratic detrend. Five sites were detrended: L2 was detrended with a linear model, and D1, D2, D4, and T1 were detrended with quadratic models. We then subsampled each surface model to 10,000 points to speed up processing time as original surface models were 200 approximately 100,000,000 points. We observed no significant difference in results from the original surface model based on our subsampling routine.

Hummock delineation and validation
We classified the final surface model into two elevation categories: hummocks and hollows. We first classified hollows using a combination of normalized elevation and slope thresholds; hollows have less than average elevation 205 and less than average slope. This combined elevation and slope approach avoided confounding hollows with the tops of hummocks since the tops of hummocks are typically flat or shallow sloped. We removed hollows and used the remaining area as our domain of potential hummocks.
Within the potential hummock domain, we segmented hummocks into individual features using a novel approach -TopoSeg (Stovall et al., 2019) and thereby created a hummock-level surface model for each site. We first used the 210 local maximum (Roussel and Auty, 2018) of a moving window to identify potential microtopographic structures for segmentation. The local maximum served as the "seed point" from which we then applied a modified watershed delineation approach (Pau et al., 2010). The watershed delineation inverts convex topographic features and finds the edge of the "watershed", which in our case are hummock edges. The defined boundary was used to clip and segment hummock features into individual hummock surface models. 215 For each delineated hummock within all sites, we calculated perimeter length, total area, volume, and height distributions relative both to local hollow datum and to a site level datum. To calculate area, we summed total number of points in each hummock raster multiplied by the model resolution (1 cm 2 ). We calculated volume using the same method as area, but multiplied by each points' height above the hollow surface. Perimeter was conservatively estimated by converting our raster-based hummock features into polygons and extracting the edge length from each 220 hummock. We estimated lateral hummock area by modelling each hummock as a simple cone, and calculating the lateral surface area from previously estimated volume and height. We believe this conical estimation method to be a conservative representation of the average height around the perimeter of the hummock because real hummock shapes are more undulating and complex than simple cones. We elected not to use a cylindrical model because we observed some tapering of hummocks from their base to their top. We note that a cylindrical model would increase lateral 225 surface area estimation by approximately 15% compared to the conical model and therefore may provide an upper bound on our conservative estimates.
To validate the hummock delineation, we compared manually delineated and automatically delineated hummock size distributions at one depression site (D2) and one transition site (T1), both with clearly defined hummock features. We omitted using a lowland site for validation because none of these sites had obvious hummock features that we could 230 manually delineate with confidence. We manually delineated hummocks for the D2 and T1 sites with a qualitative visual analysis of raw TLS scans using the clipping tool in CloudCompare (2018). Stovall et al. (2019) found no significant differences between the manual and automatically segmented hummock distributions, and feature geometry had an RMSE of less than ~20%.
After the automatic delineation procedure and subsequent validation, we performed a data cleaning procedure by 235 manually inspecting outputs in the CloudCompare software. We eliminated clear hummock mischaracterization that was especially prevalent at the edges of sites, where point densities were low. We also excluded downed woody debris from further hummock analysis because, although these features may serve as nucleation points for future hummocks, they are not traditionally considered hummocks and their distribution does not relate to our broad hypotheses. Finally, we excluded delineated hummocks that were less than 0.1 m 2 in area because we did not observe hummocks less than 240 this size during our field visits. This delineation and manual cleaning process yielded point clouds of hummocks and hollows for every site that could be further analyzed.

Surface model performance
Validation of surface models using the validation spheres indicated that surface models were precise (RMSE = 3.67 ± 1 cm) and accurate (bias = 1.26 ± 0.1 cm) across all sites (Stovall et al. 2019). The gently sloping lowland sites (L) 245 had substantially higher RMSE and bias than the transition (T) and depression (D) sites. The relatively high error of lowland site validation points resulted from either low point density or a complete absence of LiDAR returns. We observed overestimation of the surface model when TLS scans were unable to reach the ground surface, leading to the greatest overestimations in sites with dense grass cover (lowland sites). Overestimation was also common in locations with no LiDAR returns, such as small hollows, where the scanner's oblique view angle was unable to reach. 250 Nonetheless, examination of the surface models indicated clear ability of the TLS to capture surface microtopography ( Figure S1).

Hummock delineation performance
Hummocks delineated from our algorithm were generally consistent in distribution and dimension with manually delineated hummocks. However, the automatic delineation located hundreds of small (<0.1 m 2 ) "hummock" features 255 that were not captured with manual delineation, which we attribute to our detrending procedure. We did not consider automatically delineated hummocks less than 0.1 m 2 in further analyses, as we did not observe hummocks smaller than this in the field. Both area and volume size distributions from the manual and automatic delineations were statistically indistinguishable for both t-test (p-value = 0.84 and 0.51, respectively) and Kolmogorov-Smirnov test (pvalue = 0.40 and 0.88, respectively). Automatically delineated hummock area, perimeter:area, and volume estimates 260 had 23%, 19.6%, and 24.1% RMSE, respectively, and the estimates were either unbiased or slightly negatively biased (-9.8 %, 0.2 %, and -11.9 %, respectively). We consider these errors to be well within the range of plausibility, especially considering the uncertainty involved in manual delineation of hummocks, both in the field and on the computer. Final delineations showed clear visual differences among site types in the spatial distributions of hummocks ( Figure S2). 265

2.3
Field data collection

Hydrology
To address our hypothesis that hydrology is a controlling variable of microtopographic expression in black ash wetlands, we instrumented all 10 sites to continuously monitor water level dynamics and precipitation. Three sites (L1, L2, and L3; Slesak et al., 2014) were instrumented in 2011 and seven in June 2016 following the same protocols. 270 At each site, we placed a fully-slotted observation well (schedule 40 PVC, 2-inch diameter, 0.010-inch-wide slots) at approximately the lowest elevation; at the flatter L sites, wells were placed at the approximate geographic center of each site. Ground surface at the well served as each site's datum (i.e., elevation = 0 m). We instrumented each well with a high-resolution total pressure transducer (HOBO U20L-04, resolution = 0.14 cm, average error = 0.4 cm) to record water level time series at 15-minute intervals. We dug each well with a hand auger to a depth associated with 275 the local clay mineral layer and did not penetrate the mineral layer, which ranged from 30 cm below the soil surface to depths greater than 200 cm. We then backfilled each well with a clean, fine sand (20-40 grade). At each site, we also placed a dry well with the same pressure transducer model to measure temperature-buffered barometric pressure and frequency for barometric pressure compensation (McLaughlin and Cohen, 2011).

Mineral layer depth measurements 280
To quantify the control that underlying mineral layer microtopography has on surface microtopography, we conducted synoptic measurements of mineral layer depth and thus organic soil thickness at each site. Within each of the 10 m diameter plots used for TLS at each site, we took 13 measurements (co-located with the randomly established validation spheres) of depth to mineral layer using a steel 1.2 m rod. At each point the steel rod was gently pushed into the soil with consistent pressure until resistance was met and the depth to resistance was recorded (resolution = 1 285 cm) as the depth to mineral layer. We then associated each of these depth-to-mineral-layer measurements with a soil elevation based on TLS data and the site-level datum (i.e., elevation at the base of each site's well).

2.4
Data analysis

Hydrology
We calculated simple hydrologic metrics based on the three years (2016)(2017)(2018) of water level data for each site. For 290 each site, we calculated the mean and variance of water level elevation relative to ground surface at the well, where negative values represent belowground water levels and positive values indicate inundation. We also calculated the average hydroperiod of each site by counting the number of days that the mean daily water level was above the soil surface at the well each year, and averaging across years.

Elevation distributions 295
Our first line of inquiry was to evaluate the general spatial distribution of elevation at each site. We first calculated site-level omni-directional and directional (0°, 45°, 90°, 135°) semivariograms using the gstat package in R (Pebesma 2004 andGräler, 2016). We calculated directional variograms to test for effects of anisotropy (directional dependence) of elevation. Semivariogram analysis is regularly used in spatial ecology to determine spatial correlation between measurements (Ettema and Wardle, 2002). The sill, which is the horizontal asymptote of the semivariogram, is 300 approximately the total variance in parameter measurements. The nugget is the semivariogram y-intercept, and it represents the parameter variance due to sampling error or the inability of sampling resolution to capture parameter variance at small scales. The larger the difference between the sill and the nugget (the "partial sill"), the more spatially predictable the parameter. If the semivariogram is entirely represented by the nugget (i.e., slope = 0), the parameter is randomly spatially distributed. The semivariogram range is the distance where the semivariogram reaches its sill, and 305 it represents the spatial extent (patch size) of heterogeneity, beyond which data are randomly distributed. When spatial dependence is present, semivariance will be low at short distances, increase for intermediate distances, and reach its sill when data are separated by large distances. We used detrended elevation models for this analysis to assess more directly the importance of microtopography on elevation variation as opposed to having it obscured by site-level elevation gradients. From these semivariograms we calculated the best-fit semivariogram model among exponential, 310 Matérn, or Matérn with Stein parameterization model forms (Minasny and McBratney, 2005). We also extracted semivariogram nuggets, ranges, sills, and partial sills.
Our second line of inquiry was to evaluate the degree of elevation bimodality in these systems, which is indicative of a positive feedback between hummock growth and hummock height (Eppinga et al., 2008). Based on the classification into hummock or hollow from our delineation algorithm, we plotted site-level detrended elevation distributions for 315 hummocks and hollows and determined a best-fit Gaussian mixture model with Bayesian Information Criteria (BIC) using the mclust package (Scrucca et al., 2016) in R (R Core Team, 2018), which uses an expectation-maximization algorithm. Mixture models were allowed to have either equal or unequal variance, and were constrained to a comparison of bimodal versus a unimodal mixture distribution.

Subsurface topographic control on microtopography 320
We assessed the importance of mineral layer microtopography on soil surface microtopography by comparing the depth-to-mineral-layer measurements with the soil surface elevation TLS measurements. We first calculated the elevation of the mineral layer relative to each site-level datum by subtracting the depth-to-mineral-layer measurement from its co-located soil elevation measurement estimated from the TLS campaign. We then plotted the depth-tomineral-layer measurement (hereafter referred to as "organic soil thickness") as a function of this mineral layer 325 elevation, noting which points were on hummocks or hollows as determined from the TLS delineation algorithm. We fit linear models to these points and compared the regression slopes to the expected slopes from: 1) a scenario where surface microtopography is simply a reflection of subsurface microtopography (slope = 0, or constant organic soil thickness), and 2) a scenario of flat soil surface where organic soil thickness negatively corresponds to varying mineral layer elevation (slope = -1, or varying soil thickness). The first scenario would indicate that surface microtopography 330 mimics subsurface microtopography, whereas the second would indicate organic matter/surface soil accumulation and smoothing over a varying subsurface topography. Observations above the -1:1 line would indicate surface processes that increase elevation above expectations for a flat surface.

Hydrologic controls on hummock height
To test our hypothesis that hydrology is a broad, site-level control on hummock height, we first regressed site mean 335 hummock height against site mean daily water level. We also conducted a within-site regression of individual hummock heights against their local mean daily water level. To do so, we first calculated a local relative mean water level for each delineated hummock location by subtracting the elevation minimum of the hummock (i.e., the elevation at the base of the hummock) from the site-level mean water level elevation. This calculation assumes that the water level is flat across the site, which is likely valid for the high permeability organic soils at each site, low slopes (<1%), 340 and relatively small areas that we assessed. This within-site regression allowed us to understand more local-scale controls on hummock height.

Hummock spatial distributions
To test whether there was regular spatial patterning of hummocks at each site, we compared the observed distribution of hummocks against a theoretical distribution of hummocks subject to complete spatial randomness (CSR) with the 345 R package spatstat (Baddeley et al., 2015). We first extracted the centroids and areas of the hummocks using TopoSeg (Stovall et al., 2019) and created a marked point pattern of the data. Using this point pattern, we conducted a nearestneighbor analysis (Diggle, 2002), which evaluates the degree of dispersion in a spatial point process (i.e., how far apart on average hummocks are from each other). If hummocks are on average further apart (using the mean nearest neighbor distance, μNN) compared to what would be expected under CSR (μexp), the hummocks are said to be 350 overdispersed and subject to regular spacing; if hummocks are closer together than what CSR predicts, they are said to be underdispersed and subject to clustering. We compared the ratio of μNN and μexp, where values greater than 1 indicate overdispersion and values below 1 indicate clustering, and calculated a z-score (zANN) and subsequent p-value to evaluate the significance of overdispersion or clustering (Diggle, 2002, Watts et al., 2014. Z-scores were computed from the difference between μNN and μexp scaled by the standard error. We also evaluated the probability distribution 355 of observed nearest neighbor distances to visualize further the dispersion of wetlands in the landscape.

Hummock size distributions
To test the prediction that hummock sizes are constrained by patch-scale negative feedbacks, we plotted site-level rank-frequency curves (inverse cumulative distribution functions) for hummock perimeter, area, and volume. These curves trace the cumulative probability of a hummock dimension (perimeter, area, or volume) being greater than or equal to a certain value (P[X≥x]). We then compared best-fit power (P[X≥x] = αX β ), log-normal (P[X≥x] = βln(X)+ β0), and exponential (P[X≥x] = αe βX ) distributions for these curves using AIC values. Power-scaling of these curves occurs where negative feedbacks to hummock size are controlled at the landscape-scale (i.e., hummocks have approximately equal probability to be found at all size classes). Truncated scaling of these curves, as in the case of exponential or lognormal distributions, occurs when negative feedbacks to hummock size are controlled at the patch-365 scale (Scanlon et al., 2007, Watts et al., 2014.

Hydrology
Hydrology varied across sites, but largely corresponded to hydrogeomorphic categories (Table 2). Depressions sites were the wettest sites (mean daily water level = -0.01 m), followed by transition sites (-0.04 m), and lowland sites (-370 0.32 m). Lowland sites also exhibited significantly more water level variability than transition or depression sites, whose water levels were consistently within 0.4 m of the soil surface. Although lowland sites exhibited greater water level drawdown during the growing season, they were able to rapidly rise after rain events.

Elevation distributions
Semivariograms demonstrated much more pronounced elevation variability at depression and transition sites than at 375 lowland sites (Figure 4). In general, lowland sites reached overall site elevation variance (sills, horizontal dashed lines) within 5 meters, but best-fit ranges (dotted vertical lines in Figure 4) were less than 1 m. In contrast, best-fit semivariogram ranges for depression and transition sites were several times greater. Therefore, depression and transitions sites have much larger ranges of spatial autocorrelation for elevation than lowland sites. Semivariograms were all best fit with Matérn models with Stein parameterizations, and nugget effects were extremely small in all cases 380 (average <0.001), which we attribute to the very high precision of the TLS method. As such, partial sills were quite large (i.e., the difference between the sill and nugget), indicating that very little elevation variation is at scales less than our surface model resolution (1 cm); the remaining variation is found over site-level ranges of autocorrelation.
We did not observe major differences in directional semivariograms compared to the omnidirectional semivariogram, implying isotropic variability in elevation, and do not present them here. 385 We observed bimodal elevation distributions at every site, with hummocks clearly belonging to a distinct elevation class separate from hollows ( Figure 5). Bimodal mixture models of two normal distributions were always better fit to the data than unimodal models based on BIC values. Differences in mean elevations between these two classes ranged from 12 cm at the lowland sites to 20 cm at depression sites, and hummock elevations were more variable than hollow elevations across sites. Across sites, 27±10% of all elevations did not fall into either a hummock or a hollow category, 390 with lowland sites having considerably more elevations not in these binary categories (36-44%) compared to depression (22-27%) or transition sites (16-22%). However, we emphasize that even when considering the entire site elevation distribution (i.e., including elevations that did not fall into a hummock or hollow category), bimodal fits were still better than unimodal fits, but to a lesser extent for lowland sites ( Figure S3). Delineated hummocks varied in number and size across and within sites. We observed the greatest number of hummocks in the depression and 395 transition sites, with approximately an order of magnitude less hummocks found in lowland sites ( Figure 5).

Subsurface topographic control on microtopography
Across sites, organic soil thickness varied and was greatest at the lowest mineral layer elevations, indicating that surface microtopography is not simply a reflection of subsurface mineral layer topography with constant overlying organic thickness (as illustrated with 0-slope line in Figure 6). In contrast, at most sites, except for D1 and L2, there 400 was a strong negative linear relationship between soil thickness and mineral layer elevation, with five sites exhibiting slopes near -1, which we define as the smooth surface model of soil elevation (dashed -1:1 line in Figure 6). If only hollows (open circles; Figure 6) were used in the regression, then D1 also exhibited a significant (p<0.001) negative slope in this relationship (-0.4, R 2 = 0.52). A majority of depth to mineral layer measurements at D3 were below detection limit with our 1.2 m steel rod, and all but one measurement at T1 were below detection limit. At sites D2 405 and L2, there was indication that some hollows were actually better represented by the subsurface reflection model (i.e., slope = 0). However, at all sites, though to a lesser extent at lowland sites (e.g., L1 and L3), hummocks (closed circles; Figure 6) tend to plot above hollows and above the -1:1 line, indicating that their elevation is greater than would be expected for a smooth surface model.

Hydrologic control on hummock height 410
We observed a significant (p<0.001) positive linear relationship between site level mean hummock height and site level mean daily water level (Figure 7, top panel). Because lowland sites were clearly influential points on this linear relationship, we also conducted this regression excluding the lowland sites and still found a significant (p = 0.007) positive linear trend between these variables with reasonable predictive power (R 2 =0.8)wetter sites have on average have taller hummocks than drier sites. We found very little variability in average hummock heights across sites when 415 relative to site-level mean water level elevation (mean normalized hummock height = 0.31±0.06 m), indicating that hummocks were generally about 30 cm higher than the site mean water level.
Within sites, we also observed clear positive relationships between individual hummock heights and their local mean daily water level (Figure 7, bottom panel). At all but two of the sites (D4 and L1), individual hummock heights within a site were significantly (p<0.01) taller at wetter locations than drier locations. Slopes for these individual hummock 420 regressions varied among sites, ranging from 0.4-1.1 (mean±sd = 0.7±0.2), and local hummock mean water level was able to explain 12-56% (mean±sd = 0.36±0.14) of variability in hummock height within a site.

Hummock spatial distributions
All sites characterized as depressions or transitions exhibited significant (p <0.001) overdispersion of hummocks compared to what would be predicted under complete spatial randomness ( Figure 8). For these sites, the nearest 425 neighbor ratios (μNN:μexp) indicated that hummocks are 25-30% further apart than would be expected with complete spatial randomness, with spacing ca. 1.5 meters, as evidenced by the narrow distributions in nearest neighbor histograms (Figure 8). In contrast, all lowland sites, while having hummock nearest neighbor distances 2-3 times as far apart as depression of transition sites, were not significantly different than what would be predicted under complete spatial randomness (p = 0.129, 0.125, 0.04 for sites L1, L2, and L3, respectively). 430

Hummock size distributions
Hummock dimensions (perimeter, area, and volume) were strongly lognormally distributed across sites (Figure 9), though exponential models were typically only slightly worse fits. For each hummock dimension, site fits were similar within site hydrogeomorphic categories, but drier lowland site distributions were clearly different from wetter depression and transition site distributions, which were more similar (Figure 9). Lowland sites had significantly lower 435 (p < 0.05) coefficients for hummock property model fits than depression or transition sites, with slopes approximately 20% more negative on average, indicating more rapid truncation of size distributions. Across sites, average hummock perimeter was 4.2±0.8 m, average hummock area was 1.7±0.5 m 2 , and average hummock volume was 0.17±0.06 m 3 .
Hummock areas were typically less than 1 m 2 in size at all sites (Figure 9). Similar to hummock spatial density, hummock area per site (ratio of hummock area to site area) was lower at drier lowland sites (2-5%) compared to Discussion We tested our hypothesis that microtopography in black ash wetlands self-organizes in response to hydrologic drivers ( Figure 1) using an array of commonly used diagnostic tests from landscape ecology, including analyses of multimodal elevation distributions, spatial patterning, and patch size distributions. We further analyzed the influence of hydrology 445 on these diagnostic measures and tested a potential null hypothesis that surface microtopography was simply a reflection of subsurface microtopography. Diagnostic test results of elevation bimodality, hummock spatial overdispersion, and truncated hummock areas along with clear hydrologic influence on microtopographic structure provide strong support for our hypothesis.

Controls on microtopographic structure 450
Bimodal soil elevation distributions at all sites suggest that the microsite separation into hummocks and hollows is a common attribute of black ash wetlands. Soil elevation bimodality was most evident at the wetter depression and transition sites, where hummocks were more numerous and occupied a higher fraction of overall site area (15-20%).
Sharp boundaries between hummocks and hollows were not always observed in soil elevation probability densities ( Figure 5), which may be indicative of weak positive feedbacks between primary productivity and elevation (Rietkerk 455 et al., 2004; Figure 1). On the other hand, modeling predictions indicate that if evapoconcentration feedbacks (i.e., that hummocks harvest nutrients from hollows through hydraulic gradients driven by hummock-hollow ET differences) are strong, boundaries between hummocks and hollows will be less sharp (Eppinga et al., 2009), possibly implicating hummock evapoconcentration as an additional feedback to hummock maintenance (Figure 1). Greater levels of soil chloride in hummocks relative to hollows in these systems may be an additional layer of evidence for 460 this mechanism (Diamond et al., 2019).
We also observed clear evidence of decoupling between surface microtopography and mineral layer microtopography at all of our sites. Hollows were best represented by a smooth surface model, with a relatively constant surface elevation despite variable underlying mineral soil elevation. Importantly, we also observed that regardless of underlying mineral layer, hummocks had greater soil thickness than hollows ( Figure 6). That is, irrespective of mineral 465 layer microtopography, hummocks are maintained at local elevations that are higher than would be predicted for a smooth soil surface. Moreover, drier lowland (L) sites had less clear patterns in this regard than the wetter depression (D) or transition (T) sites, supporting our hypothesis for hydrology driven hummock development. We also note that some measurement locations had deeper organic soils than we could measure with our rod (particularly at our wettest sites) and that this is likely further evidence for our contention that hummocks are self-organized mounds on a smooth 470 surface of organic soil, rather than an argument against it. Smoothing of soil surfaces relative to variability in underlying mineral layers or bedrock is observed in other wetland systems where soil creation is dominated by organic matter accumulation (e.g., the Everglades, Watts et al., 2014). This implies that deviations from this smooth organic soil surfaces are related to other surface-level processes, such as spatial variation in organic matter accumulation resulting from hypothesized elevation-productivity feedbacks. 475 Hummock heights relative to mean site-level water level were approximately 30 cm, aligning with field observations of relatively constant hummock height within sites. Generally consistent hummock height across sites in conjunction with clear bimodality in soil elevations supports the contention that hummocks and hollows are discrete, self-organized ecosystem states (sensu Watts et al., 2010). However, variability in site-level hummock heights-especially at depression and transition sites-may partially be attributable to hummocks in non-equilibrium states. From our 480 feedback model (Figure 1), it seems reasonable that within a site, some hummocks may be in growing states (e.g., increasing in height over time via the elevation-productivity positive feedback) and some may be in shrinking states if hydrologic conditions have recently become drier (e.g., decreasing in height via the elevation-respiration negative feedback), the combination of which may result in a distribution of hummock heights centered around an equilibrium hummock height. Future efforts could leverage time-series observations of hummock properties (e.g., area, height and 485 volume), but we note the likely decadal time-scales required to detect hummock growth or shrinkage (Benscoter et al., 2005;Stribling et al., 2007).
Local hydrology exhibited clear control on hummock height, providing evidence for our hypothesis that hummocks are a biogeomorphic response to hydrologic stress in wetlands. We found support for this contention at both the site level and at the hummock level. The tallest hummocks were consistently located in the wettest sites and in the wettest 490 zones within sites. At the site-scale, 85% of the variance in average hummock height could be explained by mean water level alone. Within sites, local mean water level explained on average 35% of the variability in hummock height ( Figure 7); prevalence of non-equilibrium hummock states may explain much of the additional variability. The considerable variation in the ability of local water levels to explain hummock height within sites (adjusted R 2 =0.12-0.56), and in the strength of that relationship (linear regression slopes=0.4-1.1) may be attributed to two factors: 1) 495 the across-site flat water level assumption, and 2) lack of long trends for hydrology. The flat water level assumption is likely to be a minor effect in transition sites with deep organic wetland soils (e.g., Nungesser, 2003;Wallis and Raulings, 2011;Cobb et al., 2017) but could be significant at depression and lowland sites with shallower O horizons.
Lack of sufficient data to characterize mean water level may also be an issue at several of our sites, because hummocks likely develop over the course of decades or more, whereas our hydrology data only span three years. To our 500 knowledge, this study represents the first empirical evidence of the positive relationship between hummock height and hydrology in forested wetlands. These results are consistent with previous research on tussocks of northern wet meadows (Peach and Zedler, 2006;Lawrence and Zedler, 2011) and shrub hummocks in brackish wetlands (Wallis and Raulings, 2011). The concordance in hydrologic control in these disparate systems suggests a common mechanism of (organic) soil building and accumulation on hummocks that may result from increased vegetation growth from 505 reduced water stress and/or from transport and accumulation of nutrients (Eppinga et al., 2009;Sullivan et al., 2011;Heffernan et al., 2013, Harris et al., 2019.

Controls on microtopographic patterning
We found clear support for our hypothesis that hummocks are non-randomly distributed in our wettest study sites.
Hummocks exhibited spatial overdispersion at all sites, but this overdispersion was only significant at depression and 510 transition sites (Figure 8). Significant spatial overdispersion indicates regular hummock spacing in contrast to clustered distributions or completely random placement. Regular patterning of landscape elements is observed across climates, regions, and ecosystems (Rietkerk and van de Koppel, 2008), and is indicative of negative feedbacks that limit patch expansion (Quinton and Cohen 2019). Our results indicate similar patterning for forested wetland microtopography and, importantly, demonstrate the hydrologic controls on that patterning. Hydrology appears to be a 515 common driver in regular pattern formation in wetlands (Heffernan et al., 2013) and drylands (Scanlon et al., 2007).
Thus, water stress-both too much (Eppinga et al., 2009) and too little (Deblauwe et al., 2008;Scanlon et al., 2007)appears to be an important regulator of patch distribution across the landscape.
We observed lognormal hummock size distributions, suggesting that some hummocks may attain very large areas (i.e., over 10 m 2 ), but the majority of hummocks (~80%) are less than 1 m 2 (Figure 9). This finding aligns with field 520 observations, where most hummocks were associated with a single black ash tree, but some hummocks appeared to have merged to create large patches. Truncated patch size distributions are common in other systems as well, such as the stretched exponential distribution for geographically isolated wetlands (Watts et al., 2014) or the lognormal distribution for desert soil crusts (Bowker et al., 2013). These types of distributions have fewer large patches than would be expected for systems without patch-scale negative feedbacks, and have a central tendency towards a common 525 patch size. Hence, truncation in hummock size distributions comports with hypothesized patch-scale negative feedbacks (i.e., tree competition for light and/or nutrients) that inhibit expansion. Hummocks at drier lowland sites did not conform to size distributions for wetter depression and transition sites, supporting our hypothesis that the feedbacks that control hummock maintenance and distribution are governed by hydrology and amplified in wetter conditions.. This work adds to recent efforts across climates and systems to use patch size distributions to infer drivers 530 of ecosystem self-organization and response to environmental conditions (Kefi et al., 2007;Maestre and Escudero, 2009;Weerman et al., 2011;Schoelynck et al., 2012;Tamarelli et al., 2017).
Characteristic hummock sizes in association with overdispersion in black ash wetlands suggest that hummocks are laterally limited in size by negative feedbacks on the scale of meters (Manor and Shnerb, 2008). We posit that there are two patch-scale negative feedbacks: 1) overstory competition for nutrients and 2) understory and overstory 535 competition for light. Hummocks associated with black ash trees, which account for more than 85% of measured hummocks, are likely limited in area by the radial growth of the tree's root system. Evapoconcentration feedbacks bring nutrients to the tree roots, limiting the degree to which roots must search for them (Karban, 2008), and therefore limiting root lateral expansion. Indeed, evidence suggests that a majority of fine tree roots occur within hummocks in forested wetland systems (Jones et al., 1996;Jones et al., 2000). Moreover, finite nutrient pools may lead to 540 development of similarly sized nutrient source basins for each hummock, further limiting lateral hummock expansion Eppinga et al., 2008). Black ash trees must also compete for light with other ash trees, but leaf area is typically low in these systems (Telander et al., 2015). Low LAI and observed crown shyness (sensu Long and Smith, 1992) in black ash wetlands may imply less competition among individuals than would be expected in mixed stands (Franco, 1986). On the other hand, less-than-expected canopy competition for light in the overstory may 545 increase light availability for understory hummock species, and therefore allow subsequent hummock expansion from the understory. Therefore, based on evidence and observations presented here and in Diamond et al. (2019), we suggest that a major difference between microtopography in forested versus non-forested wetland systems will be the size distributions and spacing of hummocks. In other forested systems, hummocks associated with trees will likely be limited in size, exhibiting characteristic sizes and spacing due to local negative feedbacks from the crown competition. 550 In contrast, non-forested wetland hummocks may have a much wider distribution of size classes, where negative feedbacks to hummock expansion may be largely due to local nutrient competition effects (e.g., Eppinga et al. 2008).

Evidence for patch self-organization
In this work, we used common landscape ecology diagnostics to characterize microtopographic pattern and infer responsible reinforcing processes, including analyses of multimodal distributions of elevation, spatial patterns of 555 hummock patches, and hummock size distributions. Other recent work has used nearly identical diagnostic measurements to infer self-organization of depressional wetland features (~100 m wide) in a karst landscape (Quinton and Cohen 2019), demonstrating the broad utility of the approach and the various spatial scales that pattern may manifest. However, we note that this diagnostic approach alone does not directly implicate hypothesized mechanisms of hummock persistence, and that more measurements are required to support inferences made here. To that end, in 560 complementary work we observed support for the elevation-productivity feedback, where we found hummocks to be loci of higher tree occurrence and biomass, more understory diversity, and greater phosphorus and base cation soil concentrations (Diamond et al., 2019). Further, these associations were most evident in the wettest sites, concordant with the hydrologic controls observed here for hummock height, pattern, and size distributions. Together, these multiple lines of evidence lend strong support for the hydrologically driven self-organization hypothesis of hummock 565 growth and persistence (Figure 1).

Broader implications
The consequences of wetland microtopography are clear at small scales, but can also scale to influence site-and regional-scale processes. For example, microtopographic expression results in a drastic increase in surface area within wetlands. We conservatively estimate an average of 22% and up to 42% relative increase in surface area due to the 570 presence of hummocks (i.e., additional surface area provided by the sides of hummocks; Table 3). These estimates comport with studies in tussock meadows, where tussocks (ca. 20 cm tall) increased surface area by up to 40% (Peach and Zedler, 2006). Further, increases in the diversity of biogeochemical processes occurring at the individual hummock or hollow scale (Deng et al., 2014) likely aggregate to influence ecosystem functioning at large scales. For example, microtopographic niche expansion allows for local material and solute exchange between hummocks and 575 hollows, creating coupled aerobic-anaerobic conditions with emergent outcomes for denitrification (Frei et al., 2012) and carbon emission (Bubier et al., 1995;Minick et al. 2019ab).
While our results implicate hydrology as a major determinant of microtopographic structure and pattern, microtopography can reciprocally influence system-scale hydraulic properties. Results from our hummock property analysis indicate that hummock volume displacement may be a significant factor in water level dynamics of wetlands. 580 Specific yield, which governs water level response to hydrologic fluxes, is commonly assumed to be unity when wetlands are inundated. However, inclusion of microtopography may render this assumption invalid, with hummock volumes up to 30% of site volumes (Table 4). These observations are supported in other studies of microtopographic effects of specific yield (Sumner, 2007;McLaughlin and Cohen, 2014;Dettmann and Bechtold, 2016). Therefore, while hydrology exerts clear control on the geometry of hummocks, hummocks may exert reciprocal control on 585 hydrology by amplifying small hydrologic fluxes into large water level variations.
Last, black ash hummocks provide unique microsite conditions that support increased vegetation growth and diversity (Diamond et al. 2019), aligning with observations in other wetland systems (Bledsoe and Shear, 2000;Peach and Zedler, 2006;Økland et al., 2008). Accordingly, recent wetland restoration efforts have begun to use microtopography as a strategy to promote seedling success and long-term project viability (Larkin et al., 2006;Bannister et al., 2013;590 Lieffers et al., 2017). Specific to our focal system, there are increasing efforts to mitigate potential black ash loss due to the emerald ash borer and possible regime shifts to marsh-like states (Diamond et al., 2018). We posit that hummock presence and persistence may allow for future tree seedlings to survive wetting up periods following this ash loss (Slesak et al., 2014), and for consequent resilience of forested ecosystem states.
Overall, this study adds to the growing body of evidence that the structure and regular patterning of wetland 595 microtopography is an autogenic response to hydrology. Although the imprint of biota on landscapes may be masked by the signature of larger scale physical processes (Dietrich and Perron, 2006), we show clear evidence here for a microtopographic signature of life.

5
Code and data availability The authors will provide code and data upon request, and will upload code to Github upon acceptance of the 600   Figure 1. Conceptual model for autogenic hummock maintenance in wetlands. Incipient mechanisms create small-scale variation in soil elevation that is amplified by autogenic feedbacks, which grow and maintain elevated hummock structures. Solid lines indicate positive feedback loops and dashed lines indicate negative feedback loops. Font in italics refer to feedback processes hypothesized to only affect lateral hummock extent (thus hummock area), whereas non-italic font indicates mechanisms that affect both vertical and lateral hummock extent. Processes in blue indicate that these 850 mechanisms are influenced by hydrology. Soil mass refers to the amount of (organic) soil in a hummock, which can include roots, leaves, and decaying organic matter.      best-fit normal distributions (red lines). Text indicates the mean nearest-neighbor distance (μNN ± standard error); the ratio of the measured mean nearest-neighbor distance and the expected nearest neighbor distance for complete spatial randomness (μexp); and the p-value for a z-score comparison between μNN and μexp. p-values less than 0.001 indicate that hummocks are significantly overdispersed. Figure 9. Inverse cumulative distributions of hummock dimensions (perimeter, area, and volume) across sites (points), split by hummock dimension and site type. The y-axis is the probability that a hummock dimension value is greater than or equal to the corresponding value on the x-axis. Best-fit lognormal distributions are shown for each site as lines. All fits were highly significant (p<<0.001). Text indicates mean (±sd) within-group coefficient for a model of the form P(X≥x)=β*ln(dimension_value) + β0.