Quantifying Topographic Characteristics of Wetlandscapes

Topography underpins natural processes ranging from incident solar radiation to overland flow and water pooling. Despite the influence of topography on natural processes, especially in wetland ecosystems reliant on uplands for water inputs, topography has not been adequately incorporated into reclamation planning. Instead, wetland reclamation projects are typically guided by height-to-length ratios that produce little resemblance to natural wetlands. We present a methodology to quantify the topographic characteristics in landscapes with an abundance of wetlands to guide the reclamation of naturally appearing and self-sustaining wetlandscapes. Topographic characteristics in over 3000 sample landscapes were quantified using terrain roughness and landform element composition and configuration. A large set of metrics was reduced to a statistically independent subset that was applied and compared across three natural regions and a gradient of disturbance. Our results showed that surface roughness and landform element patterns significantly differ among natural regions and that high disturbance landscapes significantly differ from other disturbance levels. To ensure reclaimed wetlandscapes are naturally appearing and self-sustaining, they should replicate the topographic characteristics found within the distribution of surrounding natural landscapes by applying topographic characteristic benchmarks to reclamation design. The presented methodology can be used as a step towards achieving this goal.


Introduction
Characteristics of the surface of Earth (e.g., slope, aspect) and variation in elevation affect erosion (Martín-Duque et al. 2010), hydrology (Los Huertos and Smith 2013), and other natural processes that subsequently influence the distribution of biotic communities across local-to-regional landscapes (Nellemann and Fry 1995;Hofer et al. 2008). Through terrain analysis (Wilson and Gallant 2000;Hengl and Reuter 2009), topographic characteristics and variation in elevation are quantified to better understand how topography and landform elements underpin natural processes and the spatial distribution of ecosystems.
Wetland ecosystems are an archetype of natural systems often reliant on surrounding topography as a driver of wetland hydrology (Los Huertos and Smith 2013). These ecosystems provide a range of services that include: important waterfowl habitat (Conly and van der Kamp 2001) and environmental regulation including storm water mitigation, water purification, and carbon sequestration (Zedler and Kercher 2005). Despite the benefits of wetland ecosystem services, global wetland loss is estimated at 64-71% since 1900 AD (Davidson 2014) and these ecosystems are difficult to restore back to natural biological structure and function (Moreno-Mateos et al. 2012).
Ecological restoration is practiced to restore important ecosystem services to a region, of which wetlands have been a prevalent focus (Wortley et al. 2013). However, it is difficult to assess the success of restoration projects (Wortley et al. 2013) and not all projects result in healthy and resilient wetlands (e.g., Kauffman-Axelrod and Steinberg 2010) or even functional wetlands (Moreno-Mateos et al. 2012). Unsuccessful restoration projects often occur by focusing on restoring wetland biota without considering the underlying hydrology and terrain conditions in a region (Kauffman-Axelrod and Steinberg 2010). Similarly, wetland restoration focused on individual wetlands, referred to as at the wetland scale (e.g., Newcomer et al. 2013), is likely to fail when it is not integrated into the broader landscape (Van Meter and Basu 2015).
The success of restoration projects is dependent on the degree to which the restored wetland function is integrated with the surrounding landscape (White D, Fennessy S 2005). A key driver of wetland function is the underlying topography of the landscape as it affects hydrological processes and water availability (Los Huertos and Smith 2013). For example, the topographic position (e.g., toeslope or depression) of a wetland within a landscape can influence the source (overland flow, ground water flow, or precipitation) and amount of water received.
Where topography and ecosystems have been degraded by anthropogenic disturbance, solely restoring ecosystem biota without reshaping the landscape to create a natural topography is not sufficient to achieve a stable and functionally similar ecosystem that existed prior to disturbance can be an unacievable goal (Suding 2011;Wortley et al. 2013;Rooney et al. 2015). Instead, the movement of soil to reshape the landscape and improve degraded land to facilitate future land use and ecosystem function, which we refer to as reclamation (Burton 1991), can improve the chances of creating a stable and functional ecosystem. However, reclamation guidelines for wetlands typically advise for the creation of gentle slopes surrounding a wetland at a 2:1 or 3:1 height-to-length ratio and irregular shorelines for the wetland itself (e.g., Online Resource 1). These guidelines lead to the homogenization of wetlandscapes with evident anthropomorphic history that lack natural topographic variation and pattern. Instead, to ensure reclaimed wetlands or more broadly wetlandscapes are naturally appearing and self-sustaining, guidelines should facilitate the replication of topographic characteristics found within the distribution of surrounding natural wetlands or wetlandscapes. By implementing mature landforms, reclaimed landscapes can establish natural channel and surface morphology that have a greater resistance to erosion from precipitation and discharge events while avoiding expensive, continuous maintenance and regrading (McKenna 2002;Martín-Duque et al. 2010).
While the reference condition approach has been used to create benchmarks for indicators of biological integrity (Bailey et al. 2004), river restoration (Nestler et al. 2010), and the composition and configuration of land cover (Evans et al. 2017) to evaluate the health of ecosystems and guide ecosystem reclamation and regulations (Stoddard et al. 2006), the approach has its limitations. The reference condition approach is sensitive to variations in sampling conditions between the reference state and future comparisons. For example, seasonal or spatial variation between the reference condition and current sampling may introduce error into an analysis, which could be compounded by instrument, measurement, or human error (Reece and Richardson 1999;Bailey et al. 2004).
In addition, the variables used to determine the pristine nature of an ecosystem are a function of the dominant social and cultural values (Rapport 1989) and therefore are likely to include some level of bias since not all attributes of an ecosystem can be measured.
Unlike other applications of the reference condition approach, no benchmarks have been established for topographic landforms or variation in wetlandscapes. However, recent research has demonstrated the applicability of topography to wetlandscapes, where it was found that the statistical distribution of real wetlands' size and shape were similar to the size and shape distributions of potential wetlands identified using topographic depression analysis (Bertassello et al. 2018). The adoption of topographic benchmarks by reclamation planners would increase the potential for designed landscapes, with diverse topographic features, to integrate with the surrounding natural landscape, be naturally appearing and self-sustaining, and lead to successful reclamation.
To improve our understanding of the relationship between wetlands and their surrounding landscapes and to assist resource managers and policy makers in reclaiming wetland ecosystems, we quantified the spatial pattern and characteristics of topography within wetlandscapes to answer the question: to what degree do topographic characteristics and variation in elevation differ by natural region and along a gradient of human disturbance?

