Assessing Optimal Digital Elevation Model Selection for Active River Area Delineation Across Broad Regions

The Active River Area (ARA) is a spatial approach for identifying the extent of functional riparian area. Given known limitations in terms of input elevation data quality, ARA studies to date have not achieved effective computer-based ARA components delineation, limiting the efficacy of the ARA framework in terms of informing riparian conservation and management. To determine the optimal input elevation data for future ARA studies, this study tested a novel digital elevation model (DEM) smoothing algorithm and assessed ARA outputs derived from a range of DEMs for accuracy and efficiency. It was found that the tested DEM smoothing algorithm allows the ARA framework to take advantage of high-resolution LiDAR DEM and considerably improves the accuracy of high-resolution LiDAR DEM derived ARA results; smoothed LiDAR DEM in 5-m spatial resolution best balanced ARA accuracy and data processing efficiency and is ultimately recommended for future ARA delineations across large regions. The scientific findings provided by this study further enhance the efficacy of the ARA framework, and ultimately the confidence in modelled ARA outputs for application in riparian conservation and management contexts across broad geographic regions.


Introduction
Riparian areas are the interface between aquatic and terrestrial systems along inland watercourses (Naiman and Décamps 1997). They cover a small proportion of the landscape mosaic but play important roles in both terrestrial and aquatic ecosystems. Movement of water and organic materials within riparian areas makes freshwater ecosystems one of the most biologically diverse and rich systems (Reid et al. 2019). Riparian areas are also among the most degraded ecosystems due to human disturbances (Cao and Natuhara 2020). Given the important ecological role of riparian areas and the degradation caused by human disturbances, effective riparian management is warranted (Singh et al. 2021). The initial stage of riparian management is to delineate the functional riparian area, which is important for both ecological and managerial reasons (De Sosa et al. 2018). By understanding the extent of riparian area and how rivers interact with adjacent lands, managers can better recognize the processes involved and how to design effective protection and restoration measures.
Delineating the functional riparian area extent, however, is not a straightforward task and may require the understanding of local-scale environmental conditions since riparian functional attributes depend on site-specific variations such as landscape composition and environmental setting. Given these challenges, riparian management projects have typically applied a "fixed-width buffer guideline" (Kuglerová et al. 2020;Singh et al. 2021). Unfortunately, fixed-width buffers may ultimately fail in meeting management goals as they do not account for variations in site-specific attributes and likely underestimate the full extent of riparian areas (Kuglerová et al. 2014).
The Active River Area (ARA) framework (Smith et al. 2008) was developed to address challenges associated with shortcomings in standard buffer-width guidelines. To facilitate delineation of key functional riparian areas, the framework separates the area into five components-floodplain, terrace, meander belt, riparian wetland, and material contribution zone-based on the valley setting and geomorphic stream type. The most extensive application of the ARA framework to date was collaboratively conducted by The Nature Conservancy (TNC) and the Nature Conservancy of Canada (NCC) for the Northern Appalachian/Acadian Ecoregion, based on 30-m Shuttle Radar Topography Mission (SRTM) DEM and national scale hydrologic network (1:100,000 U.S. National Hydrography Dataset and 1:50,000 Canadian National Hydrographic Network). Although laudable for applying ARA delineation across a large geographic area, its limitations are acknowledged in accuracy of component outputs associated with the input DEM quality. Failure to differentiate terraces from modern floodplains arises from the vertical imprecision of SRTM DEMs and limits the utility of the results for riparian management purposes (Smith et al. 2008). Refinements of methods and input elevation data are therefore required for more precise yet technically feasible ARA application in large geographies.
Theoretically, the spatial extent of ARA is heavily controlled by stream and valley morphology. The use of DEMs and digitized hydrological networks should achieve effective ARA delineation by synthesizing river geomorphology into a spatial framework. Previous studies reveal that DEM horizontal resolution and vertical accuracy play important roles in topographic and hydrologic modelling (Goyal et al. 2018;Roostaee and Deng 2020). The accuracy of ARA delineation therefore relies heavily on the properties of the input DEM. Generally, high quality DEMs, such as Light Detection and Ranging (LiDAR) derived DEMs, provide a detailed representation of topography (Zhao et al. 2010;Wu and Lane 2017), which can benefit river valley morphology modelling. High resolution DEMs, however, have their own limitations and are not necessarily better than lower resolution DEMs in terms of topographic representation and hydrologic modelling (Tan et al. 2018). A very high-resolution LiDAR DEM, for instance, could result in the representation of a terrain surface that is much more detailed than necessary for the process being modelled, potentially creating errors in hydrologic modelling and ARA delineation, and imposing challenges with increased data storage and computational intensity requirements for processing and manipulation (Lin et al. 2010).
Studies focused on selecting optimal DEMs that satisfy both data processing and storage capabilities and representation of spatial variability for ARA modelling are necessary. To date, efforts have been made to identify optimal input DEMs for several topographically based flood-inundation models (e.g., TOPMODEL, Bathtub, HAZUS-MH). For instance, Tate et al. (2015) indicated that DEM accuracy greatly influences the uncertainty of HAZUS-MH flood modeling, based on global sensitivity analysis. Fereshtehpour and Karamouz (2018) used both deterministic and probabilistic approaches to assess the performance of various resampled LiDAR DEMs in terms of flood inundation modelling and suggested that resampled LiDAR DEM is more appropriate than publicly available DEMs for Bathtub-based flood inundation analysis. However, no ARA-related studies have carried out accuracy-efficiency analyses to identify optimal input DEM, limiting the efficacy of the ARA framework and confidence in modelled outputs in riparian management contexts. Moreover, the advantages of high-resolution LiDAR DEMs in terms of topographic and hydrologic modelling have been hindered by excessive surface roughness resulting from over-detailed topographic representation (Van Nieuwenhuizen et al. 2021). To allow the ARA framework to take advantage of high-resolution LiDAR DEMs, an appropriate DEM smoothing algorithm must be applied to subdue surface roughness before commencing ARA delineation.
Accordingly, the primary objectives of this study were to achieve ARA methodology improvement and identify optimal input DEM for future ARA delineation across large regions. To achieve these objectives, we (i) tested whether the application of high-resolution LiDAR DEM can achieve floodplain and terrace separation; (ii) delineated the ARA using the refined framework and various DEM scenarios; (iii) assessed the performance of various DEM scenarios on ARA delineation; (iv) identified and tested an appropriate DEM smoothing algorithm to reduce the surface roughness of high-resolution LiDAR DEM; (v) conducted error analyses for DEMs of varying resolutions and sources and analyzed DEM-induced topographic and hydrologic modelling uncertainties; and, (vi) applied an accuracy-efficiency tradeoff analysis to select the optimal DEM for future ARA studies, especially in large regions.