Study Area
Wetlands have been decreasing in number and areal extent globally (Davidson 2014) and Canada is no different (Ducks Unlimited Canada 2008). Within Canada an estimated 70% of wetlands have been lost (Ducks Unlimited Canada 2008) and the Prairie Pothole Region of the Province of Alberta (Shaw et al. 2013) has experienced wetland consolidation from agriculture and, while estimates vary, an approximate wetland area loss of between 68% (Serran et al. 2017) and 90% of historic wetland area (Van Meter and Basu 2015). To ameliorate wetland loss, the Government of Alberta has developed the 'Water for Life' action plan, with one among several goals of ensuring healthy aquatic ecosystems, which requires the conservation and restoration of Alberta's wetlands through an improved understanding of wetland functions (Alberta Environment 2003;Alberta Water Council 2008).
We situate our research within the Grassland, Parkland, and Dry Mixedwood subregion of the Boreal natural regions of Alberta, Canada, which are within the prairie and boreal plains Canadian ecozones. Within these natural regions, our research is specifically situated within two wetland inventories that cover 27.8% of the province (Fig. 1). The inventories are located within the more populated southern and central regions of the province, the majority of which lies within the Prairie Pothole Region of North America. This region is dominated by mineral wetlands that are reliant on seasonal snow melt as a primary water input (Shaw et al. 2013).

Methodology Overview
To quantify topographic characteristics and variation in elevation in wetlandscapes, our analysis is divided into 6 conceptual steps (Fig. 2). We begin with a selection of sample landscapes and subsequently group sample landscapes by natural region and proportion of disturbance. These steps are followed by quantifying topographic characteristics and landform elements, calculating spatial configuration metrics, reducing the calculated metrics to a parsimonious and orthogonal set, and finally testing for differences in metric values across sample landscape groups (Fig. 2).
Key methodology and analysis steps are outlined to provide a clear overview of the process to classify and characterize the topography of wetlandscapes. The sample landscapes with our study area were grouped based on the natural region they are located within (grassland, parkland, or boreal) as well the landscape's percentage of human disturbance (0-20%, 20-40%, 40-60%, 60-80%, and 80-100%) to create 15 landscape groups. Topographic roughness and the spatial configuration of landform Fig. 1 Study location in the Boreal, Parkland, and Grassland regions of southern Alberta. The study encompassed the majority of southern Alberta, Canada east of the Rocky Mountains and extended north into the boreal natural region north of Edmonton. A gap in the study area existed east of Edmonton due to GIS data availability. The study area was randomly sampled using 3434 1 km 2 study sites that provided the basis for the comparative topographic analysis elements were calculated for each sample landscape using a variety of roughness metrics and landscape metrics, respectively. The large set of calculated metrics were reduced to an orthogonal subset using pairwise correlation grouping and principal component analysis. The resultant topographic roughness and landform element pattern metrics for each group of sample landscapes were statistically compared between landscape groups to determine if landscape-scale topography statistically varied by human disturbance or natural region

Data and Sample Design
Terrain analysis theoretically only requires elevation data. Elevation data were acquired from AltaLIS Ltd. (AltaLIS 2015) with a 10 m spatial resolution, which was based on grid points, break lines, and spot heights that were obtained from 1:60,000 aerial photography flown between 1980 and 1995. AltaLIS performed quality control and periodic updates to the DEM and flattened the elevation of all lakes (AltaLIS 2015). After obtaining the DEM, elevation values were smoothed using two passes of a 5 × 5 cell mean filter to reduce noise in the data (e.g. MacMillan et al. 2003;Buchanan et al. 2014).
While the DEM does not contain water body bathymetry, its 10 m resolution is fine enough to enable the representation of topographic changes that may influence wetland ecosystems in the region (MacMillan et al. 2003).
Wetlands do not exist as an ecosystem in isolation from the surrounding landscape. The land cover surrounding wetlands, particularly in hydrologically 'upstream' regions, are influential to the water availability and subsequent health of wetland ecosystems (Conly and van der Kamp 2001). Two types of land cover that can negatively impact a wetland's upstream region, and are dominant in southern Alberta, are agricultural and urban lands, referred to as disturbance in this research. Land-cover data for Alberta from Agriculture and Agri-Food Canada were acquired and the proportion of disturbance (areas classified as 'agriculture', 'pasture', or 'urban') was calculated for each sample landscape (Evans et al. 2017).
Within the study area, 4000 1 km 2 square sample landscapes were previously generated at random and only those that contained wetlands were retained, resulting in 3434 sample landscapes (Evans et al. 2017). The 1 km 2 square sample landscapes are sized to match the average size of the surrounding landscape that impacts a wetland (500 m radius around a

Landscape Grouping
Landscapes with similar characteristics have been shown to behave functionally similar and identifying landscapes that have similar form and function can assist in landscape planning (de Groot et al. 2010). For example, landscapes with similar characteristics have measurable impacts on water quality (Huang et al. 2013), can be used to predict vertebrate species abundance (Mazerolle and Villard 1999), and pest abundance (Zaller et al. 2008). To generalize these landscapes, develop a parsimonious set of criteria that can be used by regulators and industry, and simplify their application, landscapes can be grouped together based on a variety of landscape-level characteristics.
Previous research demonstrated that the composition and configuration of land cover in our sample landscapes differed by natural region and disturbance (Evans et al. 2017). To coincide with Evans et al. (2017), 15 discrete landscape groups were created by partitioning the landscapes by five disturbance levels in each of three natural regions of Alberta: Grassland, Parkland, and Boreal Dry Mixedwood subregion (Table 1).

Terrain Roughness
Terrain analysis methods are used for relational and descriptive purposes. Relational analyses seek to improve our understanding of the relationship between topography and environmental variables. For example, the relationship between elevation and environmental variables such as soil depth and type (Florinsky et al. 2002), soil wetness (Murphy et al. 2011), snow depth (Lapen and Martz 1996), and potential for ecosystem restoration (Kauffman-Axelrod and Steinberg 2010) have been quantified. While relational analysis is not the focus of this paper, correlation and regression analyses have been used to quantify the relationship between these two ratio datasets to test if a relationship exists between topography and an environmental variable (e.g., Murphy et al. 2011;Kraft et al. 2019).
When the terrain analysis is conducted for descriptive purposes, which is the purpose of the presented research herein, a different set of methods is required to quantify the variation in topography across an area. Despite a wide range of metrics used in terrain analysis (e.g. Olaya 2009;Florinsky 2012) to quantify elements of topography, from individual cells (elevation value), to local cell windows (slope), and regional relationships (flow accumulation), none are designed specifically for application at a landscape scale. Quantifying the topographic characteristics of a landscape is a challenge (MacMillan et al. 2000;Deng et al. 2007), which is partly due to the continuous surface of ratio data representing elevation. To overcome this challenge, variation in elevation has been 1) described aspatially using measurements of terrain roughness (Hengl and MacMillan 2009) or 2) transformed into nominal data describing landform elements (Reuter et al. 2006). The link between terrain roughness and ecological process emphasizes the need to design heterogeneous landscapes that mimic natural surface roughness to improve the success rate of reclamation projects.
To find suitable terrain analysis methods to quantify terrain roughness, a literature review of terrain roughness metrics (TRMs) was conducted by the authors that identified 8 metrics designed to quantify the variation in elevation that could be conceptually applied at a landscape scale (Table 2) using spatially continuous data (Grohmann et al. 2011). In total, we used seven TRMs to quantify surface roughness in each sample landscape (see Online Resource 2 for a detailed description of metric calculations and modifications). We excluded the topographic position index from our analysis as it provides a less accurate representation of topographic variation than the deviation from mean elevation (DEV) metric (De Reu et al. 2013). For five of the implemented seven metrics, DEV, standard deviation (SD) elevation, SD profile curvature, SD slope, and slope variability, two variations on the metric calculations were performed: 1) the metric was calculated based on all cells in the entire sample landscape (i.e., global), and 2) the metric was calculated using a moving 5 × 5 focal window and then the average of the window values for the entire sample landscape was recorded. Two calculation scales were used due to the conceptual difference in calculation methodology, which resulted in a noticeable difference in metric values between the two scales. The focal window method calculates the average of small-scale variations across a landscape, which may not capture larger-scale topographic changes. Conversely, the global method calculates the variation in metric values across a landscape without requiring values to be averaged, resulting in larger topographic roughness values. The focal 2D:3D Area Ratio and Inverse Vector Dispersion metrics could not be directly calculated for the entire landscape so they were only calculated using the focal window method, resulting in 12 total metric calculations. Summary statistics of the average and standard deviation for each TRM were calculated by natural region and disturbance level (Online Resource 3). Each landscape group's TRM distribution was tested for normality using the Shapiro-Wilks test and visually assessed using QQ plots, where it was determined that all TRMs had a non-parametric distribution for all landscape groups.
Next, the 12 TRMs were reduced by grouping metrics that shared a Pearson's correlation coefficient > 0.9 with all other members of the group (Riitters et al. 1995;Moreno-Mateos et al. 2008). A representative metric was selected from each correlation group and the process repeated until all representative metrics fell below the correlation threshold (Riitters et al. 1995). Through the iterative comparison of each metric pair's correlation value, a final set of five TRMs was selected.
A Kruskal-Wallis statistic was used to test for significant differences between landscapes grouped by natural region and disturbance. In the case where a significant difference was found, a post-hoc Dunn's test was used to perform pairwise comparison tests between landscape groups. A Bonferroni correction was also applied to align with previous research (Evans et al. 2017); control for Type 1 errors, which occur when a test incorrectly reports a significant difference when none exist (Dunn 1961); to increase analytical rigor since we did not measure spatial autocorrelation; and to create statistically-sound guidelines for industry.