Study Watersheds
To delimit the scope while maintaining a hydrologically meaningful boundary, two watersheds-Lower St. John River Watershed and Miramichi River Basin (Fig. 1)were selected from the Canadian Hydrological Unit (CHU) system at the CHU-6 watershed level (Supplementary Information [SI] 1). The watersheds are located within the province of New Brunswick (NB), Canada, where full 1-m LiDAR DEM coverage is available. The Lower St. John River Watershed and Miramichi River Basin represent topographical variation within the study region, with elevational ranges above mean sea level of − 5-553 m and − 4-763 m, respectively; and comprise areas of similar sizes (i.e., 14,343 km 2 and 12,069 km 2 , respectively). The main stems are the Lower St. John River, which is fed by many significant tributaries and flows into the Bay of Fundy, and the Miramichi River, which comprises two main branches, and drains into Miramichi Bay in the Gulf of St. Lawrence.

Data Acquisition and Preparation
Input datasets included 1-m LiDAR DEM, 30-m SRTM DEM, NB Hydrographic Network (NBHN), provincial forest soils data, and satellite derived flood extent (SI 2). To produce multi-resolution LiDAR DEM for test purposes, 1-m LiDAR DEM was resampled into five lower resolution DEMs (i.e., 3-m, 5-m, 10-m, 15-m, 30-m) using Bilinear interpolation (bilinear interpolation justification provided in SI 3). The NBHN is a digital hydrological network with a general horizontal accuracy of ± 2.5 m. Strahler stream order (Strahler 1957) was applied to separate watercourses into five size classes (i.e., Great, Medium, Small Rivers, Creeks and Headwaters) based on their Strahler order numbers for further ARA delineation (size class justification provided in SI 4).