Landform Classification Methodology
Describing topographic change across a landscape is difficult to do quantitatively. The use of surface roughness and other summary statistics are aspatial and lose both information about topography at a specific location (i.e., site characteristics) and the relationship of topographic change between the location and neighbouring locations (i.e., situation). To quantify the pattern of topography in our sample landscapes, we classified the topography into different landforms and then used landscape metrics to describe the pattern of those landforms.
A landform is 'a physical feature of the Earth's surface having a characteristic, recognizable shape and produced by natural causes' (pg 228, MacMillan and Shary 2009). Landform classification has been widely applied to different regions (e.g., Turkey, Tagil and Jenness 2008; Germany, D'Oleire-Oltmanns et al. 2013) and scales (e.g., microtopographic hummocks and hollows (Nagamatsu and Miura 1997) to continental-scale mountains and plains (Schmidt and Hewitt 2004)). The purpose of the classification is, typically, to describe the composition of landforms in a study area in relation to an environmental variable of interest (e.g. Florinsky 2012). Apart from Sărăşan and Ardelean (2015) who used landscape metrics to differentiate landform patterns, the authors are unaware of research quantifying the spatial configuration of landforms. Although the spatial configuration of landforms has not been quantified, characterizations of wetlandscapes in our study area of southern Alberta has been completed to better understand the hydrological functions of prairie wetland (Hayashi et al. 2016) and the intermittent hydraulic connectivity between individual wetland depressions (Shaw et al. 2013).
The use of 1 km 2 sample landscapes situate our classification at a scale for defining landform elements (Schmidt and Hewitt 2004) rather than traditional landforms such as drumlins, plains, or mountains. Landform elements are small regions with homogeneous terrain metric values that occur along a hillslope (MacMillan et al. 2000;Schmidt and Hewitt 2004). A hillslope can be composed of multiple landform elements, such as a peak, shoulder, toeslope, and depression (Reuter et al. 2006).
The presented research uses a discrete (or rule-based) approach to delineate landforms based on classification rules for   (Pennock et al. 1987;Reuter et al. 2006). While the classification rules need to be adapted for the region in which the landform classification is occurring (Reuter et al. 2006), their use enables 1) transparency and ease of use, 2) comparison across different regions, and 3) builds a compendium of examples for meta-analysis. Our discrete landform classification used geometric terrain metrics of slope, profile curvature, plan curvature (Pennock et al. 1987;MacMillan et al. 2000) and deviation from mean elevation (DEV) (Lindsay et al. 2015). These four terrain metrics were combined to create 12 distinct landform elements, 11 that correspond to those by Reuter et al. (2006) and one taken from (MacMillan et al. 2000; Table 3, Online Resource 4).