Active River Area Delineation
To improve the efficacy of the ARA framework, methodological refinements were made for separating floodplains and terraces and delineating material contribution zones. Given that environmental conditions vary along meandering watercourses and are difficult to predict without the support of an interdisciplinary research team, this study did not quantify the extent of meander belts. To test the performance of different DEMs on the remaining four components of ARA delineation, the refined framework was first applied using various input DEMs. The spatial extents of resulting ARA outputs were then compared to illustrate differences generated by variations in DEM scenarios.

Floodplain and Terrace Delineation
To delineate the extent of riparian basezones (i.e., floodplains and terraces), watercourse grids were treated as source cells from which the least accumulated cost distance for every output cell location was calculated. Friction surfaces were generated from slope grids (in degrees) of different DEMs. Satellite-derived flood extent was then applied to calibrate cost thresholds for each river size class to delineate riparian basezones (SI 5). Once riparian basezones were determined, efforts were made to separate terraces from floodplains (SI 6). The proposed separation framework aims to detect the cost distance value underneath the transition lines between floodplain and terrace in slope and planform over high-resolution LiDAR DEM surfaces (i.e., 3 m and 5 m). Locations of transition lines were determined by analyzing river valley profiles extracted from LiDAR DEMs. Based on the location of the first level terrace, the floodplain-terrace transition line can be identified over the cost distance surface, and thereafter the underlying cost distance value of the transition line can be extracted. To get a representative cost threshold for terraces and floodplains separation, several cross-sectional lines were created for different rivers in the same size class; the average of the extracted cost distance values was then calculated and used as the separation threshold for that river size class.

Riparian Wetland Delineation
Riparian wetlands are low-gradient areas with inundated or hydric soils from high surface water and/or groundwater levels. To detect land areas with relevant features, topography and flow accumulation patterns were modelled separately. Topography patterns were modeled by analyzing slope changes across the watersheds; flow accumulation patterns were modelled by creating flow accumulation grids (Wu and Lane 2017;Stengård et al. 2020). 'Topography Moisture Index' (TMI) (SI 7) was then calculated based on the following formula.
A high TMI value indicates a high flow accumulation value and a low slope value, and therefore areas with high TMI values may be considered potential wetlands. The extent of riparian wetland was then determined by applying a wet-dry threshold calibrated by provincial wetland and hydric soil data.

Material Contribution Zone (MCZ) Delineation
Based on a review of buffer sizes for riparian conservation regulation, MCZ width was determined as 60 m. However, this was found to be inadequate for delineating MCZs for rivers located in confined valleys with steep slopes that extend beyond 60 m (Fig. SI 8A). If slopes adjacent to streams have a 15% or greater rise, sediments and nutrients will enter from areas beyond 60 m to the top of the slope. To address this, a component named 'ARA steep slope area' (i.e., from stream to top of ≥ 15% slope) was delineated to supplement the 60-m-buffer MCZ (SI 8).
The ARA components then were combined based on their flooding likelihood to delineate the final ARA extents (the detailed description is provided in SI 9).

High-Spatial Resolution LiDAR DEM Smoothing
Surface roughness in high-resolution DEM surfaces has been found to inevitably affect the characterization of macro-scale topography (Lindsay et al. 2019), ultimately affecting the accuracy of ARA results. DEM smoothing algorithms, commonly low-pass filters (i.e., mean, median and Gaussian filters), are frequently applied to LiDAR DEMs to lessen the surface roughness before commencing topographic and hydrologic modelling (Wu and Lane 2017;Stengård et al. 2020). These techniques, however, tend to blur edges of important topographic and drainage features (e.g., ditches, gullies) (Van Nieuwenhuizen et al. 2021). The ability of high-resolution LiDAR DEMs to represent these local-scale features provides the justification for their application in many studies (Zhao et al. 2010;Tan et al. 2018). Thus, it is counter-productive to smooth high-resolution DEMs if the outcome is the removal of important topographic information.
The DEM smoothing algorithm applied in this study is the "Feature-Preserving DEM Smoothing" algorithm (FPDEMS) (Lindsay et al. 2019), specifically designed to smooth DEMs while preserving important topographic features. The smoothing process was automatically carried out by the 'FeaturePreservingSmoothing' tool, available in the open-source geospatial analysis platform, Whitebox GAT. Three mail tool parameters were used-11 * 11 kernel, 15° threshold, and 3 iterations for 5-m LiDAR DEM; and more aggressively, 17 * 17 kernel, 25° threshold, and 5 iterations for 3-m LiDAR DEM-representing a compromise between processing time and performance (Lindsay et al. 2019).