Quantifying Topographic Composition and Configuration
To quantify the composition and configuration of landform elements and characterize topographic patterns by natural region and disturbance level, landscape metrics were applied to our sample landscapes. While a variety of approaches to conceptualizing space are possible, landscape metrics offer one approach to statistically quantifying the pattern of categorical data. Although 100 s of landscape metrics have been developed, they have been categorized into those that describe areaedge, shape, aggregation, and the diversity of landscape categories (McGarigal et al. 2012), which may be applied to an individual patch, a class, or at a landscape scale (McGarigal et al. 2012). Many landscape metrics quantify a relationship between ecological processes and land cover patterns (e.g., core area and contrast, McGarigal et al. 2012). Since these relationships are irrelevant to landforms, all core area and contrast metrics were excluded from our analysis. Thirty-two landscape metrics were deemed applicable to landform elements and were calculated at the landscapescale (as opposed to patch-or class-scale) for each sample landscape (Online Resource 5) using FRAGSTATS (McGarigal et al. 2012). Sample landscapes were grouped based on natural region and disturbance level (Table 1). Each landscape metric was tested, by group, for normality using the Shapiro-Wilks test and visually assessed using QQ plots, which determined that all metrics had a nonparametric distribution. The thirty-two metrics were reduced to a subset of non-correlated landscape metrics using the same correlation analysis described above for terrain roughness. Further metric reduction was conducted using an approach borrowed from bioclimate envelope modelling (Metzger et al. 2013), whereby the remaining metrics were included in a principal component analysis (PCA) that generated a new set of orthogonal variables that are known as components. The components were not explicitly used, instead, the factor loadings of the remaining metrics were evaluated and only the metrics with the highest absolute factor loading for each principal component that had an eigenvalue greater than one were retained (similar to Mladenoff et al. 1997 andEvans et al. 2017). Using the PCA factor loadings, a final subset of landscape metrics was identified and used in subsequent analyses. A Kruskal-Wallis test was used to identify significant differences between metrics by natural region and disturbance level. Where significant differences were found, a post-hoc Dunn's test with Bonferroni correction was used to perform pairwise landscape metric tests between landscape groups.