DEM Error and Uncertainty Assessment
Three types of DEM errors were considered, arising from resampling, data sources and smoothing. Errors were estimated by computing Root Mean Square Error (RMSE) and Mean Error (ME) for each DEM with respect to the reference elevation dataset (Eqs. 2, 3), where 'Zt' refers to the elevation value for the 'ith' point obtained from the reference dataset, 'Zp' is the elevation value for the 'ith' point obtained from the lower quality DEMs, and 'n' is the total number of points for which elevation values were retrieved. The reference dataset comprised elevations obtained from the 1-m LiDAR DEM. To calculate the RMSE and ME, points (i.e., 1000) were randomly generated and evenly distributed across landscape types (e.g., toe slope, mid-slope, shoulder slope) as ground control points (GCPs) in each watershed. Slopes and TMIs derived from different DEMs were selected and applied to analyze the effect of DEM qualities on topographic and hydrologic modelling, as these are major indices for riparian basezone and riparian wetland delineation, respectively. (

Computer-Based Accuracy-Efficiency Tradeoff Analysis
Tradeoff analyses were conducted by separately measuring accuracy and efficiency components of different parallel ARA analyses. The measurement of ARA output accuracy involved the derivation of Kappa Coefficients based on comparison of delineated and reference ARA components for a set of specific locations. Given the absence of field verification data and high-resolution aerial imagery, and the potential error of 3-m LiDAR DEM derived ARA caused by over-detailed topographic representation, the reference dataset used was the 3-m smoothed DEM-derived ARA. The 1000 GCPs were then applied to extract the underlying ARA components from the reference DEM and other DEM-derived ARAs for confusion matrix creation. For efficiency, computational effort was expressed as time spent to run key geoprocessing tools. Among various geoprocessing tools applied for ARA delineation, "Cost Distance" and "Flow Accumulation" are the most time-consuming. Runtimes for these two geoprocessing tools considerably increase with higher DEM resolutions. Thus, runtimes of these two tools were considered as efficiency components and applied for comparison against accuracy assessment results.

Variations in Spatial Extent of Different DEM-Derived ARAs
Riparian basezone is the dominant component of the active river area in all spatial resolutions and sources of DEM-derived ARA results (Table 1, Figs. 2, 3). The proportion of riparian basezone located on the "wet" area is greater than that on the "non-wet" area in all DEM-derived ARA results (Table 1). This is not surprising, since riparian basezones are generally low-gradient areas with hydric soils, having the same features as riparian wetlands. However, riparian basezones derived from different LiDAR-DEMs vary in their spatial extent (Table 1, Figs. 2,3). A negative relationship between the LiDAR-DEM spatial resolution and the total area of the riparian basezone can be detected in both study watersheds. For instance, the total basezone areas (wet and non-wet) in the Miramichi River Basin range from 1711.66 km 2 (30-m resolution) to 3238.17 km 2 (3-m resolution) (Table 1b). Additionally, while similar boundary extents were detected for riparian basezones of great, medium, and small rivers, more blank areas (i.e., no riparian basezone delineated) were found within the extents of high-resolution LiDAR DEM derived riparian basezones (Fig. 2). By visually inspecting the aerial photo of those blank areas, it was evident that most of the blank areas are human infrastructural features with sharp boundaries. Because high-resolution LiDAR DEMs can generate more detailed slope grids and model subtle topographic transitions, human infrastructure with sharp boundaries can be represented and will receive high-cost distance values, which will ultimately exclude them from the extent of riparian basezone. Accordingly, it is reasonable to suspect that low-resolution LiDAR DEMs will overpredict the riparian basezone extent at certain areas (e.g., human encroachment areas).
LiDAR DEMs in 3-m and 5-m spatial resolution were found to provide sufficiently detailed elevational and horizontal distance changes in river valleys to allow differentiation of terraces from floodplains for great, medium and small rivers. The floodplain and terrace separation results show that the floodplain is the dominant component within the riparian basezone for rivers in those size classes (Fig. SI 10A). Similar to floodplain distribution, larger proportions of terraces are located on "wet" areas than "non-wet" areas in both study watersheds (Table 1). However, due to the appearance of excessive local surface roughness (Fig. SI 10B), accurate valley morphology information for small streams (i.e., creeks and headwaters) was not able to be extracted, as causes (i.e., presence of terraces or surface roughness) of abrupt slope changes were indistinguishable. By visually comparing the riparian wetland extents derived from different LiDAR-DEMs, it was found that the extents of high-resolution riparian wetland match better than those of low-resolution with the provincial wetland inventory (Fig. SI 10C). Lowresolution LiDAR DEMs and SRTM DEM based analyses were found to underpredict the extent of riparian wetland and overlapping areas of riparian basezone and riparian wetland (Table 1). Such overlapping areas are considered important in conservation planning because highly productive and diverse riparian plant communities tend to establish themselves in areas with rich alluvial soils (Reid et al. 2019). The underprediction issue in low quality DEM-based analyses can therefore be problematic in riparian conservation applications.
Steep slopes are concentrated in relatively small portions of both study watersheds (Fig. SI 10D). The area of steep slope MCZs increases as the resolution of LiDAR DEM becomes higher, with the same trend in both watersheds (Table 1). For example, the area of steep slope MCZs increases from 671.7 km 2 (30-m LiDAR) to 1003.6 km 2 (3-m LiDAR) in the Lower St. John River Watershed. This result is not surprising, since high-resolution LiDAR DEM can generate more detailed slope grid and pick up subtle topographic changes on bare Earth, while low-resolution LiDAR DEMs and SRTM DEM will eliminate subtle topographic transitions when converting into slope grids. Given this limitation, only extremely steep slope areas will be identified from low-resolution slope grids, while

High-Spatial Resolution LiDAR DEM Smoothing
The FPDEMS method effectively subdued topographic complexity at local scales in both 3-m and 5-m LiDAR DEMs, while retaining the complexity of macro-scale landforms (Fig.  SI 10E). The rough appearance of raw high-resolution LiDAR DEMs was successfully removed, while the boundaries of important topographic features (i.e., channel edges) were preserved (Fig. SI 10E c, d). Moreover, the run time of the FPDEMS method is acceptable-~ 5 min and 25 min for 5-m LiDAR DEM and 3-m LiDAR DEM, respectively, in both study watersheds-with the aid of a high-end laptop.
Smoothed high-resolution LiDAR DEMs generated more realistic creek and headwater basezones in comparison to raw LiDAR DEMs (Fig. SI 10F). Raw high-resolution LiDAR DEMs exhibited many cut-off planforms or so-called "padi terraces" at the local scale (Fig.  SI 10E). These areas are typical of closed contours where all the surrounding pixels exhibit the same slope value. Accordingly, all pixels inside the closed contours were assigned the same cost distance value and thus the determined cost thresholds did not necessarily cut off accurate basezone extents for creeks and headwaters. Through visual inspection, the widths of creeks and headwaters modelled by raw high-resolution LiDAR DEMs appear much greater than the general widths of creek and headwater floodplains, reported as ranging from 5 to 50 m (Benda et al. 2007). Furthermore, many creek and headwater basezones are connected to each other, thereby magnifying the problem of overprediction (Figs. 2,  3). The FPDEMS method successfully removed these cut-off planforms; pixels therefore received accurate cost distance values given their location to source cells. As a result, the determined cost distance thresholds effectively cut off areas that are not likely to be dynamically linked to creeks or headwaters (Fig. SI 10F).
The area of floodplain derived from smoothed DEM was considerably smaller than those derived from raw LiDAR DEMs. For instance, in the Lower St. John River Watershed, total floodplain areas (wet and non-wet) derived from raw and smoothed 5-m LiDAR DEMs are approximately 3130.16 km 2 and 1552.98 km 2 , respectively, a decrease of 1577.18 km 2 (> 50%) ( Table 1). The area of riparian wetland shows a similar trend in both study watersheds: in the Miramichi River Basin, the area of riparian wetland decreases from 1200.37 to 889.88 km 2 (25%) with smoothing of 5-m LiDAR DEM (Table 1). This result is not surprising, since the FPDEMS algorithm removed many cut-off planforms (i.e., areas of low slope likely to be misidentified as riparian wetlands). An increasing trend in area of steep slope with smoothing can be found in both watersheds (Table 1), also explained by removal of cut-off planforms.

DEM Error and Uncertainty
For all resampled and smoothed LiDAR DEMs, MEs were in the centimeter level, while MEs calculated for SRTM DEMs were considerably greater than those for LiDAR DEMs of the same spatial resolution, indicating that biases introduced by DEM resampling and smoothing were negligible as compared to those introduced by DEM source (Table SI 10A). RMSEs increased with decreasing resolution of LiDAR DEMs, indicating lower accuracy of lower-as compared to higher-resolution LiDAR DEMs, with the same trend found in both study watersheds. For instance, in the Lower St. John River Watershed, RMSEs range from 0.2588 to 1.0281 m, increasing for all resampled LiDAR DEMs as the resolution decreases from 3-m to 30-m (Table SI 10A). Additionally, RMSEs for SRTM 30-m DEMs are significantly higher than for resampled 30-m LiDAR DEMs, by approximately 4 times in both study watersheds, indicating overall low accuracy of SRTM DEMs as compared to LiDAR DEMs of the same resolution. It is also interesting to note that RMSEs of smoothed high-resolution DEMs and raw high-resolution DEMs differ only marginally (Table SI 10A), indicating the ability of the FPDEMS method to preserve original elevation values.

Effects of DEM Error on Topographic Modelling
Maximum slope values increase dramatically as the resolutions of resampled LiDAR DEMs increase (or RMSE decrease) (Fig. SI 10G). For example, maximum slope value increases from approximately 60° to 80° as LiDAR DEM resolution increases from 30-m (RMSE = 1.0281) to 3-m (RMSE = 0.2588) in the Lower St. John River Watershed (Fig. SI  10G). Because smoothed and raw high-resolution LiDAR DEMs have similar RMSE, the ranges of their slope values are similar (Fig. SI 10G), reflecting an ideal characteristic of the FPDEMS in terms of important topographic feature preservation. The inner value distributions of slope grids, however, differ between smoothed and raw high-resolution LiDAR DEMs (Fig. SI 10G). For instance, in both study watersheds, interquartile ranges (i.e., 25th to 75th quartiles) of derived slopes from smoothed high-resolution LiDAR DEMs are smaller than those from raw high-resolution LiDAR DEMs, as the FPDEMS method removed several planforms with low-slope values.

Effects of DEM Error and Uncertainty on Hydrologic Modelling
TMIs derived from DEMs with low RMSEs tend to follow a normal distribution; in contrast, TMIs derived from DEMs with high RMSEs are positively skewed, though it is difficult to comment on the accuracy of one distribution over another in the absence of reference data. In the Miramichi River Basin, for instance, the difference between median (744) and mean (747) values of 5-m LiDAR DEM (RMSE = 0.3550) derived TMI is 3; differences between mean and median TMI value increase to 162.32 and 217.89, respectively, for the TMIs derived from 30-m LiDAR DEM (RMSE = 0.8854) and 30-m SRTM DEM (RMSE = 3.8066) ( Table  SI 10B). A declining trend can be seen in both mean and median TMIs as DEM RMSEs increase; a similar trend can be found in both study watersheds (Table SI 10B), indicating that more pixels received lower TMI values in the low-quality DEM-based parallel analyses, possibly explaining in part the underprediction issue of riparian wetland. Additionally, the DEM smoothing algorithm lowered the mean and median values of TMI in both study watersheds, which partially explains the declining trend in the area of riparian wetland shown in Table 1, as more pixels received lower TMI values and were considered dry areas.

Accuracy Assessment
The results of accuracy assessment for different DEM-derived ARAs, including user's accuracy, producer's accuracy, and Kappa Coefficients, are presented in SI 11. A power trendline can best explain the relationship between Kappa Coefficient and RMSE, with the correlation of determination values (R 2 ) at 0.76 and 0.93 in the Lower St. John River Watershed and Miramichi River Basin, respectively (Fig. 4a, b). The relationship between RMSEs and Kappa Coefficients of raw high-resolution LiDAR-DEM derived ARAs, however, cannot be explained by this power trendline (Fig. 4a, b), as the creek and headwater basezone extents are grossly inaccurate due to over-detailed topographic representation (Fig. SI 10F). In addition to the reference ARA, the 5-m smoothed LiDAR DEM-derived ARAs achieved highest overall accuracies in final ARA output, with Kappa Coefficients near 0.8 in both study watersheds (Table SI 10C). Since the SRTM-DEM has the highest RMSE, the Kappa Coefficients of 30-m SRTM-DEM derived ARAs are smaller than those of 30-m LiDAR DEM derived ARAs in both study watersheds, reflecting the weak ability of SRTM-DEM in terms of ARA delineation as compared to LiDAR DEM of the same resolution.

Efficiency Assessment
Using consistent computational power (i.e., laptop with 6-core 2.2 GHz processor and 24 GB of memory), there are moderate differences in total runtime for key geoprocessing tools among low-resolution LiDAR DEM-based analyses (Table SI 10D). In contrast, a significant increase in runtime is evident once LiDAR DEM resolution reaches 5-m. Accordingly, the 5-m LiDAR DEM can be assumed as a turning point in runtime Fig. 4 The relationship between RMSE and Kappa Coefficient (a, b) and relationship between runtime and Kappa Coefficient of different ARA parallel analyses (c, d). Raw 3-m and 5-m LiDAR DEM derived ARAs (triangle) were excluded from the correlation of determination calculation to reduce bias, given raw 3-m and 5-m LiDAR DEM-derived ARAs contain considerable errors in creek and headwater basezones in both study watersheds. The overall accuracy of 5-m LiDAR DEM derived ARA, however, is relatively low due to the inaccurate delineation of creek and headwater basezones. It is also interesting to note that the DEM smoothing algorithm saves ~ 30% in data processing time by smoothing noisy microtopographic detail. For instance, in the Miramichi River Basin, total runtime decreases from 591 to 414 min for raw 3-m and smoothed 3-m DEM, respectively (Table SI 10D).

Accuracy-Efficiency Trade-Off
Accuracy-efficiency curve for different ARA parallel analyses was created using runtime and Kappa Coefficient as variables. An exponential trendline can best fit the data (Fig. 4c, d). The inflection point exists at the 5-m smoothed LiDAR DEM, with Kappa Coefficient at approximately 0.8 while data processing time remains low (Fig. 4c, d). As indicated by Cohen (1960), Kappa Coefficient greater than "0.8" represents perfect agreement between delineation and observation (i.e., 3-m smoothed DEM-derived ARA). The 5-m smoothed LiDAR DEM therefore achieves high accuracy in ARA results. Although 3-m smoothed LiDAR DEMs can achieve the most accurate ARA results in both watersheds, runtimes to process 3-m and 5-m smoothed LiDAR DEMs differ greatly (by almost 3 times), indicating that 3-m smoothed LiDAR DEM is less efficient than 5-m smoothed LiDAR DEM. To balance accuracy and efficiency, 5-m smoothed LiDAR DEMs are recommended for future ARA delineations, especially across large spatial extents or multiple watersheds.
This recommendation should be applied carefully, however, since the accuracy-efficiency trade-off analysis is not without limitations. Although it has been well acknowledged in the literature that high-resolution LiDAR DEMs can achieve higher accuracy in results in terms of floodplain delineation (Zhao et al. 2010;Tan et al. 2018), there is no reference literature to demonstrate that ARA studies also fit this trend. The reliability of the reference ARA needs further examination; and, the results of trade-off analyses may differ if reference ARAs change. Nevertheless, the accuracy assessment framework provides a unique lens through which to effectively assess ARA accuracy. In future studies, ground reference surveys are recommended to verify and enhance the reliability of the reference ARA. Ground reference surveys can enhance the accuracy of RMSEs calculated for different DEMs by replacing the reference elevation values extracted from 1-m LiDAR DEM with elevation values measured in the field.
The RMSE was applied to indicate DEM quality in this study, which assumes error to be aspatial. Error related to the DEM, however, varies spatially and cannot be sufficiently assessed by a global metric such as RMSE (Hawker et al. 2018). Indeed, scientists have observed that local DEM error can be large and spatially correlated, though global RMSE is small (Bater and Coops 2009), and thus it is possible that the optimal DEM resolution for ARA delineation is landscape specific. More comprehensive ways to consider DEM error from a spatial perspective (or landscape perspective), such as Classification and Regression Tree (CART) analysis (Bater and Coops 2009) and Monte Carlo-based Sequential Gaussian Simulation (SGS) analysis (Fereshtehpour and Karamouz 2018), are therefore required. Such analyses can allow the user to better understand the relationship between DEM quality and accuracy of ARA. Optimal DEM for ARA delineation may also vary with differences in landscape features, and thus, ideally, or eventually, it is anticipated that optimal DEM will be identified for specific landscape contexts in future applications.

Conclusion
This study (1) tested whether high-resolution LiDAR DEM can achieve floodplain and terrace separation; (2) verified the applicability of a novel DEM smoothing algorithm on ARA delineation; (3) analyzed the relationship between DEM qualities and topographic and hydrologic modelling uncertainties; and (4) identified optimal input DEM for future ARA delineation across large regions by assessing various DEMs for ARA output accuracy and data processing efficiency. Results show that a refined ARA framework based on high-resolution LiDAR DEM was successful at floodplain and terrace separation, representing an important methodological improvement. Comparison of ARA outputs from various input DEMs reveals differences in performance of DEMs on ARA delineation, suggesting that high-resolution LiDAR-DEM can more precisely model the extent of riparian basezone and steep slope areas and detect areas of overlap between riparian basezone and riparian wetland, which are considered to be of conservation value. Analyzed topographic and hydrologic modelling uncertainties were shown to be subject to input DEM quality, suggesting that high-quality input DEM can reduce the topographic and hydrologic modelling uncertainties and enhance ARA output accuracy. The selected FPDEMS algorithm was found to be successful at smoothing the high-resolution LiDAR DEM surface roughness caused by over-detailed topographic representation at local spatial scales, while preserving topographic and hydrologic features important to accurate ARA delineation. Finally, the accuracy-efficiency tradeoff analysis suggests that 5-m smoothed LiDAR DEM represents the optimal DEM for future ARA studies, especially those applied across large study sites.
The findings address important gaps and needed refinements in ARA methodologies, enhancing the efficacy of the ARA framework, and ultimately the confidence in modelled ARA outputs for application in riparian conservation and management contexts across broad geographies. Important questions remain about the reliability of the reference ARA applied for the accuracy assessment and the identification of landscape-specific optimal input DEMs. Further analyses are needed to enhance the reliability of the reference dataset for accuracy assessment and identify optimal input DEM for specific landscape contexts. Nonetheless, the study provides two novel and important methodological contributions: (1) it describes and tests a refined approach to ARA delineations based on high-resolution LiDAR DEM that is successful at floodplain and terrace separation; and, (2) it demonstrates through an accuracy-efficiency tradeoff analysis that 5-m smoothed LiDAR DEM represents an optimal DEM for ARA analyses applied across multiple watersheds or large study sites. Researchers and practitioners concerned with watershed and riparian conservation and management can benefit from the scientific findings of this study to more effectively and efficiently identify functional riparian area extent at the initial stage of riparian conservation projects.