Terrain Roughness Results
Five representative terrain roughness metrics (TRMs) were selected based on our iterative correlation analysis (Online Resource 6). The correlation analysis was repeated for each natural region separately to determine if the Resource 3) and analysis of variance tests were applied across natural regions and disturbance levels (Table 4). We tested for significant differences among natural regions within the same disturbance level using a Kruskal-Wallis test for each of the five representative TRMs (Table 4, comparison #1). Of the 25 tests, 12 (48%) were significantly different at p < 0.001, 4 (16%) were significantly different at p < 0.01, 3   .05, which with a bonferroni correction applied were deemed as not statistically significant in our analysis (12%) were significantly different at p < 0.05, and 6 (24%) had p values >0.05, which for the presented analysis with bonferroni correction were deemed as not significantly different. Focal slope variability, global slope variability, and area ratio metrics were significantly different among natural regions at p < 0.05 for all five disturbance levels, whereas global DEVand inverse vector dispersion were significantly different for only two of the five disturbance levels (Table 4, comparison #1). The remaining two metrics, global DEV and inverse vector dispersion, were not congruent but still demonstrated sensitivity to natural region. Our results, that 76% of our tests have significantly different metric values, suggest that surface roughness differs among natural regions regardless of the level of disturbance in the landscape.
To identify specific differences in TRMs between natural regions across all levels of disturbance, a pairwise Dunn's test was conducted. Of the five TRMs, focal and global slope variability as well as area ratio were significantly different between all three natural regions (Table 4, comparison #2). The remaining two TRMs, global DEV and vector dispersion, were only significantly different between Parkland and Grassland landscapes. These results confirm that significant differences in surface roughness exist between natural regions; however, these differences are less pronounced between the Grassland and Parkland regions.
To evaluate if surface roughness differed with disturbance level, a Kruskal-Wallis test was conducted separately for each natural region for the five representative TRMs based on landscapes grouped by proportion of disturbance. The Kruskal-Wallis test identified that a significant difference existed between landscapes at different disturbance levels for all five TRMs at p < 0.001, except for the inverse vector dispersion metric in the Parkland (significant at p < 0.01) and Grassland (not significant) regions (Online Resource 7). These results suggest that landscapes with different levels of disturbance within the same natural region have significantly different surface roughness.
To identify specific differences in surface roughness between disturbance levels, a Dunn's posthoc pairwise comparison test was performed for each of the three natural regions. For the Boreal region, the test identified that high disturbance landscapes (80-100% disturbance) were significantly different than all other landscapes for all five representative TRMs (Table 4, comparison #3). There was also a significant difference between landscapes with 0-20% and 20-40% disturbance for 4 of the 5 metrics and focal slope variability was significantly different between 20 and 40% and 60-80% disturbance landscapes (Table 4, comparison #3).
Results for Parkland landscapes found significant differences between high disturbance landscapes (80-100% disturbance) and all other landscapes for at least 3 of the 5 TRMs (Table 4, comparison #4). Landscapes with 0-20% in disturbance were also found to be significantly different than 40-60% landscapes for 4 of the 5 TRMs (Table 4, comparison #4).
For the Grassland natural region, the pairwise comparisons found significant differences between high disturbance landscapes and all other landscapes at p < 0.001 for all TRMs except for focal inverse vector dispersion (Table 4, comparison #5). There were no significant differences found for any of the other pairwise comparisons (Table 4, comparison #5).
Differences in surface roughness were distinguishable between different levels of disturbance for the Boreal natural region, suggesting that the topography of the Boreal differs greatest from that modified by disturbance. The ability to distinguish differences in surface roughness among disturbance levels decreased in the Parkland and was least differentiated in the Grassland (Table 4, comparison #6). However, our pairwise comparisons across disturbance levels within each of the natural regions demonstrate congruence that surface roughness of 80-100% disturbance landscapes are significantly different than landscapes with other levels of disturbance regardless of natural region. Among the five TRMs evaluated, focal slope variability was the strongest differentiator across disturbance levels, showing significant differences in 15 of the 30 comparisons, which corroborates its influence in differentiating between natural regions (Table 4, comparison #6).

Landform Results
Each sample landscape's topography was classified into 12 discrete landform elements and were grouped based on their expected overland flow characteristics (Fig. 3). In general, low and high level landform elements, which are flat ground located below or above the average landscape elevation, are the most prevalent landform elements across all landscapes, while convergent backslope and convergent footslope elements are the least prevalent. We quantify and compare the composition and configuration of these landform elements across natural regions and disturbance level.
Landform elements in each sample landscape were classified based on a defined rule set that incorporated four terrain metrics: topographic slope, profile curvature, plan curvature, and deviation from mean elevation. Each terrain metric was reclassified into 3 distinct levels (i.e. no slope (0°), low slope (0°-3°), and high slope (<3°)) and the four metrics were combined to create 12 distinct landform elements in a sample landscape (i.e. a divergent shoulder landform element consists of low or high slopes, convex profile curvature, convex plan curvature, and any deviation from mean elevation value)

Metric Reduction
Thirty-two landscape metrics were reduced using an iterative correlation analysis, which grouped all metrics with a paired Pearson correlation coefficient greater than 0.9 and then a representative from each group was selected. The correlation analysis resulted in 15 representative landscape metrics for the Boreal and Grassland regions, 14 metrics for the Parkland region, and 16 metrics when considering the entire region (Online Resource 8).
To further refine the set of landscape metrics, the representative metrics identified in the correlation analysis were included in a Principal Component Analysis (PCA). The metrics that retained the highest factor loading for each principal component, with an eigenvalue greater than one, yielded a set of five landscape metrics for the entire study area, four landscape metrics for each natural region, and a total of seven unique metrics (Table 5). Each natural region retained the same two metrics (SHEI and ENN_AM), while Boreal and Parkland both retained GYRATE_AM and Parkland and Grassland both retained PROX_AM (Table 5).
Out of the seven unique representative landscape metrics identified, three landscape metrics measured landform element aggregation, two metrics measured landscape diversity or evenness, one measured landform element area or edge, and one measured landform element shape (Online Resource 5). The seven representative metric's mean and standard deviation values were calculated to compare the distribution of metric values across landscapes based on both natural region and proportion of disturbance (Online Resource 3).

Relating Metrics to Natural Region and Human Disturbance
The representative landscape metrics for each sample landscape group were tested for normality using a Shapiro-Wilks test and visually compared using QQ-plots; each landscape category reported a p value of <0.01, suggesting the data were non-normally distributed. A Kruskal-Wallis test was conducted to determine if significant difference existed among the three natural regions for landscapes with the same proportion of human disturbance (Table 6, comparison a). Of the 35 total comparisons made, 17 (48.6%) were significant at p < 0.001, 5 (14%) were significant at p < 0.01, 5 (14%) were significant at p < 0.05, and 8 (22.9%) had a p value >0.05, which for the presented analysis we determine as having no significant Fig. 3 Landform elements within a sample landscape selected from the Boreal natural region with 60-80% of its total area disturbed difference (Table 6, comparison a). Two of the seven metrics were significantly different for all five disturbance groups (GYRATE_AM and SIDI) ( Table 6, comparison a).
Landscapes in different natural regions were statistically compared using Dunn's pairwise comparison test for the seven representative landscape metrics. The test found significant differences in landscape metrics between all three natural regions, and the greatest number of differences was found between Parkland and Boreal where six out of the seven representative metrics were significantly different (Table 6, comparison b). Parkland and Grassland had the least number of differences, with three out of the seven representative metrics significantly different (Table 6, comparison b). These results suggest that the configuration of landscape elements differ significantly between natural regions. Given significant differences between natural regions, we tested for significant difference between landscape disturbance levels within each natural region. Results showed significant differences between landscapes with different proportions of disturbance for all seven representative metrics apart from the CONNECT metric in the Boreal and Grassland regions and the SHEI metric in the Grassland region (Table 7).
Pairwise comparisons between disturbance classes in the Boreal using an adjusted post-hoc Dunn's test found significant difference between high disturbance landscapes and all other landscape groups for the GYRATE_AM and SIDI metrics. Furthermore, six of the seven representative metrics are significantly different between landscapes with 0-20% and 80-100% disturbance (Table 6, comparison c). Significant differences were also found for three of the seven metrics between 0 and 20% and 20-40% landscapes (Table 6, comparison c).
The Parkland natural region results also found that the high disturbance landscapes are significantly different from landscapes at all other disturbance levels ( Table 6, comparison d  74). At least four of the seven representative landscape metrics were significantly different between high disturbance landscapes and all other landscapes and six of the seven metrics were significantly different between landscapes with 60-80% and 80-100% disturbance. The only other difference detected was for the SIDI metric at p < 0.05 between 0 and 20% and 40-60% disturbance landscapes (Table 6, comparison d).
Our comparison between levels of disturbance within the Grassland region identified significant differences between high disturbance landscapes and all other landscapes, with no other significant differences found. Between three and five representative landscape metrics were significantly different between 80 and 100% disturbance landscapes and all other landscapes (Table 6, comparison e). Three metrics, GYRATE_AM, SHAPE_AM, and SIDI were significantly different across all four comparisons for high disturbance landscapes (Table 6, comparison e).
The results of our landform element analysis were similar to the results of the surface roughness analysis. Specifically, landscapes with 80-100% disturbance held a composition and configuration of landform elements that differed significantly from landscapes with all other levels of disturbance regardless of natural region. Landform elements were least differentiated across disturbance levels in the Grassland, suggesting that similar spatial patterns of landform elements exist in grassland regions with medium to low levels of disturbance. 14 of the 30 comparisons (Table 6) showed SIDI to be significantly different, and the most pronounced differentiator of landform element pattern across disturbance levels and natural regions.

Discussion
The presented results show that our representative roughness metrics were statistically different between natural regions in Alberta, Canada. Results were less systematic across disturbance levels, such that high disturbance landscapes (i.e., 80-100% disturbed) significantly differed from all other disturbance levels in each natural region and the degree of differentiation in surface roughness among disturbance levels was less pronounced in the Parkland and least pronounced in the Grassland. In all three natural regions, highly disturbed  (Yeh and Li 1998), which, when combined with conventional tillage practices in our study area, reduces the topographic variation. Similarly, development, and land-use change in general, decreases topographic variability (Lóczy and Gyenizse 2010). This flattening of topography by anthropogenic land-uses reaches a threshold that significantly differs for highly disturbed landscapes relative to all other disturbance levels. For example, differences in landscape element pattern are visually apparent when comparing landscapes at different disturbance levels (Fig. 4).
Landform elements in each sample landscape were classified based on a defined rule described for Fig. 3. All six landscapes are in the Boreal natural region and show examples of  landform element diversity across all five disturbance levels; a sixth landscape is added to show landform element diversity for a landscape with close to 0% disturbance. Shannon's diversity index values, from the landform element configuration metrics, are provided for each wetlandscape, where a value close to 1 exemplifies a landscape with a diverse range of landform elements. A roughness metric, focal slope variability, is also provided for each wetlandscape, where a higher value indicates that a landscape is topographically rougher and has greater elevation variability Our characterization of landscape topography showed different patterns in landform elements across natural regions and disturbance levels. From 32 metrics, we identified an orthogonal set of seven unique landscape metrics that can be used to quantify landform patterns. The seven representative landscape metrics were systematically different among natural regions, with only GYRATE_AM and SIDI statistically different across all five disturbance levels (Table 7). Within natural regions, high disturbance landscapes had significantly different landform element spatial pattern from landscapes at all other disturbance levels and no significant differences were observed between landscapes with less than 80% disturbance in the Grassland and Parkland, with one exception. Landscapes in the Boreal region were the exception, where 50% of landscape metric comparisons for landscapes with less than 80% disturbance were significantly different (Table 7). Landscape diversity, which measures landform element composition (SIDI), and an area/edge metric measuring landform element configuration (GYRATE_AM) were the only two metrics that differed significantly between high disturbance landscapes and all other landscapes across all natural regions. Which again was likely due to land leveling, which results in a low number of large landform elements and less landform element diversity. The landform element spatial pattern analysis corroborates the findings from the roughness analysis that high disturbance has less topographic variability, less complexity, and overall flatter landscapes.
Differences in topographic roughness, or diversity, has an impact on landscape-scale ecological function and individual ecosystem diversity. A landscape that is topographically rough (i.e., has greater elevation variation) has been demonstrated to impact species richness (Hofer et al. 2008), vegetation growth patterns, and wildlife behavior (Nellemann and Fry 1995). High species richness is an important characteristic of ecosystems that are resistant to disturbance (Hofer et al. 2008;Nellemann and Fry 1995), which suggests that designing heterogeneous landscapes that mimic natural surface roughness can improve the success rate and ecological stability of reclamation projects.
While landscape-scale topographic roughness has impacts on a landscape's ecological function, our quantification of landform elements and their diversity in wetlandscapes allows a specific comparison between topography and wetland function. Landform elements can be related to wetland function based on the hydrological position of wetlands within a landscape, which impacts wetland size, permanence, and salinity, which in turn can impact the diversity of vegetation growth (Hayashi et al. 2016). For example, if a wetland is located within an area classified as a high, level landform in our research, it likely relies primarily on precipitation for water inputs, with a small, low salinity, and less permanent pond (Hayashi et al. 2016). Conversely, a wetland located in an area we classify as a low, level landform receives more water inputs from upland areas during spring snowmelt runoff or during fill-spill events (Shaw et al. 2013;Hayashi et al. 2016;Wu and Lane 2016), which would increase a wetland's salinity and vegetation diversity (Hayashi et al. 2016). A range of wetlands located in different hydrological positions is necessary for an ecologically diverse wetlandscape (Hayashi et al. 2016). Our quantification of the spatial distribution of landform elements allows us to quantify the diversity of topographic features and a wetlandscape with more diverse topographic features and landform elements provides greater opportunities for wetlands in a variety of hydrological positions.

Use of Terrain Analysis for Reclamation
The presented research identified an orthogonal set of five surface roughness metrics and seven landform element pattern metrics that can statistically differentiate the terrain of landscapes in different natural regions and between high disturbance landscapes and those with less than 80% disturbance. Given that, one roughness metric, focal slope variability, and two landform element pattern metrics, SIDI and GYRATE_AM, were shown to be sensitive to topographic differences between natural regions and along a gradient of disturbance, they may be used as a simplified set to quantify and compare terrain roughness and spatial pattern of landscape elements.
The observed values for each metric (by natural region and disturbance level; Online Resource 3) can be used as a lookup table to guide reclamation planning and design by industry. Alternatively, the distribution of values could be used by regulators to evaluate reclamation closure plans to ensure reclamation landscapes fall within the range of topographic characteristics found in natural landscapes. However, oversight is required to ensure that the aggregation of reclaimed landscapes applying the presented topographic characteristic benchmarks represents the distribution of our metric values. If our benchmarks are used (e.g., just the median value) to generate a collection of homogeneous landscapes then our efforts have not advanced upon the height-to-length ratio guidelines (e.g., Online Resource 1) for wetland reclamation that fail to represent natural topographic variability.
Quantification of topography is one of many characteristics that affect wetland function and subsequently reclamation planning and acquisition of closure permits, e.g., soil and hydrological characteristics (McKenna 2002). Identifying the covariation among these different biophysical components remains a research gap and opportunity (de Groot et al. 2010) since the dependencies among topography and hydrology (Shaw et al. 2013), solar radiation (Hofer et al. 2008Mei et al. 2015), and vegetation diversity (Hofer et al. 2008) are strong. Incorporating covariation among the components of reclamation will improve the probability of reclamation success.

Spatial Resolution
Terrain analysis results have some dependency on the resolution of the input elevation data used (Thompson et al. 2001;Deng et al. 2007). For example, a coarser resolution can cause lower slope values in steep slope regions and higher slope values in flat slope regions while also decreasing the range in topographic curvature values (Thompson et al. 2001). While the presented research uses a 10 m resolution digital elevation model (DEM), the degree to which this resolution captures the true topographic variation within the prairie pothole region is unknown. Terrain analysis in geographically similar landscapes have used elevation data with resolutions ranging from 1 m (Van Meter and Basu 2015; Wu and Lane 2016) to 10 m (MacMillan et al. 2003). Recent advances in unmanned aerial vehicle (UAV) technology equipped with light detecting and ranging (LiDAR) or used in structure for motion derivation of surface elevation models may provide a cost-effective approach for quantifying the effects of spatial resolution on topographic representation. Comparison of the presented benchmarks across a range of spatial resolutions and DEM smoothing techniques (e.g. O'Neil et al. 2019) would add value to the presented research by quantifying the relative change in metric values by spatial resolution while also validating landform classification.
Similarly, the proportion of disturbance used in our study was calculated using AAFC annual crop inventory (land cover) data, which has a 56 m spatial resolution (Fisette et al. 2014). At this spatial resolution it is likely that we have both overestimated disturbance and the rectilinear representation of land cover patches in wetlandscapes. The degree to which our results are affected by this spatial resolution remains an avenue for future research.

Conclusion
Successful reclamation of wetlandscapes requires integrating reclamation projects with surrounding landscapes, of which topography is a key component. Topography is a driver of ecological processes such as water availability and sunlight potential (Los Huertos and Smith 2013), which are components in the health and location of wetland ecosystems. We quantified and described topographic patterns in wetlandscapes that yield a set of baseline topographic characteristics of wetlandscapes in Alberta.
A parsimonious set of terrain roughness and landform pattern metrics have been defined that can quantify the topographic variability of wetlandscapes and describe the spatial pattern of topographic landform elements across natural regions and at a gradient of disturbance. These metrics can be applied to landscapes to statistically compare topographic characteristics and to define baseline topographic characteristics for different landscapes. Significant differences in topographic characteristics between landscapes in different natural regions and between high disturbance landscapes and all other landscapes illustrates that the methodology and identified metrics may be used to guide reclamation design and evaluation by industry and regulators, but the range of values in topographic characteristics are specific to natural areas and different levels of disturbance.