Remote sensing-based framework for detecting and interpreting permafrost terrain hydrologic connectivity

Abstract

Thermal shifts have accelerated ice-wedge degradation and reorganized polygonal trough networks along the Utqiaġvik, Alaska coastline. Thus, quantifying their structural variability and hydrologic connectivity across spatial scales remains challenging. This study applies high-resolution remote sensing and terrain analysis to evaluate hydrologic controls on ice-wedge polygon morphology. Using a 0.5 m lidar-derived DEM, polygon boundaries were manually digitized and compared with Thiessen (Voronoi) tessellations to quantify structural divergence, boundary misalignment, and differences in area. Hydrologic influences on polygon development were assessed through compound terrain analysis, drainage network extraction, and surface flow modeling. Spatial intersection analyses reveal geometric discrepancies underscoring the limitations of automated polygon proxies in permafrost terrain. Modeled flow paths exhibit spatial congruence with mapped trough networks, indicating that surface hydrology plays a role in trough evolution. Hydrologic and morphometric parameters demonstrate high runoff potential associated with low-relief topography. Elevated Topographic Wetness Index (TWI > 12 ) and Stream Power Index (SPI > 60) values delineate zones of saturation and flow accumulation that frequently coincide with trough depressions. Despite high predicted runoff, minimal gradients produce ponded flow regimes that promote wetland formation. This integrative framework enhances detection and interpretation of permafrost terrain features, providing a scalable methodology for monitoring Arctic landscape dynamics.

1 Introduction

Permafrost landscapes occupy approximately 24% of the Northern Hemisphere’s land surface and are undergoing rapid transformation due to freeze–thaw cycles and associated changes in hydrologic and geomorphic processes (Brown et al., 1997; Romanovsky et al., 2010; Biskaborn et al., 2019). Among the most visually and functionally distinctive features throughout these landscapes are ice-wedge polygons—regularly patterned ground structures formed through thermal contraction cracking and subsequent infilling by ice and sediment (Leffingwell, 1915; Mackay, 1990). These polygons dominate Arctic lowlands and coastal plains, exerting significant influence on surface hydrology, microtopography, vegetation, and soil thermal regimes (French, 2018; Liljedahl et al., 2016). However, as permafrost thaws and surface subsidence progress, the morphology and hydrologic connectivity of polygonal terrain are increasingly altered, leading to shifts in landscape drainage, stability, and carbon cycling (Jorgenson et al., 2006; Lara et al., 2015; Nitzbon et al., 2020; Speetjens et al., 2024; Freitas et al., 2025).

Understanding the structural organization and hydrologic function of ice-wedge polygons is critical for assessing permafrost vulnerability and predicting the trajectory of Arctic landscape change. Polygon networks are composed of elevated rims, depressed centers, and troughs formed above ice wedges. These features control the direction and intensity of surface water flow, influencing thaw driven erosion, thermokarst formation, and local hydrologic feedbacks (Mackay, 2000; Kanevskiy, et al., 2017; Shur et al., 2025). Subtle microtopographic variations, often less than 1 m in relief, govern the accumulation and movement of meltwater and precipitation, creating small scale hydrologic gradients that feed back to surface stability and thaw patterns. Consequently, accurately representing these features through high-resolution terrain and hydrologic modeling is essential for quantifying permafrost degradation and geomorphic change.

Traditional field-based mapping and photo interpretation have provided foundational insights into polygon development and degradation (French, 2018; Black, 1976). However, recent advances in remote sensing technologies, particularly light detection and ranging (lidar), have enabled unprecedented detail in capturing Arctic microtopography. High-resolution digital elevation models (DEMs) derived from airborne or terrestrial lidar now facilitate fine-scale detection of polygon rims, troughs, and centers, as well as derivation of secondary topographic parameters such as slope, curvature, and aspect (Boike et al., 2013; Jones et al., 2013; Abolt and Young, 2020). When coupled with hydrologic terrain modeling, these data provide a robust means to evaluate how surface water dynamics interact with polygonal morphology to drive permafrost change (Lara et al., 2015; Liljedahl et al., 2012; Nitzbon et al., 2019; Abolt et al., 2020; Harp et al., 2021).

Recent work has further emphasized that fine-scale hydrologic connectivity plays a critical role in accelerating permafrost degradation and broader Arctic carbon feedbacks. Studies have shown that hydrologic redistribution within degrading polygonal terrain can enhance thaw settlement, alter drainage efficiency, increase dissolved organic carbon export, and accelerate greenhouse gas emissions from Arctic systems (Abolt et al., 2020; Harp et al., 2021; Speetjens et al., 2024; Freitas et al., 2025). Recent studies have also highlighted how permafrost mass wasting, thermokarst expansion, and hydrologically driven landscape instability are becoming increasingly important across ice-rich Arctic terrain (Young et al., 2026), while broader reviews of polygon evolution emphasize the role of changing surface hydrology in reshaping ice-wedge systems (Shur et al., 2025). These findings detail the growing need for high-resolution approaches capable of identifying preferential flow pathways and hydrologic connectivity patterns that may serve as early indicators of accelerating permafrost degradation and carbon cycling feedbacks.

Despite these advances, significant challenges remain in accurately delineating and characterizing polygonal networks. Manual delineation based on visual interpretation of DEMs can capture complex polygon geometries and degraded structures but is labor intensive and subject to interpreter bias. In contrast, automated methods such as Thiessen (Voronoi) polygon generation provide a computationally efficient means of partitioning space into geometrically regular units (Thiessen, 1911; Okabe et al., 2000). However, such idealized representations often fail to reproduce the irregular and hydrologically dynamic nature of real-world polygonal landscapes. The differences between these methods raise critical questions about how spatial simplification affects the representation of terrain morphology and its hydrologic implications.

Hydrologic modeling provides a complementary perspective by explicitly linking terrain geometry to water flow processes. Indices such as the Topographic Wetness Index (TWI) and Stream Power Index (SPI) (Beven and Kirkby, 1979; Moore et al., 1991) integrate slope and contributing area to represent potential saturation and erosive energy, respectively. These indices, when applied to high-resolution DEMs, can reveal fine scale zones of water accumulation, drainage convergence, and surface runoff intensity that directly correspond to polygon troughs and degradation zones (Liljedahl et al., 2016; Zhang et al., 2003). Such coupling of geomorphic and hydrologic modeling enables a process-based understanding of how surface water redistribution both reflects and reinforces the structural evolution of ice-wedge polygon terrain (Liljedahl et al., 2024).

Given the increasing rates of permafrost thaw and surface hydrologic redistributions observed across Arctic lowlands, there is a pressing need to integrate geometric and hydrologic approaches to better characterize the morphologic and functional evolution of polygonal landscapes. Previous studies have emphasized that microtopographic processes such as rim collapse, trough widening, and thermokarst ponding are strongly linked to surface and subsurface hydrologic pathways (McCarthy, 1994; Kanevskiy et al., 2017; Young et al., 2026). Yet quantitative comparisons of polygon delineation methods and their implications for hydrologic modeling remain limited, particularly at sub-meter scales.

This study addresses this gap by combining manual delineation and Thiessen polygon generation with hydrologic terrain modeling derived from a 0.5 m lidar DEM to analyze the structure and function of an Arctic ice-wedge polygon network. Specifically, we (1) compare the geometric properties and spatial overlap of manually delineated and Thiessen-derived polygons, (2) evaluate terrain derivatives and hydrologic indices that describe microtopographic and flow dynamics, and (3) examine the relationship between polygon morphology and modeled surface hydrology to infer processes of thaw and degradation. By integrating these approaches, this work aims to advance understanding of how hydrologic processes are both influenced by and contribute to the structural development of polygonal permafrost, offering new insights into the geomorphic and hydrologic evolution of Arctic lowland systems. The testable hypothesis is that hydrologic processes are spatially and functionally aligned with ice-wedge trough networks, such that areas of high modeled flow accumulation, topographic wetness, and erosive potential coincide with mapped trough depressions. Through the integrated approach, this study will enhance the detection and interpretation of permafrost terrain dynamics and provide a scalable framework for assessing Arctic hydrologic connectivity and permafrost vulnerability.

2 Study area

Utqiaġvik, Alaska—formerly known as Barrow—is located on the Arctic Ocean coast within Alaska’s North Slope Borough at approximately 71.3°N latitude and 156.8°W longitude (Figure 1). As the northernmost city in the United States, it lies above the Arctic Circle and is bordered to the northwest by the Chukchi Sea. The region is characterized by continuous permafrost extending to depths of approximately 400 m (Kanevskiy et al., 2017) and a low-relief landscape dominated by thermokarst features such as lakes, ponds, and drained thaw-lake basins. The surface geology consists of unconsolidated marine, eolian, and lacustrine sediments, including shales, clays, siltstones, and sandstones (McCarthy, 1994). Peat-rich soils with a high organic content overlay these deposits, supporting a tundra ecosystem dominated by wetland vegetation, such as sedges, grasses, mosses, and lichens. The local climate is strongly influenced by proximity to the Beaufort and Chukchi Seas. Average temperatures range from −19°F to −7°F during the winter months (November to April), and from 36°F to 47°F during summer (June to September). Synoptic-scale geostrophic winds contribute to relatively consistent surface wind speeds throughout the year (Moritz, 1979). The soil active layer, which thaws seasonally from late May through August, separates the ground surface from the underlying permafrost and plays a critical role in hydrologic and geomorphic processes in the region. Due to the proximity of the permafrost, this frozen sediment acts as a shallow confining layer limiting vertical water movement; thus, water remains at the surface with limited infiltration into the subsurface (McCarthy, 1994). This water forms lakes and ponds, as well as stream channels, due to the limited relief of the drainage basins. Across the landscape, drainage basins are determined by topographic highs and lows that contribute to the formation of ice-wedge polygon troughs that alter drainage patterns. Such surface drainage phenomena, including the formation of polygonal troughs and differential thawing of the active layer, are often pronounced and may induce corresponding alterations in subsurface drainage pathways and consequently in groundwater flow directions. Additionally, the limited vertical extent of the active layer constrains groundwater flow to shallow regimes, with flow systems operating over spatial scales typically ranging from centimeters to several tens of meters rather than across broader, regional domains (McCarthy, 1994).

FIGURE 1

3 Methodology

This study used an integrated remote sensing and geomorphometric approach to analyze the structure and hydrologic dynamics of ice-wedge polygon networks. High-resolution (0.5 m) Unoccupied Aerial System (UAS) lidar-derived DEMs were acquired using the RIEGL VUX-12023 3D lidar sensor manufactured by RIEGL Laser Measurement Systems GmbH, United States of America Inc. headquartered in North America in LAS format on 17 July 2025, during snow free peak summer conditions, providing detailed microtopography suitable for examining spatial structure and surface flow connectivity. Ice-wedge troughs were first manually delineated from hillshade and terrain derivatives, then compared to idealized Thiessen polygons generated from manually selected polygon intersection points. Terrain and hydrologic analyses—including flow direction, flow accumulation, modeled channel networks, and Topographic Wetness Index—were used to assess how hydrologic routing and terrain indices relate to trough morphology and network organization (Figure 2).

FIGURE 2

3.1 Data acquisition

Lidar data was acquired using a RIEGL VUX-12023 sensor, a high-performance airborne laser scanning system designed for high-resolution topographic mapping and capable of capturing fine-scale surface variability in low-relief landscapes (RIEGL Laser Measurement Systems GmbH, 2021; Glennie et al., 2013; Toth and Jóźków, 2016). The sensor provides high point density (∼450 pt/m²) and centimeter-scale vertical accuracy to resolve subtle microtopographic features characteristic of ice-wedge polygon terrain.

The lidar system was integrated with a high-precision GPS/IMU (Applanix APS-20 manufactured by Trimble Applanix headquartered in Canada), and differential post-processing was performed using a Trimble base station to ensure accurate georeferencing of the point cloud. Raw lidar returns were recorded in LAS format and subsequently geometrically corrected to account for sensor orientation, platform motion, and terrain induced distortions. These processing steps produced a spatially consistent, high-resolution point cloud suitable for DEM generation and hydrologic terrain analysis.

3.2 Data processing

CloudCompare v2.12.4 (64-bit) was used as the primary platform for processing and analyzing the 3D point cloud data. Each LAS dataset was evaluated using key metadata attributes, including acquisition timestamp, number of returns, scan angle, scan direction, intensity, point source ID, and flight line information. However, data fields with default or missing values were excluded from further analysis. Global shift and scale values were retained and referenced to the original coordinate system at the time of data acquisition for each survey site.

To clean and classify the raw point cloud data, the Cloth Simulation Filter (CSF) algorithm was applied to extract ground points from the discrete return lidar datasets (Zhang et al., 2016; Cai et al., 2023). The CSF algorithm is widely utilized in lidar post-processing workflows and has been validated in numerous studies (Cai et al., 2023; Saritas and Kaplan, 2023; Yilmaz, 2021; Bailey et al., 2022; Sabirova et al., 2019). For this study, a cloth resolution of 0.5 m was used to define the grid size, and the filter was run for a maximum of 500 iterations. Point classification of ground and non-ground categories was based on the vertical distance between each point and the simulated cloth surface. A classification threshold ranging from 0.1 to 0.8 was iteratively tested to determine the optimal value for each site. The classification threshold range was determined empirically through exploratory sensitivity testing of terrain and hydrologic derivatives derived from the 0.5 m DEM. Threshold values were incrementally varied to evaluate their influence on feature continuity, noise suppression, and the delineation of polygon troughs and interiors. Threshold values were incrementally varied and evaluated using three criteria: (1) preservation of continuous polygon trough networks, (2) suppression of isolated topographic artifacts and high-frequency noise, and (3) consistency between derived terrain features and visually interpreted polygon morphology from the raw point cloud. Lower thresholds (<0.1) resulted in excessive noise and fragmented features, while higher thresholds (>0.8) overly generalized the landscape and excluded coherent trough networks. In cases where these optimization targets conflicted, priority was given to preserving coherent hydrologic connectivity and recognizable polygon morphology while minimizing noise related artifacts. The selected range therefore brackets the values over which hydrologically meaningful and spatially stable patterns were consistently observed and was used to assess relative feature sensitivity rather than to define a single optimal classification boundary.

Following ground point extraction, the lidar tool by Headwall Photonics, Inc application was employed to generate digital elevation products. GPS and IMU sensor offsets were accounted for, including a roll offset of 2.2°, a pitch offset of −9°, and a yaw offset of 0°. After positional correction and filtering, digital surface models (DSMs) were produced. The process was subsequently repeated to generate digital elevation models (DEMs), using a grid spacing of 0.1 m and a mean output value for elevation interpolation.

3.3 Polygon delineation

To characterize the spatial structure of ice-wedge polygons, two delineation methods were employed: manual digitization based on visual interpretation of high-resolution topographic data, and computational approximation using Thiessen polygon generation. This dual approach was used to examine differences in geometric representation and to compare idealized polygonal partitions with morphologically expressed polygon boundaries derived from surface microtopography. The two polygon datasets are not treated as independent or hierarchical in accuracy, but rather as contrasting representations of polygon network structure for comparative geometric analysis.

3.3.1 Manual digitization of ice-wedge polygons

Ice-wedge polygons were manually delineated using lidar-derived DEMs and associated topographic derivatives, including hillshade, slope, and curvature rasters (Figure 3a). Digitization was conducted in QGIS at a consistent mapping scale to ensure positional consistency and geometric accuracy across the study domain. Ice-wedge polygon boundaries were delineated through a rule-based interpretation of terrain and hydrologic surfaces. Polygon edges were defined primarily by continuous trough features identified using a combination of relative elevation, slope, curvature, and flow accumulation layers. Boundaries were traced along coherent linear depressions corresponding to trough centers and margins, while polygon interiors were constrained by convex microtopographic high and low connectivity.

FIGURE 3

Polygon boundaries were identified through visual interpretation of surface microtopography, with troughs recognized as continuous or semi-continuous concave linear depressions associated with ice-wedge networks and hydrologic connectivity. Although field validation was not available, delineated polygon boundaries were indirectly validated through spatial alignment with independently derived hydrologic flow networks and terrain convergence metrics. Because independent field validation data were unavailable, manually delineated polygons are interpreted as terrain-informed reference features rather than ground-truth boundaries.

To reduce bias, delineation followed consistent decision criteria applied uniformly across the study area, including: (1) alignment with persistent trough depressions visible across multiple terrain derivatives, (2) continuity of flow accumulation pathways, and (3) consistency with known geometric characteristics of ice-wedge polygon networks. Polygons were only delineated where boundaries were clearly expressed, and ambiguous features were excluded.

Manual delineation allows incorporation of geomorphic context, such as irregular trough geometry, asymmetric degradation, and locally discontinuous rims, that are difficult to parameterize using fully automated geometric methods. The resulting vector geometries represent an interpretive, terrain-constrained depiction of polygon boundaries rather than a validated ground-truth dataset. No field-based validation or inter-operator consistency assessment (e.g., Kappa statistics) was conducted; therefore, manual delineations are not assumed to be more accurate than automated representations but are used as a morphologically informed reference for evaluating systematic geometric differences among polygon models.

Although manual interpretation inherently involves expert judgment, subjectivity was minimized by using high-resolution lidar data and standardized delineation rules rather than visual interpretation of a single surface. The resulting polygons therefore represent conservative, reproducible representations of well-defined ice-wedge structures. Subsequent comparison with automated Thiessen polygon generation provides a complimentary geometric benchmark for evaluating systematic differences between idealized and terrain-informed polygon representations.

3.3.2 Thiessen polygon generation

To provide an idealized geometric reference for ice-wedge polygonal spatial organization, Thiessen (Voronoi) delineated vectorized polygons were generated using manually identified intersection points as seed locations (Figure 3b). Seed points were placed at inferred polygon centers where three or more trough vertices converged, based on consistent topographic indicators observed in the lidar-derived digital elevation model, including concave curvature patterns and local elevation minima.

The Thiessen polygon method partitions space into non-overlapping regions such that each polygon contains locations closest to its associated seed point under the assumption of isotropic spatial distribution (Thiessen, 1911; Okabe et al., 2000). In this study, Thiessen polygons are not treated as an independent or objective control dataset. Because both Thiessen polygons and manually delineated terrain-based polygons involve human interpretation, strict independence is not claimed.

Instead, Thiessen polygons represent an idealized geometric partitioning of space and serve as a conceptual baseline for comparison with terrain constrained polygonal boundaries derived from lidar microtopography. The comparison is therefore intended to evaluate systematic geometric divergence between idealized and morphologically expressed structures rather than to serve as an accuracy assessment (Equations 1–3).

Formally, given a set of seed points S = {s₁, s₂, … , s_n}, where each s_i= (x_i, y_i), the Voronoi polygon V_i associated with s_i is defined as the locus of all points x = (x, y) that satisfywhere the Euclidean distance function d (·) is

The boundaries between adjacent polygons correspond to the set of points equidistant to s_i and s_j, satisfyingwhich defines the perpendicular bisector between the two input points (Okabe et al., 2000; Aurenhammer, 1991).

3.4 Hydrologic terrain modeling

Hydrologic and terrain modeling was performed with the lidar-derived DEM as the primary input. The analysis focused on deriving terrain attributes and hydrologic indices to evaluate surface water flow dynamics and their relationship with ice-wedge polygon morphology. Four core processes were conducted: DEM preprocessing, terrain parameter extraction, channel network, and hydrologic flow modeling.

The study area is underlain by continuous permafrost exceeding several hundred meters in depth, with a shallow active layer that thaws seasonally (Woo, 2012). Accordingly, vertical percolation was assumed negligible, and permafrost was treated as an effective aquiclude. Surface flow was therefore restricted to the land surface and shallow active layer, and subsurface routing was not represented. The hydrologic analysis is thus intended to characterize relative surface connectivity and flow organization rather than to simulate full permafrost hydrologic processes. Under these conditions, surface topography exerts first-order control on lateral water redistribution, making DEM-based flow routing a suitable proxy for assessing hydrologic organization in polygonal tundra (Woo, 2012; Bui et al., 2020).

3.4.1 DEM pre-processing

The DEM was hydrologically conditioned using the Fill Sinks (Wang & Liu) algorithm, which enforces a minimum downslope gradient while minimizing distortion of surrounding elevations (Wang and Liu, 2006). This approach is well suited for high-resolution DEMs with complex microtopography because it preserves realistic flow paths compared to traditional sink filling. After sink removal, aGaussian smoothing filter with a 6 m kernel radius and 2σ standard deviation was applied to reduce high-frequency noise while retaining coherent polygon scale features (Zhang et al., 2003). The kernel radius reflects a horizontal smoothing scale smaller than the typical spacing between polygon rims and trough, ensuring that meter scale elevation contrasts remain intact. These pre-processing steps improve slope, curvature, and flow direction calculations by stabilizing gradient surfaces while maintaining meaningful microtopographic patterns. Although smoothing slightly attenuates sharp local gradients, terrain derivatives are interpreted based on relative spatial variation. All procedures were performed in floating point precision to minimize rounding errors and ensure numerical stability in subsequent hydrologic modeling.

3.4.2 Terrain parameter extraction

The terrain analysis computed geomorphometric parameters based on spatial derivatives of the pre-processed DEM (Equations 4–10). The following key terrain metrics were extracted:

Slope (S): Calculated as the gradient magnitude using finite difference approximations of the first-order partial derivatives of elevation in the x and y directions (Wang and Liu, 2006; Zhang et al., 2003).where (rate of change in elevation in the x (east–west) direction) and (rate of change in elevation in the y (north–south) direction) are the first-order partial derivatives of elevation z with respect to spatial coordinates x and y.

Aspect (A): Azimuthal direction of maximum slope, derived from the arctangent of the ratio of elevation gradients along orthogonal axes (Burrough, et al., 2015; Wilson and Gallant, 2000).

Adjusted to compass directions (0°–360°) based on quadrant.

Plan and Profile Curvature: Second-order derivatives quantifying the convexity or concavity of the terrain surface in horizontal (plan) and vertical (profile) planes, computed via local polynomial fitting of elevation surfaces (Zevenbergen and Thorne, 1987; Shary et al., 2002).

Curvatures describe how the surface bends locally, based on second derivatives:

This equation assumes true derivatives p, q, r, s, t, as they are approximated from elevation values. Due to the derivatives already containing the necessary length factors, the length factor is defaulted to the grid discretization. As a result, curvature values are scale consistent regardless of DEM resolution.

Plan Curvature (curvature perpendicular to the slope direction):

Profile Curvature (curvature in the slope direction):

The exponents m and n in Equation 10 are empirical parameters derived from experimental observations of soil loss under varying slope lengths and gradients (Moore and Burch, 1986; Wischmeier and Smith, 1978). Constant m controls sensitivity to slope length (flow accumulation or contributing area). Constant n controls sensitivity to slope steepness (slope gradient effect).

3.4.3 Channel network and drainage basin

To characterize surface drainage patterns and analyse trough connectivity within ice-wedge polygon networks, a compound channel network and drainage basin approach was performed. This approach integrated topographic metrics to delineate drainage features from a high-resolution DEM (Equations 11–20). As a result, surface flow pathways based on geomorphic and hydrologic indicators were produced. The following outputs were computed:

This is defined within a geometric domain, specifically the Euclidean planar domain (ℝ²). This domain defines a scalar distance field over the domain, where each raster cell is assigned the shortest horizontal (planimetric) distance to the nearest channel cell.

This provides an estimate of gravitational potential energy available for flow or erosion.

Flow Accumulation: Quantifies the number of upstream cells or areas that contribute flow to a given cell based on the terrain defined flow direction. In this assessment, the algorithm used was the Multiple Flow Direction (MFD). This algorithm allows the fractional distribution of flow from a grid cell to multiple downslope neighbors, rather than assigning all flow to the single steepest direction, as in Deterministic 8 (D8). The D8 algorithm assigns flow direction to one of eight neighboring cells based solely on the steepest downslope gradient; this approach can exaggerate channelization and grid alignment in low-relief terrain, where flow divergence and diffuse redistribution are common.

The MFD algorithm was therefore selected to better represent dominant diffuse flow patterns across low-slope polygon interiors characteristic of ice-wedge polygon landscapes (Quinn et al., 1991; Seibert and McGlynn, 2007). It is recognized, however, that MFD algorithms may underestimate flow concentration along locally steep microtopographic features, such as trough margins, where slopes can reach several degrees. Accordingly, the use of MFD represents a scale dependent modeling compromise intended to capture dominant diffuse surface flow organization across heterogeneous polygonal terrain rather than resolve concentrated flow at individual trough edges.

Only neighbors where s_i > 0 are considered.

The value β = 1.3 comes from empirical studies focused on modeling hydrology (Quinn et al., 1991). For values greater than 1, flow is often concentrated along steeper paths but is not forced entirely along a single steep slope. This is ideal for low-relief terrains, where flow spreads out gradually but tends to follow subtle slope gradients. This value avoids over diffuse flow, which would occur with β = 1. Additionally, unrealistic channelization would be avoided with β = 2.0 − plus.

The value starts at 1 to count the cell itself and then accumulates flow from contributing neighbors.

The base level is the elevation at the outlet or nearest downstream channel cell.

Drainage Basin: Derived from flow direction and accumulation modeling. Each basin B_j is defined as the set of all cells that drain to a specific channel mouth (O'Callaghan and Mark, 1984; Jenson and Domingue, 1988; Freeman, 1991):

This approach allows distributed flow accumulation and provides a better representation of low-relief, convergent landscapes compared to single direction methods. The MFD approach was selected because the study area is characterized by predominantly low-relief polygonal tundra and surface flow commonly dispersed across broad polygon interiors before converging into trough networks. The localized trough margins with steeper gradients are spatially limited and occur primarily along narrow trough edges rather than across the broader landscape. MFD was therefore considered appropriate for capturing dominant lateral flow partitioning across polygon networks while preserving flow convergence into downslope trough systems. In contrast, single-flow algorithms such as D8 may artificially over-concentrate flow paths across low-gradient landscapes.

Channel Heads and Mouth: These parameters were determined by analyzing the flow paths and topography. The Channel Head is identified by a cell on the channel network with no inflowing upstream channel cells, while the Channel Mouth is a channel cell at the lowest elevation with no downstream channel cell (Peckham, 2009).

3.4.4 Hydrologic indices

The analysis performed used a high-resolution DEM to derive hydrologic indices related to surface water flow and terrain driven saturation potential. Total and specific catchment areas were computed using the Multiple Flow Direction algorithm (Quinn et al., 1991). The Topographic Wetness Index (TWI) and Stream Power Index (SPI) were calculated from catchment and slope metrics, following the formulas derived from (Beven and Kirkby, 1979; Moore et al., 1991). Closed depressions were extracted to identify local surface storage zones, and a calibrated flow accumulation threshold was used to define channel initiation points. These hydrologic terrain metrics were used to assess trough connectivity, flow routing, and hydrologic controls on polygon evolution (Equations 21–24).

Closed Depressions: Depressions were identified as local minima in the elevation surface that lacked a downslope flow path. These features were detected post-processing using a sink detection algorithm applied after pit removal (Wang and Liu, 2006). Cells within a depression were identified where no downslope neighbor existed under the flow routing algorithm (MFD).

Total Catchment Area (TCA): Computed using the MFD algorithm (Quinn et al., 1991), which distributes surface flow to downslope neighbors based on relative slope weights. For each grid cell, TCA represents the total contributing area from all upslope cells. TCA is computed the same as flow accumulation; however, TCA represents the total contributing area from upslope terrain, while flow accumulation is the number of upslope cells draining to a given cell (Freeman, 1991). These metrics are related but are not interchangeable and are reported separately throughout this study to avoid unit inconsistency. In knowing the grid cell size of 0.5 m resolution for the given DEM, the conversion is written as:

This metric improves the hydrologic nature by accounting for flow concentration per unit width, especially in microtopographic depressions and troughs.

TWI values increase where large catchment areas coincide with gentle slopes, representing potential surface saturation zones such as ice-wedge polygon troughs.

Higher SPI values indicate areas of concentrated flow and greater downslope energy, often associated with incised troughs or flow concentration zones.

Channel Initiation Threshold: Channel networks were delineated by applying a threshold to the flow accumulation raster using TCA (Montgomery and Dietrich, 1989). A minimum contributing area (A_min) was empirically determined by iteratively testing thresholds against known trough locations:

The channel initiation threshold was set to 25 m², corresponding to 100 contributing cells in the 0.5 m resolution DEM. This threshold was selected empirically to optimize correspondence between modeled flow paths and observed trough networks in the ice-wedge polygon terrain. The 25 m² threshold provided the best balance between eliminating noise and preserving the observed hydrologic connectivity. Through iterative sensitivity testing across lower and higher thresholds and evaluated against mapped trough geometry, thresholds below 25 m² produced fragmented and noise driven flow paths across polygon interiors, while thresholds above 25 m² excluded primary trough segments observed in the DEM. The selected threshold aligns with the characteristic spacing and contributing geometry of primary polygon troughs in the study area, where flow typically converges over short distances before entering laterally connected trough networks.

3.5 Polygon structure and hydrologic comparison

To evaluate structural divergence and hydrologic correspondence between computationally derived and terrain-informed polygon representations, a two-part comparative analysis was conducted: (1) geometric comparison between manually digitized and Thiessen polygons, and (2) hydrologic alignment assessment examining the spatial correspondence between polygon networks, terrain derivatives, modeled flow pathways, and hydrologic connectivity.

3.5.1 Spatial overlay analysis

Geometric comparison between manually delineated ice-wedge polygons and Thiessen-derived polygons was performed using spatial overlay and intersection analyses. Each Thiessen polygon was intersected with its corresponding manually mapped counterpart to quantify differences in area, perimeter, and boundary alignment. Statistical metrics were used to evaluate structural congruence and geometric variability across the study area. Manual polygons generally captured irregular and non-equidistant geometries representative of natural microtopographic variability, whereas Thiessen polygons provided an idealized, isotropic framework with uniform vertex spacing. Discrepancies in alignment and area were interpreted as indicators of natural deformation and thermokarst driven variability within the polygon network.

3.5.2 Hydrologic alignment evaluation

To assess qualitative hydrologic correspondence of polygon delineations, spatial correspondence between polygon trough networks and modeled hydrologic pathways was evaluated. The degree of alignment between polygon trough networks and modeled hydrologic pathways was then evaluated to produce qualitative evidence of spatial correspondence between the two approaches. Served as a proxy for assessing the role of surface hydrology in polygon evolution. Strong spatial correspondence between modeled hydrologic pathways and trough depressions indicated that surface water routing exerts a primary control on ice-wedge degradation and trough expansion.

3.5.3 Data interpretation and relationship analysis

The integrated geometric and hydrologic analyses were used to evaluate relationships between surface water dynamics and ice-wedge polygon morphology. Spatial comparisons among hydrologic indices were conducted to assess how hydrologic connectivity corresponded with polygon structure and potential zones of geomorphic change. Areas displaying elevated TWI and flow accumulation values were examined in relation to trough intersection, polygon boundaries, and closed depressions to identify patterns of surface water concentration and potential subsidence pathways. Conversely, areas with lower hydrologic connectivity were evaluated to characterize regions with limited surface drainage and reduced potential for flow development.

4 Results

Hydrologic modeling based on the principles of gravitational flow, water accumulation, and erosion potential was employed to simulate the routing of surface water across the landscape. Key indices, such as the TWI and SPI, were used to identify zones of saturation and erosion, providing insights into the thermal and mechanical processes driving the thawing and evolution of ice-wedge polygons. The results reveal strong correlations between topographic features, hydrologic flow paths, and polygon characteristics, underscoring the role of hydrological processes in shaping the physical structure of the landscape and influencing the degradation of permafrost.

4.1 Ice-wedge polygon delineation

Manually delineated ice-wedge polygons (Figure 3a) were compared with Thiessen (Voronoi) polygons derived from intersection points (Figure 3b) to evaluate representational differences, geometric assumptions, and the effects of simplifying natural terrain into idealized spatial partitions. This analysis assesses relative geometric correspondence rather than formal accuracy, with manually mapped polygons treated as terrain informed references.

A total of 1515 manual polygons and 1509 Thiessen polygons were analyzed (Table 1). Thiessen polygons covered nearly twice the total area (226,236 m² vs. 116,419 m²) and had larger mean and median areas (149.92 m² and 96 m²) compared to manual polygons (76.84 m² and 39 m²), indicating Thiessen polygons exhibited consistently larger mean areas and different variability patterns relative to manually delineated polygon features, reflecting systematic geometric divergence between idealized and terrain informed representations. Shape metrics also differed as Thiessen polygons were more compact and regular (mean Shape Index 0.774), whereas manual polygons showed broader morphological diversity (mean 0.725), reflecting natural variability in ice-wedge troughs.

The percent overlap and geometric agreement were calculated to quantify how much of a given polygon area was retained in the intersection between the manually delineated and Thiessen-derived polygons. Results indicated a strong agreement for manual polygons, with a mean overlap of 88% and a median of 100%, indicating that most were fully contained within their intersected Thiessen polygons (Table 2). In contrast, Thiessen polygons had lower and more variable overlap (mean 68%, median 74%), with some showing 0% overlap, revealing clear spatial mismatches or overgeneralizations, especially in areas of degraded polygon morphologies. These results highlight the limitations of idealized Thiessen partitions for capturing trough alignment, hydrologic pathways, and degraded polygon morphology.

4.2 Terrain parameter analysis

Terrain derivatives were extracted to characterize microtopographic variability to inform hydrologic and geomorphic interpretations of ice-wedge polygon networks (Figure 4). Elevation values extracted from the DEM ranged from 0 to 3.3 m (0–4 m scale) above mean sea level (Figure 4a). The landscape exhibits minimal vertical relief, characteristic of Arctic coastal plains dominated by permafrost and polygonal tundra. Despite the narrow elevation range, subtle microtopographic variations (often <1 m) were sufficient to delineate polygon rims, centers, and troughs, as well as to influence surface water routing and pond formation. The DEM displayed elevated values highlighting polygonal thermokarst heaved features.

FIGURE 4

Aspect was calculated in radians and converted to degrees (0°–359°) for directional analysis (Figure 4b). While the overall aspect distribution appeared isotropic at the landscape scale, directional clustering at the microtopographic scale indicated trough alignment trends along the southwest and southeast axes, with a mean aspect value of 181°, suggesting a dominant trend of south facing slopes across the study area.

Slope values ranged from 0 to a maximum of 7°. The landscape is characterized by low-relief terrain, consistent with ice-wedge polygon morphology (Horn, 2005) (Figure 4c). Flatter zones (0°–2°) generally corresponded to polygon centers of ponded trough intersections, while localized increases in slope (up to 7°) were observed along polygon rims, trough margins, and incised microchannels. These subtle yet significant slope gradients influence surface flow accumulation and may promote differential thaw and erosion along trough networks. The slope gradients were enhanced visually by producing a hillshade generated using a 315° azimuth and 45° altitude to enhance visualization of subtle topographic features (Figure 4d). This visualization layer improved manual delineation accuracy of polygon boundaries and highlighted microtopographic heterogeneity (Figures 4e–h).

Plan and Profile curvature values were derived to understand microtopographic formation and flow behavior in the permafrost landscape (Figure 5). Plan curvature values ranged from −0.312 (convergent) to +0.341 (divergent), highlighting the alternating structure of polygonal troughs and rims. Convergent features corresponded closely with mapped trough networks and modeled flow paths, indicating zones of potential water accumulation and thaw. Divergent features, such as rims or elevated centers where surface water tends to spread outward, were also captured. Profile curvature ranged from −0.312 (concave) to +0.341 (convex), with concave slopes generally aligned with troughs, promoting flow acceleration and potential ice-wedge degradation. Negative (concave) values indicate downslope facing troughs or channel heads that contribute to increased surface water flow or thawing potentials. Positive (convex) values were observed predominantly along polygon centers, rims or around a heaved surface, where the slope begins to transition. These values suggest a symmetric distribution of curvature in the landscape, which is common in polygonal terrain, alternating convex polygon centers/rims and concave troughs.

FIGURE 5

To enhance the interpretation of surface routing, thaw dynamics, and geomorphic differentiation within the polygon network relative slope position (RSP), closed depressions, and LS-Factor terrain derivatives were applied (Figure 6). RSP values, which represent the relative elevation of a cell along a slope profile (ranging from 0 = low-lying tundra to 1 = ridge top), revealed a clear stratification of polygon components (Figure 6a). The spatial distribution of RSP also correlated with areas of surface flow convergence, indicating potential pathways for lateral water movement during thaw and precipitation events.

FIGURE 6

Closed depressions were extracted to identify hydrologically isolated basins where water may accumulate, leading to surface ponding, thermokarst formation, or localized ground ice melt (Figure 6b). These depressions, primarily occupying polygon centers or trough networks, were typically shallow (depth <0.6 m) but significant in extent due to the high spatial resolution. Their distribution aligned with low RSP values and regions of minimal slope, indicating a correlation between microtopographic lows and permafrost vulnerability.

The LS-Factor, integrates both slope length and steepness to estimate relative erosional susceptibility. Although the study area is characterized by very low local slope gradients, resulting in generally modest LS-Factor values, localized LS-Factor peaks were detected along polygon rims, incised trough networks, and channelized flow paths where long, laterally continuous flow paths and flow convergence increase the slope-length component of the LS formulation (Figure 6c). Under these low-relief conditions, LS-Factor variability is therefore driven primarily by contributing area and flow path length rather than slope angle alone. Values ranged from 0 (low) to 5 (very high), reflecting spatial variability in flow accumulation and microtopographic convergence. Areas of high LS-Factor often occurred with negative plan curvature and concave profile curvature, reinforcing the geomorphic linkage between flow convergence, enhance runoff concentration, and thaw related surface degradation.

To assess how microtopography controls surface hydrology and polygon organization, cumulative distribution functions (CDFs) of TWI, LS-Factor, and RSP were evaluated (Figure 7). TWI shows a pronounced rightward shift in closed depressions relative to depression terrain, with a ponding to runoff transition around TWI ≈ 6-7 (Figure 7a). Median TWI is higher in closed depressions (6.49 vs. 5.83), indicating greater contributing area and enhanced saturation potential in enclosed polygon centers and troughs. A cutoff of CD = 0.05 distinguishes hydrologically meaningful depressions from effectively non-depressional areas, consistent with threshold-based filtering approaches in similar studies (Huang et al., 2019; Jiang et al., 2023).

FIGURE 7

The LS-Factor CDF indicates generally low erosional potential, with a median of ∼0.03 and only a small portion of the landscape showing elevated values (>1.5) that correspond to continuous troughs and concentrated flow paths (Figure 7b). RSP distributions reveal strong topographic stratification where ∼75% of the surface has RSP < ∼0.13, characteristic of troughs and low-centered polygons, while values >0.6 delineate rims and elevated features (Figure 7c). These thresholds were empirically identified from inflection points in the CDF and were further supported through visual comparison with mapped polygon morphology and terrain derivatives. Values between these thresholds represent transitional slope positions between trough bottoms and elevated polygon rims. The shaded threshold zones were selected based on distribution inflection points, percentile breaks, and consistency with observed polygon morphology rather than arbitrary thresholding. Collectively, these distributions show that hydrologic connectivity and geomorphic sensitivity in this low-relief terrain are governed by subtle variations in relative elevation, flow accumulation, and slope position rather than slope magnitude alone.

4.3 Channel network analysis

Surface hydrology emerged as a key control on evolution of ice-wedge polygons by influencing trough development, thaw dynamics, trough deepening or widening, and lateral water redistribution. Polygon troughs often act as micro-channels that concentrate and direct overland flow, creating pathways that may enhance thermal erosion and modifying surface connectivity. Although regional gradients are minimal, the terrain exhibits high runoff potential once local surface storage thresholds are exceeded. Shallow micro-depressions and polygon centers frequently promoted localized surface ponding; however, limited infiltration due to the shallow active layer and underlying permafrost facilitated rapid lateral spillover into polygon troughs, which function as efficient surface drainage pathways. Hydrologic terrain metrics revealed these patterns through distinct spatial distributions of horizontal and vertical overland flow distance, flow accumulation, and channel network distance, which collectively highlighted highly connected drainage pathways across the polygonal tundra landscape (Figure 8).

FIGURE 8

Horizontal overland flow distance was derived from the DEM to quantify the planimetric distance that surface water travels before converging into a defined flow path (Figure 8a). To better capture microtopographic variability across the ice-wedge polygon landscape, the analysis scale was refined to a range of 0–16 m. This adjustment improved sensitivity to subtle changes in surface drainage pathways within low-relief terrain. Results indicate that areas with shorter horizontal flow distances (0–4 m) correspond closely with mapped polygonal troughs and networks, where minimal topographic resistance facilitates rapid flow convergence. In contrast, longer distances (12–16 m) were primarily observed in the interior of larger, elevated, or less degraded polygons, where flow is more dispersed and hydrologic connectivity is limited.

Vertical overland flow distance was calculated to quantify the cumulative elevation change that surface water experiences along its downslope flow path (Figure 8b). This metric provides a measure of the vertical component of surface runoff potential and is particularly useful for assessing subtle elevation gradients in low-relief permafrost regimes. Initial output values ranged from 0 to 0.9799 m; however, to enhance the detection of microtopographic variation relevant to polygon trough morphology, the analysis was rescaled to a range of 0–0.22 m. The adjusted result reveals that low vertical flow distances (0–0.05 m) are concentrated along well-defined troughs and topographic depressions, indicating minimal elevation change and efficient lateral water movement. Higher vertical flow distances (0.15–0.22 m) occur predominantly within elevated polygon interiors or poorly connected polygons, where flow paths traverse broader elevation gradients before entering drainage features.

To further quantify surface runoff efficiency across the polygonal tundra, horizontal and vertical overland flow distances were jointly evaluated and colored by RSP (Figure 9a). The results reveal a monotonic relationship between planimetric flow distance and cumulative elevation loss, reflecting the constrained routing behavior of surface water in low-relief terrain. Low RSP values (≤0.2), corresponding to polygon troughs and microtopographic lows, are characterized by short horizontal travel distances (<4 m) and minimal vertical change (<0.05 m), indicating rapid lateral convergence into drainage pathways. In contrast, higher RSP values (>0.6), associated with polygon rims and elevated interiors, exhibit longer horizontal flow distances (up to ∼16 m) and greater cumulative elevation loss (up to ∼0.22 m), reflecting more dispersed overland flow prior to channel entry.

FIGURE 9

Flow accumulation was derived to quantify the number of upslope contributing cells draining to each grid cell, effectively modeling potential surface runoff concentration across the ice-wedge polygon landscape (Figure 8c). Total catchment area was subsequently computed by converting flow accumulation into contributing area (m²) using the 0.5 m grid resolution. Output values ranged from 0 to 2413 cells (equivalent to 603 m² TCA). Higher accumulation values were concentrated along well-developed troughs and drainage networks, indicating zones of hydrologic convergence where surface water is most likely to concentrate under gravitational flow.

High flow accumulation values (>1000) with an approximate contributing area of 250–603 m² were aligned with the deepest and most continuous polygon troughs, particularly at intersections where multiple flow paths converged. These locations represent primary drainage conduits across the polygon network and suggest persistent hydrologic activity that may enhance thermal erosion and trough deepening over time. Moderate values (250–1000) with an approximate contributing area of 62.5–250 m² were distributed along secondary troughs and shallow microchannels, while low values (<250) of approximately 0–62.5 m² were more concentrated near polygon interiors and elevated rims, indicating minimal upslope contribution and limited hydrologic connectivity.

The joint density of flow accumulation and RSP further illustrates the strong control of microtopography on channeled flow development (Figure 9b). High cell densities are connected at low RSP values (<∼0.2), indicating that the largest contributing areas are preferentially associated with polygon troughs and low-centered terrain. Flow accumulation decreases systematically with increasing RSP and elevated polygon rims (RSP > ∼0.6) exhibiting consistently low accumulation values. This pattern demonstrates that, in this low-relief permafrost terrain, surface drainage organization is governed primarily by relative elevation position, with subtle microtopographic gradients exerting control on flow convergence and channel formation.

Channel network distance was calculated to quantify the relative distance of each grid cell from the modeled flow network (Figure 8d). This metric provides insight into the lateral extent and symmetry of surface hydrologic influence relative to polygon troughs. Output values ranged from −0.52 to 0.77 m, with zero values corresponding to cells located directly along the modeled troughs. Negative values represent distances to the left bank of the drainage channel, while positive values correspond to the right bank of the drainage channel.

To isolate zones of active trough influence and reduce the effect of broader terrain variability, the analysis was rescaled to a focus range of −0.08 to 0.2 m, which approximately represents a 30 cm corridor around each channel. Within this refined range, values close to 0 indicate strong alignment with the trough network. Cells with values between −0.08 and −0.01 m were located on the left banks of the troughs, within approximately 8 cm of the modeled channel, while those between 0.01 and 0.2 m fell on the right banks, extending up to 20 cm away. This 28 cm wide corridor around the troughs represents zones of greatest hydrologic convergence and microtopographic control.

The distribution of channel network distance values highlights consistent trough alignment with modeled flow paths across the landscape, indicating that surface drainage processes are tightly coupled with polygon structure. Minor asymmetries in the spread of positive and negative values were observed in degraded or transitional polygons, suggesting localized thaw-induced surface deformation.

4.4 Hydrologic indices and polygon interaction

To further evaluate surface hydrology and its interaction with ice-wedge polygon networks, a drainage basin analysis was performed. Individual drainage basins were delineated based on flow direction and accumulation, enabling the spatial partitioning of the landscape into discrete runoff contributing units. Within each basin, hydrologic indices were computed, including total catchment area, specific catchment area, Topographic Wetness Index (TWI), Stream Power Index (SPI), and channel initiation threshold (Figures 10, 11). These metrics collectively characterize the distribution and intensity of runoff, the potential for water accumulation, and the erosive capacity of surface flow. When analyzed in conjunction with mapped polygon networks, the results provide a detailed assessment of how microtopographic and hydrologic processes influence trough formation, connectivity, and permafrost vulnerability across the polygonal landscape.

FIGURE 10

FIGURE 11

TCA was derived from flow accumulation by converting contributing cell counts to area using the 0.5 m cell size (Figure 10a). Using the DEM, TCA values ranged from 0 to 87,341 grid cells, equivalent to 0 to 21,835.25 m2 of the contributing area. This indicates that some drainage basins in the study area receive runoff from approximately >20,000 m2 of upslope terrain, likely corresponding to well-developed trough networks or larger basin outlets. High TCA values were concentrated in major trough intersections and low-lying drainage corridors, indicating extensive upslope connectivity and runoff convergence. Conversely, low values occurred across polygon interiors and elevated microtopography, reflecting limited contributing area and minimal hydrologic input. These patterns align with surface flow paths modeled in the channel network analysis and highlight zones of increased saturation and potential vulnerability to thaw related degradation.

SCA was created from the total catchment area and slope to identify zones of potential water accumulation across the landscape (Figure 10b). The maximum value of 173,260 corresponds to over 43,315 m2 per unit contour width, indicating low slope and high accumulation zones likely to experience surface ponding and prolonged saturation. Additionally, the high SCA zones were concentrated in broad trough intersections and near surface depressions. Low values were identified on steeper microtopographic features such as polygon rims or heaved ground between troughs.

TWI was developed using the natural logarithm of specific catchment area divided by local slope (Figure 10c). These values ranged from 0 to 18, with low values (0–5) indicating elevated, well drained surfaces such as polygon rims or heaved ground surfaces, and high values (>12) identifying low-lying areas with potential for surface saturation and water accumulation. These areas often aligned with polygon troughs and low-centered basin polygon features, exhibiting the highest saturation potential, indicative of frequent or prolonged surface wetness. In contrast, lower TWI values were concentrated along elevated polygon rims or heaved ground features, reflecting well drained microtopography.

SPI was derived from specific catchment area and local slope to estimate the erosive potential of surface flow (Figure 11a). Raw SPI values ranged from 0 to 20,774.85, with peak values corresponding to localized zones of steep terrain combined with large upslope contributing areas. To capture relevant microtopographic variation, values were rescaled to a range of 0–80. Ranges from 0 to 20 displayed minimal erosive force, most notably observed near polygon interiors, ranges from 20 to 60 captured moderate erosion potential which were observed in the troughs and shallow channels, and ranges >60 were high erosion potential located in trough intersections and slope depressions. Areas with high SPI values were concentrated along polygon trough intersections and sloped basin margins, indicating zones of enhanced hydrologic energy and potential thermal erosion. Lower SPI values were associated with flat, saturated polygon centers and well drained rims, where overland flow lacks sufficient energy to mobilize sediment or deepen troughs.

Channel initiation threshold values represent the upslope contributing area, in raster cells, required for the onset of concentrated surface flow (Figure 11b). In this analysis, values ranged from 0 to 2630 cells, equivalent to 0–657.5 m². Lower thresholds (<500 cells) capture minor or ephemeral flow paths across flat terrain, while higher thresholds (>1000 cells) indicate persistent flow convergence zones, particularly along polygon trough intersections and degraded polygonal features. These values provide insight into where hydrologic energy is sufficient to initiate channel development and influence ice-wedge degradation.

The relationship between SPI and flow accumulation reveals a nonlinear increase in erosive potential with increasing contributing area across the polygonal tundra (Figure 12). SPI values remain low across much of the landscape despite increasing flow accumulation, reflecting the dominance of low slope gradients. However, once accumulation exceeds the threshold associated with trough convergence and basin outlets, SPI increases, indicating localized zones where hydrologic energy is sufficient to promote trough incision and thermal erosion. The upper envelope of the distribution captures high energy flow paths, which correspond closely with mapped polygon trough intersections and major drainage corridors, representing critical points of geomorphic sensitivity within the drainage network.

FIGURE 12

Drainage basin analysis was performed on a hydrologically conditioned DEM using the Wang & Liu sink filling algorithm to ensure continuous flow paths for accurate modeling (Wang and Liu, 2006). Flow direction and accumulation were computed using a multiple flow algorithm to delineate runoff pathways and accumulation zones.

The 62 basins exhibit low relief and limited topographic variation (Table 3), with a mean area of 2991 m², perimeter of 224.45 m, and average flow path of 76.46 m (maximum 783.5655 m) reflect basin geometries controlled by polygonal microtopography rather than channelized systems, while the low relief ratio (0.0042) is consistent with a flat Arctic coastal plain. Hydrologically, high drainage density (4.58 km/km²) and short mean channel lengths (13.7 m) indicate a finely dissected, highly branched trough network typical of ice wedge polygon terrain, with numerous short flow paths and ephemeral channels.

To quantify the relationship between delineation and modeled hydrologic connectivity, total catchment area outputs were thresholded to extract concentrated flow pathways, which were subsequently converted to vector features and compared to manually delineated polygon boundaries using nearest-distance analysis. Results indicate a strong spatial correspondence between modeled flow pathways and mapped polygon boundaries, with a mean distance of 1.81 m and a median distance of 1.03 m. Additionally, 75% of modeled flow-path points occurred within 2.32 m, of delineated boundaries, demonstrating that concentrated surface flow is strongly associated with polygon trough networks rather than polygon interiors. Visual overlay of TCA and manually derived polygon boundaries further supports this relationship by displaying that high contributing area values align with major trough intersections and low-lying drainage corridors (Figure 13).

FIGURE 13

The average overland flow length (109.2 m) suggests rapid routing of water into channels, though low gradients likely produce slow velocities and promote surface ponding. Although independent ground-truth data were unavailable, these results provide additional process-based agreement between polygon delineation and modeled hydrologic behavior. Spatial agreement and hydrologic coherence showed strong correspondence, including ∼99.5% intersection coverage and a moderate Jaccard Index (0.51). Because independent field validation data were unavailable, these results should not be interpreted as formal validation. Instead, strong spatial correspondence between mapped troughs, modeled flow accumulation pathways, channel networks, and elevated TWI/SPI zones suggests that the delineated polygons capture hydrologically meaningful terrain structures.

5 Discussion

The structural and hydrologic analyses presented in this study provide new insight into the representational limits of computationally derived polygon models and their implications for understanding permafrost surface dynamics. By comparing manually delineated ice-wedge polygons with Thiessen (Voronoi) polygons derived from intersection points, and integrating these spatial datasets with hydrologic terrain indices, this study elucidates how microtopography and flow connectivity govern the morphology and degradation of polygonal tundra surfaces. The significant differences in area, shape, and overlap underscore the limitations of idealized geometric models in representing the natural complexity of periglacial terrain. This distinction is particularly important because geometric oversimplification can lead to misinterpretation of geomorphic processes and landscape evolution, a concern noted in spatial partitioning studies such as Okabe et al. (2000), who emphasized that Voronoi-based representation, while computationally efficient, rarely reproduces true natural boundaries. While Thiessen polygons approximate the spatial arrangement of polygonal networks with moderate fidelity, they oversimplify natural variability that is critical to understanding thermokarst processes and hydrologic feedbacks in permafrost landscapes.

5.1 Structural representation and methodological implications

The geometric comparison revealed that Thiessen polygons produced larger mean areas and more regular geometric patterns relative to manually delineated polygon features, resulting in a nearly two-fold increase in mean polygon area and a higher shape compactness index. These differences reflect the inherent assumptions of the Thiessen method, which partitions space into equidistant zones around input points, thereby enforcing geometric uniformity that does not exist in natural polygonal terrain. In contrast, manually delineated polygons reflected greater irregularity, elongation, and asymmetry associated with visually interpreted terrain morphology and potential degradation features. The higher variability and wider range of Shape Index values in the manual dataset indicate that natural polygons evolve through differential subsidence and lateral thawing, processes that may not be fully represented by equidistant geometric constructs. Accordingly, the comparison between manually delineated and Thiessen polygons evaluates systematic geometric divergence between terrain expressed and idealized representations rather than classification accuracy.

Despite these simplifications, the high intersection coverage (≈99.5% of the manually delineated area) and moderate global Jaccard Index (0.51) suggest that Thiessen polygons maintain a general but incomplete spatial correspondence with observed polygonal structures. This indicates that automated polygon generation may serve as a reasonable first-order geometric approximation for mapping polygonal networks at large spatial scales or in data sparse regions. However, the reduced percent overlap and greater variability observed for Thiessen polygons underscore the limitations of such models in accurately representing degraded or hydrologically modified polygons. Rather than indicating that one method is inherently more accurate, these differences highlight how alternative delineation approaches emphasize different aspects of polygon geometry and may influence subsequent hydrologic interpretations. Thiessen polygons may underrepresent fine-scale trough asymmetry and localized surface irregularity, suggesting that manual or semi-automated approaches may be better suited when high-resolution geomorphic detail is required. This supports findings by Kanevskiy et al. (2017), who emphasized that accurate representation of polygon boundaries is crucial for assessing surface stability and predicting hydrological connectivity in permafrost dominated landscapes.

5.2 Microtopographic controls on polygon morphology

Terrain metrics derived from the lidar-based DEM revealed that even subtle elevation gradients (<1 m) exert strong controls over polygon structure and surface water routing. Slope gradients and relative elevation together govern the distribution of surface water, flow accumulation, and potential thaw zones. Low slopes combined with concave plan curvature delineate trough networks and depressions that promote flow convergence, corresponding to modeled flow paths and areas with high TWI values. Variations in slope and curvature define distinct geomorphic zones, with concave plan and profile curvature values corresponding to trough networks, and convex features marking elevated polygon centers and rims. These alternating convex-concave patterns reflect the inherent feedback between ice-wedge growth and surface drainage: troughs promote flow convergence and water accumulation, which in turn enhances thermal erosion and further trough deepening. Similar observations have been made in modeling studies of ice-wedge degradation under hydrologic influence (Liljedahl et al., 2016; Nitzbon et al., 2020; Zhang et al., 2003).

The Gaussian filtering approach applied in DEM preprocessing was particularly effective in improving the numerical stability of terrain derivatives in terms of relative spatial patterns rather than absolute gradient magnitudes, similar to the morphological filtering methods used by Zang et al. (2016) to enhance micro-relief detection in low-impact terrain systems. In addition, Gaussian smoothing of the DEM introduces a tradeoff between noise reduction and preservation of fine-scale microtopographic variability. While necessary for stable computation of slope and curvature metrics, smoothing may attenuate sharp gradients and reduce sensitivity of second order derivatives (e.g., plan and profile curvature) used to identify early-stage trough incision and subtle microtopographic transitions. Therefore, curvature-based interpretations should be considered conservative estimates of terrain roughness and flow convergence. RSP and closed depression mapping further corroborate this interpretation, identifying low-lying, saturated polygon centers as primary zones of thaw susceptibility. The inverse relationship between closed depression frequency and RSP further confirms that hydrologically isolated basins are preferentially concentrated in low-lying polygon centers and troughs, reinforcing the role of relative elevation in governing surface ponding. High LS-Factor values along trough margins and intersections indicate zones where flow convergence and extended lateral flow paths amplify surface runoff concentration, despite low local slope gradients. Under these conditions, LS-Factor variability is driven primarily by slope length and contributing area rather than by slope steepness alone. Collectively, these metrics demonstrate that polygon morphology and microtopographic relief, although subtle, are fundamental to the initiation and evolution of thaw induced hydrologic trough pathways. These findings align with broader permafrost studies detailing that fine-scale microtopography influences hydrologic and thermal feedbacks in polygonal tundra (Godin et al., 2016).

5.3 Hydrologic coupling and polygonal connectivity

Because the study area is underlain by continuous permafrost that effectively restricts vertical percolation, hydrologic patterns discussed here are interpreted in terms of surface and near surface flow organization rather than subsurface hydrologic fluxes (Bui et al., 2020; Wang and Liu, 2006). Under these conditions, surface microtopography exerts first-order control on lateral water redistribution across polygonal tundra landscapes. Due to flow routing being performed using a multiple flow direction (MFD) algorithm, hydrologic metrics emphasize relative patterns of diffuse surface connectivity across low-slope polygon interiors and interconnected trough networks, while potentially underrepresenting highly localized concentrated flow at the steepest trough margins.

The coexistence of ponded flow regimes and high runoff potential reflects thresh old controlled hydrologic behavior governed by polygonal microtopography. Under low energy conditions, minimal regional gradients and shallow micro-depressions promote surface ponding within polygon interiors. As precipitation or thaw inputs increase, vertical drainage remains restricted by the shallow active layer and underlying permafrost, such that local storage capacity is exceeded and lateral runoff is rapidly routed into polygon trough networks. The strong coupling between horizontal and vertical overland flow distance further indicates that low-RSP surfaces facilitate short, efficient flow paths with minimal elevation loss, whereas elevated polygon interiors require longer planimetric travel distances to achieve comparable downslope connectivity. These networks concentrate flow and facilitate efficient lateral drainage despite low overall relief, enabling rapid transitions between ponded and runoff dominated states.

Hydrologic indices revealed a strong coupling between polygonal troughs and surface flow networks. Flow accumulation and channel distance analyses indicated that major troughs correspond closely with modeled drainage paths, confirming their role as primary conduits for overland flow and thermal erosion. This spatial correspondence is also qualitatively supported by the total catchment area and specific catchment area maps, where modeled flow pathways visibly align with trough networks and ponded surface features (Figure 10). While this does not constitute independent field validation, it provides additional process-based agreement between modeled hydrologic pathways and observable surface expressions of water redistribution. Joint density analysis of flow accumulation and RSP further reveals that high contributing areas are confined to low-RSP terrain, while elevated polygon rims remain connected even where flow paths exist. This pattern highlights RSP as an organizer of hydrologic connectivity within polygonal tundra, constraining where surface runoff can effectively concentrate. This close correspondence between mapped trough networks and modeled surface hydrology provides functional validation of polygon delineations, consistent with prior studies where field access is limited and process-based agreement is used as a validation proxy. The drainage density (4.58 km/km2) and short mean flow paths (∼13.7 m) indicate a finely dissected landscape where water movement is highly localized yet spatially continuous through interconnected troughs. These findings support previous work by McCarthy (1994), where results identified surface drainage as a key control on near surface permafrost hydrology and thaw behavior along the Arctic coastal plain. Subsequent studies have similarly shown that trough networks enhance lateral drainage and accelerate thermal degradation by concentrating surface water and increasing ground heat flux (Liljedahl et al., 2016; Nitzbon et al., 2020). These characteristics are consistent with Arctic polygonal lowlands, where shallow gradients and low relief promote extensive surface connectivity despite limited vertical relief (Abolt et al., 2018; Wales et al., 2020).

The distribution of hydrologic indices further clarifies the spatial coupling between microtopography and hydrologic processes. High TWI values along troughs and basin intersections correspond to zones of persistent saturation, while elevated SPI values mark localized areas of enhanced erosion potential at trough junctions and basin outlets where flow convergence and accumulated discharge increase flow energy following a nonlinear scaling relationship between contributing area and erosive potential, despite generally low local slope gradients. These results support the conceptual models of Beven and Kirkby (1979) and Moore et al. (1991), which link TWI and SPI to flow accumulation and local slope as proxies for runoff generation and erosion potential (Beven and Kirkby, 1979; McCarthy, 1994). These spatial patterns support the hypothesis that ice-wedge degradation and thermokarst expansion are intensified where hydrologic convergence coincides with geomorphic concavity. Conversely, polygon interiors with low TWI and SPI values exhibit limited connectivity and minimal flow energy, favoring persistence of ponded conditions or stable rims. The low channel initiation thresholds identified in these regions suggest that small perturbations in surface hydrology, such as increased precipitation or active layer deepening, could trigger new flow paths and accelerate thaw-driven redistributed polygon networks. These inferences align with recent modeling, which emphasizes lateral drainage and hydrologic feedbacks in ice-wedge polygon systems (Abolt et al., 2018; Wales et al., 2020; Zlotnik et al., 2020).

5.4 Implications for permafrost degradation and modeling

The combined geometric and hydrologic results underscore that hydrologically mediated feedbacks are central to the evolution of ice-wedge polygon networks. The correspondence between modeled flow convergence and polygon trough morphology indicates that permafrost degradation progresses preferentially along established hydrologic trough pathways, reinforcing trough connectivity and promoting lateral thaw. At the same time, the tendency of Thiessen polygons to overestimate spatial extent and homogenize shape implies that computationally derived networks may underestimate the true heterogeneity of hydrologic influence within permafrost terrain. This has implications for scaling local observations to landscape or regional models, especially in Arctic coastal lowlands where microtopographic and hydrologic processes dominate.

For modeling and monitoring applications, these findings suggest that while Thiessen-based or other automated polygon delineation methods can approximate overall polygon structure, they should be coupled with terrain and flow metrics to capture hydrologically significant variability. Incorporating high-resolution DEM derivatives such as TWI, SPI, curvature, and flow accumulation into polygon delineation and hydrologic analysis can substantially improve the true nature of the polygon network representations and support more accurate predictions of thaw susceptibility. As Arctic regions experience heightened permafrost degradation, such integrated approaches will be increasingly important for bridging geomorphic process understanding with remote sensing and modeling frameworks (Liljedahl et al., 2016; Nitzbon et al., 2020).

The hydrologic connectivity patterns identified in this study are consistent with recent evidence that preferential flow pathways can accelerate permafrost degradation across the landscape. The concentration of high flow accumulation, elevated wetness potential, and erosive energy along polygon trough networks supports findings that lateral drainage enhances thaw-driven subsidence and geomorphic instability (Abolt et al., 2020; Harp et al., 2021). Recent studies have also displayed that degrading ice wedge polygons can increase lateral water and dissolved organic carbon export (Speetjens et al., 2024), contribute to broader greenhouse gas emissions (Freitas et al., 2025), and amplify hydrologically driven mass wasting across Arctic landscapes (Young et al., 2026). Collectively, these findings further validate that fine-scale hydrologic connectivity may serve as an important indicator of future permafrost instability and carbon cycle feedbacks.

5.5 Limitations and future work

While this study provides a detailed assessment of polygonal morphology and hydrologic coupling, several limitations warrant consideration. In addition, the analysis is based on a single-time lidar acquisition, which captures polygon morphology and hydrologic organization at a specific seasonal state. While this is appropriate for evaluating spatial relationships between microtopography and surface hydrology, it does not resolve seasonal or interannual variability in thaw depth, surface subsidence, or hydrologic connectivity. As a result, the findings should be interpreted as indicators of hydrologic potential and spatial susceptibility to future degradation, rather than direct observations of temporal thaw evolution or degradation dynamics. Future studies integrating multi-temporal lidar, optical, or synthetic aperture radar (SAR) observations across seasons and years would enable direct assessment of dynamic polygon evolution and hydrologic change. Accordingly, interpretations emphasize relative geometric organization and systematic divergence between idealized and terrain constrained polygonal structures rather than absolute classification accuracy or absolute hydrologic magnitudes. Because field-based boundary validation was not available, absolute positional accuracy of individual polygon edges cannot be quantified. Accuracy is interpreted in terms of relative spatial consistency and process-based agreement rather than absolute boundary correctness. Field-based validation of polygon boundaries and hydrologic flow paths was unavailable for this study. While high-resolution lidar provides detailed representation of surface microtopography, the absence of in situ measurements limits direct assessment of absolute boundary accuracy and hydrologic fluxes; however, the results presented highlight the importance of understanding hydrologic pathways under thaw-driven geomorphic processes.

Although high-resolution lidar and hydrologic modeling captured fine-scale microtopography, the static nature of the DEM restricts the temporal interpretation of thaw dynamics. Future analyses integrating multi-temporal or seasonal datasets could better quantify active surface deformation and evolving hydrologic connectivity. Second, model parameters such as flow accumulation thresholds and smoothing algorithms influence the delineation of channel networks and basin boundaries. In addition, the use of a multiple flow direction routing algorithm represents a scale dependent modeling choice that emphasizes diffuse surface connectivity across low-slope polygon interiors. While appropriate for capturing dominant surface flow organization under continuous permafrost, this approach may underestimate localized concentrated flow along steeper trough margins. Similarly, although DEM smoothing improves numerical stability of terrain derivatives, some attenuation of sharp microtopographic gradients is unavoidable, and hydrologic metrics are therefore interpreted in terms of relative spatial patterns rather than absolute gradient or erosion magnitudes. Finally, field validation and coupling with thermal and soil moisture data will enhance confidence in model outputs and clarify the physical mechanisms linking hydrology and perma-frost degradation (Abolt et al., 2018; Wales et al., 2020). Continued integration of remote sensing, process-based modeling, and in situ observations will be essential to refine our understanding of how hydrologic feedbacks drive landscape change in Arctic polygonal tundra.

6 Conclusion

This study developed and applied a remote sensing and hydrologic modeling framework to analyze the spatial structure and connectivity of ice-wedge polygon networks within Arctic permafrost terrain. Using 0.5 m lidar-derived DEM, detailed terrain derivatives, and hydrologic indices, this study characterized how microtopography and surface water processes interact to shape polygon morphology and degradation. The findings demonstrate that even subtle topographic gradients exert strong control over water routing, thaw dynamics, and the ongoing evolution of polygonal landscapes.

The geometric comparison between manually delineated and Thiessen-derived polygons revealed that, although Thiessen polygons broadly capture the spatial footprint of polygon networks, they fail to reproduce the irregular and heterogeneous morphology characteristic of natural periglacial terrain. Manual delineations displayed greater variability in shape and size, reflecting natural deformation from thaw subsidence, trough widening, and localized hydrologic feedbacks. These results underscore the limitations of idealized geometric models such as Thiessen polygons in representing terrain where microtopography and hydrologic processes are tightly coupled, a limitation also highlighted by Okabe et al. (2000). By oversimplifying irregular boundaries, geometric models may obscure key relationships between polygon structure, water flow, and surface stability.

Hydrologic modeling effectively captured the distribution and intensity of surface flow processes that sustain and modify the polygon trough network. Modeled parameters, including horizontal and vertical overland flow distance, flow accumulation, and channel network distance, closely aligned with mapped troughs, confirming that the trough system acts as the primary hydrologic conduit across the polygon terrain. These features represent active flow corridors for surface runoff and seasonal meltwater redistribution, consistent with prior observations of hydrologic-geomorphic coupling in polygonal tundra (McCarthy, 1994; Kanevskiy et al., 2017). Zones characterized by short overland flow distances and high flow accumulation corresponded to depressions and trough intersections, emphasizing that hydrologic convergence reinforces trough incision and promotes the persistence of polygonal boundaries.

Terrain-based hydrologic indices further illuminated the coupling between water accumulation, erosion, and topographic form. The TWI and SPI both revealed spatial relationships between surface saturation and erosive potential. High TWI values corresponded to low-lying troughs and depressions indicative of frequent saturation, while elevated SPI values occurred along intersecting troughs and channelized zones, suggesting enhanced erosion and thaw. These relationships support the conceptual framework proposed by Beven and Kirkby (1979) and later refined by Moore et al. (1991), which links terrain-driven hydrologic potential to surface process intensity (Beven and Kirkby, 1979; Moore et al., 1991). Similarly, Gaussian smoothing and curvature analyses confirmed that even low-relief microtopography governs localized water convergence and flow pathways, consistent with the findings of Zhang et al. (2003). RSP emerged as a key organizing variable linking microtopography and hydrologic function, with low-RSP terrain consistently associated with closed surface depressions, high flow accumulation, and enhanced hydrologic connectivity. This confirms that relative elevation within the polygon terrain, rather than absolute slope magnitude, controls where surface water ponds, spills, and concentrates across polygonal tundra.

Drainage basin and flow network analyses demonstrated that the hydrologic connectivity of polygonal tundra is governed by both local trough geometry and broader flow organization. The high drainage density and short overland flow lengths identified in this study indicate a finely dissected microtopographic network where numerous channels facilitate rapid yet spatially constrained runoff. This structure mirrors the underlying ice-wedge system and underscores that hydrologic processes are not passive outcomes but active drivers of polygonal formation, maintenance, and degradation. Joint density and SPI analyses further indicate that hydrologic influence on polygon degradation is threshold controlled, with nonlinear increases in erosive potential occurring once contributing area exceeds levels associated with trough convergence and basin outlets. These threshold zones represent critical points of geomorphic sensitivity where modest increases in runoff or thaw input may disproportionately accelerate trough incision and ice-wedge degradation.

Overall, the presented study demonstrates that hydrologic processes—captured through high-resolution DEM analysis, terrain modeling, and spatial hydrologic indices—accurately represent the trough networks that govern polygon formation and transformation. By integrating geometric and hydrologic approaches, this study provides a robust methodological framework for quantifying the feedbacks between surface water dynamics and permafrost morphology, offering insights critical to predicting permafrost response to hydrologic processes.

References

• 1

  AboltC. J.YoungM. H. (2020). High-resolution mapping of spatial heterogeneity in ice wedge polygon geomorphology near prudhoe Bay, Alaska. Sci. Data7 (1), 87. 10.1038/s41597-020-0423-9

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  AboltC. J.YoungM. H.AtchleyA. L.HarpD. R. (2018). Microtopographic control on the ground thermal regime in ice wedge polygons. Cryosphere12 (6), 1957–1968. 10.5194/tc-12-1957-2018

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  AboltC. J.YoungM. H.AtchleyA. L.HarpD. R.CoonE. T. (2020). Feedbacks between surface deformation and permafrost degradation in ice wedge polygons, arctic coastal plain, Alaska. J. Geophys. Res. Earth Surf.125 (3), e2019JF005349. 10.1029/2019JF005349

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  AurenhammerF. (1991). Voronoi diagrams—a survey of a fundamental geometric data structure. ACM Computing Surveys (CSUR)23 (3), 345–405. 10.1145/116873.116880

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  BaileyG.LiY.McKinneyN.YoderD.WrightW.HerreroH. (2022). Comparison of ground point filtering algorithms for high-density point clouds collected by terrestrial LiDAR. Remote Sens.14 (19), 4776. 10.3390/rs14194776

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  BevenK. J.KirkbyM. J. (1979). A physically based, variable contributing area model of basin hydrology. Hydrological Sciences Journal24 (1), 43–69. 10.1080/02626667909491834

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  BiskabornB. K.SmithS. L.NoetzliJ.MatthesH.VieiraG.StreletskiyD. A.et al (2019). Permafrost is warming at a global scale. Nat. Communications10 (1), 264. 10.1038/s41467-018-08240-4

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  BlackR. F. (1976). Periglacial features indicative of permafrost: ice and soil wedges. Quat. Res.6 (1), 3–26. 10.1016/0033-5894(76)90037-5

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  BöhnerJ.KotheR.ConradO.GrossJ.RingelerA.SeligeT. (2002). in Soil Regionalization by Means of Terrain Analysis and Process Parameterization (European Soil Bureau from the European Union), 7.

  • Google Scholar
  • View reference in article

• 10

  BoikeJ.KattenstrothB.AbramovaK.BornemannN.ChetverovaA.FedorovaI.et al (2013). Baseline characteristics of climate, permafrost and land cover from a new permafrost observatory in the Lena river Delta, Siberia (1998–2011). Biogeosciences10 (3), 2105–2128. 10.5194/bg-10-2105-2013

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  BrownJ.FerriansO.HeginbottomJ. A.MelnikovE. S. (1997). Circum-arctic map of permafrost and ground-ice conditions: U.S. geological survey circum-pacific map 45. 10.3133/cp45

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  BuiM. T.LuJ.NieL. (2020). A review of hydrological models applied in the permafrost-dominated arctic region. Geosciences10 (10), 401. 10.3390/geosciences10100401

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  BurroughP. A.McDonnellR. A.LloydC. D. (2015). Principles of Geographical Information Systems. Oxford University Press.

  • Google Scholar
  • View reference in article

• 14

  CaiS.YuS.HuiZ.TangZ. (2023). ICSF: an improved cloth simulation filtering algorithm for airborne LiDAR data based on morphological operations. Forests14 (8), 1520. 10.3390/f14081520

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  DesmetP. J.GoversG. (1996). A GIS procedure for automatically calculating the USLE LS factor on topographically complex landscape units. J. Soil Water Conservation51 (5), 427–433. 10.1080/00224561.1996.12457102

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  FreemanT. G. (1991). Calculating catchment area with divergent flow based on a regular grid. Comput. & Geosciences17 (3), 413–422. 10.1016/0098-3004(91)90048-I

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  FreitasN. L.Walter AnthonyK.LenzJ.PorrasR. C.TornM. S. (2025). Substantial and overlooked greenhouse gas emissions from deep arctic lake sediment. Nat. Geosci.18 (1), 65–71. 10.1038/s41561-024-01614-y

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  FrenchH. M. (2018). The Periglacial Environment. John Wiley & Sons.

  • Google Scholar
  • View reference in article

• 19

  GlennieC. L.CarterW. E.ShresthaR. L.DietrichW. E. (2013). Geodetic imaging with airborne LiDAR: the earth's surface revealed. Rep. Prog. Phys.76 (8), 086801. 10.1088/0034-4885/76/8/086801

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  GodinE.FortierD.LévesqueE. (2016). Nonlinear thermal and moisture response of ice-wedge polygons to permafrost disturbance increases heterogeneity of high Arctic wetland. Biogeosciences13 (5), 1439–1452. 10.5194/bg-13-1439-2016

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  GruberS.PeckhamS. (2009). Land-surface parameters and objects in hydrology. Dev. Soil Science33, 171–194. 10.1016/S0166-2481(08)00007-X

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  HarpD. R.ZlotnikV.AboltC. J.BuseyB.AvendañoS. T.NewmanB. D.et al (2021). New insights into the drainage of inundated ice-wedge polygons using fundamental hydrologic principles. Cryosphere15 (8), 4005–4029. 10.5194/tc-15-4005-2021

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  HornB. K. (2005). Hill shading and the reflectance map. Proc. IEEE69 (1), 14–47. 10.1109/PROC.1981.11918

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  HuangH.ChenX.WangX.WangX.LiuL. (2019). A depression-based index to represent topographic control in urban pluvial flooding. Water11 (10), 2115. 10.3390/w11102115

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  JensonS. K.DomingueJ. O. (1988). Extracting topographic structure from digital elevation data for geographic information system analysis. Photogrammetric Engineering Remote Sensing54 (11), 1593–1600.

  • Google Scholar
  • View reference in article

• 26

  JiangA. L.HsuK.SandersB. F.SorooshianS. (2023). Topographic hydro-conditioning to resolve surface depression storage and ponding in a fully distributed hydrologic model. Adv. Water Resources176, 104449. 10.1016/j.advwatres.2023.104449

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  JonesB. M.StokerJ. M.GibbsA. E.GrosseG.RomanovskyV. E.DouglasT. A.et al (2013). Quantifying landscape change in an arctic coastal lowland using repeat airborne LiDAR. Environ. Res. Lett.8 (4), 045025. 10.1088/1748-9326/8/4/045025

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  JorgensonM. T.ShurY. L.PullmanE. R. (2006). Abrupt increase in permafrost degradation in Arctic Alaska. Geophys. Res. Lett.33 (2), , 1–4. 10.1029/2005GL024960

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  KanevskiyM.ShurY.JorgensonT.BrownD. R.MoskalenkoN.BrownJ.et al (2017). Degradation and stabilization of ice wedges: implications for assessing risk of thermokarst in northern Alaska. Geomorphology297, 20–42. 10.1016/j.geomorph.2017.09.001

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  LaraM. J.McGuireA. D.EuskirchenE. S.TweedieC. E.HinkelK. M.SkurikhinA. N.et al (2015). Polygonal tundra geomorphological change in response to warming alters future CO 2 and CH 4 flux on the barrow peninsula. Glob. Change Biol.21 (4), 1634–1651. 10.1111/gcb.12757

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  LeffingwellE. D. K. (1915). Ground-ice wedges: the dominant form of ground-ice on the north coast of Alaska. J. Geol.23 (7), 635–654. 10.1086/622281

  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  LiljedahlA. K.HinzmanL. D.SchullaJ. (2012). “Ice-wedge polygon type controls low-gradient watershed-scale hydrology,” in Proceedings of the Tenth International Conference on Permafrost (Salekhard, Russia: The Northern Publisher), 1, 231–236.

  • Google Scholar
  • View reference in article

• 33

  LiljedahlA. K.BoikeJ.DaanenR. P.FedorovA. N.FrostG. V.GrosseG.et al (2016). Pan-arctic ice-wedge degradation in warming permafrost and its influence on tundra hydrology. Nat. Geosci.9 (4), 312–318. 10.1038/ngeo2674

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  LiljedahlA. K.WitharanaC.ManosE. (2024). The capillaries of the Arctic tundra. Nat. Water2 (7), 611–614. 10.1038/s44221-024-00276-9

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  MackayJ. R. (1990). Some observations on the growth and deformation of epigenetic, syngenetic and anti‐syngenetic ice wedges. Permafr. Periglac. Process.1 (1), 15–29. 10.1002/ppp.3430010104

  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  MackayJ. R. (2000). Thermally induced movements in ice-wedge polygons, western Arctic coast: a long-term study. Géogr. Physique Quaternaire54 (1), 41–68. 10.7202/004846ar

  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  MartzL. W.GarbrechtJ. (1998). The treatment of flat areas and depressions in automated drainage analysis of raster digital elevation models. Hydrol. Processes12 (6), 843–855. 10.1002/(SICI)1099-1085(199805)12:6%3C843::AID-HYP658%3E3.0.CO;2-R

  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  McCarthyK. A. (1994). Overview of Environmental and Hydrogeologic Conditions at Barrow, Alaska: U.S. Geological Survey Open-File Report. 10.3133/ofr94322

  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  MontgomeryD. R.DietrichW. E. (1989). Source areas, drainage density, and channel initiation. Water Resources Research25 (8), 1907–1918. 10.1029/WR025i008p01907

  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  MooreI. D.BurchG. J. (1986). Physical basis of the length slope factor in the universal soil loss equation. Soil Sci. Soc. Am. J.50 (5), 1294–1298. 10.2136/sssaj1986.03615995005000050042x

  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  MooreI. D.GraysonR. B.LadsonA. R. (1991). Digital terrain modelling: a review of hydrological, geomorphological, and biological applications. Hydrol. Processes5 (1), 3–30. 10.1002/hyp.3360050103

  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  MoritzR. E. (1979). Institute of Arctic and Alpine Research University of Colorado (Boulder). Synoptic Climatology of the Beaufort Sea Coast of Alaska.

  • Google Scholar
  • View reference in article

• 43

  NitzbonJ.LangerM.WestermannS.MartinL.AasK. S.BoikeJ. (2019). Pathways of ice-wedge degradation in polygonal tundra under different hydrological conditions. Cryosphere13 (4), 1089–1123. 10.5194/tc-13-1089-2019

  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  NitzbonJ.WestermannS.LangerM.MartinL. C.StraussJ.LaboorS.et al (2020). Fast response of cold ice-rich permafrost in northeast Siberia to a warming climate. Nat. Communications11 (1), 2201. 10.1038/s41467-020-15725-8

  • CrossRef
  • Google Scholar
  • View reference in article

• 45

  O'CallaghanJ. F.MarkD. M. (1984). The extraction of drainage networks from digital elevation data. Comput. Vision, Graphics, Image Processing28 (3), 323–344. 10.1016/s0734-189x(84)80011-0

  • CrossRef
  • Google Scholar
  • View reference in article

• 46

  OkabeA.BootsB.SugiharaK.ChiuS. N. (2000). in Spatial Tessellations: Concepts and Applications of Voronoi Diagrams (John Wiley & Sons), 501.

  • Google Scholar
  • View reference in article

• 47

  PeckhamS. D. (2009). Geomorphometry and spatial hydrologic modelling. Dev. Soil Sci.33, 579–602. 10.1016/S0166-2481(08)00025-1

  • CrossRef
  • Google Scholar
  • View reference in article

• 48

  QuinnP. F. B. J.BevenK.ChevallierP.PlanchonO. (1991). The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models. Hydrol. Processes5 (1), 59–79. 10.1002/hyp.3360050106

  • CrossRef
  • Google Scholar
  • View reference in article

• 49

  RIEGL Laser Measurement Systems GmbH (2021). VUX-120-23 Airborne Laser Scanner: Technical Description and System Overview. Austria: RIEGL: Horn.

  • Google Scholar
  • View reference in article

• 50

  RomanovskyV. E.SmithS. L.ChristiansenH. H. (2010). Permafrost thermal state in the polar Northern Hemisphere during the international polar year 2007–2009: a synthesis. Permafr. Periglac. Processes21 (2), 106–116. 10.1002/ppp.689

  • CrossRef
  • Google Scholar
  • View reference in article

• 51

  SabirovaM.FedorenkoR.AfanasyevI. (2019). “Ground profile recovery from aerial 3D LiDAR-Based maps,” in 2019 24th Conference of Open Innovations Association (FRUCT) (Moscow, Russia), 367–374. 10.23919/FRUCT.2019.8711928

  • CrossRef
  • Google Scholar
  • View reference in article

• 52

  SarıtaşB.KaplanG. (2023). Enhancing ground point extraction in airborne LiDAR point cloud data using the CSF filter algorithm. Adv. LiDAR3 (2), 53–61.

  • Google Scholar
  • View reference in article

• 53

  SeibertJ.McGlynnB. L. (2007). A new triangular multiple flow direction algorithm for computing upslope areas from gridded digital elevation models. Water Resources Research43 (4), 1–8. 10.1029/2006WR005128

  • CrossRef
  • Google Scholar
  • View reference in article

• 54

  SharyP. A.SharayaL. S.MitusovA. V. (2002). Fundamental quantitative methods of land surface analysis. Geoderma107 (1-2), 1–32. 10.1016/S0016-7061(01)00136-7

  • CrossRef
  • Google Scholar
  • View reference in article

• 55

  ShurY.JonesB. M.JorgensonM. T.KanevskiyM. Z.LiljedahlA.WalkerD. A.et al (2025). Formation of low-centered ice-wedge polygons and their orthogonal systems: a review. Geosciences15 (7), 249. 10.3390/geosciences15070249

  • CrossRef
  • Google Scholar
  • View reference in article

• 56

  SpeetjensN. J.BerghuijsW. R.WagnerJ.VonkJ. E. (2024). Degradation of ice-wedge polygons leads to increased fluxes of water and DOC. Sci. Total Environ.920, 170931. 10.1016/j.scitotenv.2024.170931

  • CrossRef
  • Google Scholar
  • View reference in article

• 57

  TagilS.JennessJ. (2008). GIS-based automated landform classification and topographic, landcover and geologic attributes of landforms around the Yazoren Polje, Turkey. J. Appl. Sci.8, 910–921. 10.3923/jas.2008.910.921

  • CrossRef
  • Google Scholar
  • View reference in article

• 58

  TarbotonD. G. (1997). A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resources Research33 (2), 309–319. 10.1029/96WR03137

  • CrossRef
  • Google Scholar
  • View reference in article

• 59

  ThiessenA. H. (1911). Precipitation averages for large areas. Mon. Weather Review39 (7), 1082–1089. 10.1175/1520-0493(1911)39%3C1082b:PAFLA%3E2.0

  • CrossRef
  • Google Scholar
  • View reference in article

• 60

  TothC.JóźkówG. (2016). Remote sensing platforms and sensors: a survey. ISPRS J. Photogrammetry Remote Sens.115, 22–36. 10.1016/j.isprsjprs.2015.10.004

  • CrossRef
  • Google Scholar
  • View reference in article

• 61

  WalesN. A.Gomez-VelezJ. D.NewmanB. D.WilsonC. J.DafflonB.KneafseyT. J.et al (2020). Understanding the relative importance of vertical and horizontal flow in ice-wedge polygons. Hydrology Earth Syst. Sci.24 (3), 1109–1129. 10.5194/hess-24-1109-2020

  • CrossRef
  • Google Scholar
  • View reference in article

• 62

  WangL.LiuH. (2006). An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling. Int. J. Geogr. Inf. Sci.20 (2), 193–213. 10.1080/13658810500433453

  • CrossRef
  • Google Scholar
  • View reference in article

• 63

  WilsonJ. P.GallantJ. C. (2000). Terrain Analysis: Principles and Applications. Hoboken, NJ, USA: Wiley.

  • Google Scholar
  • View reference in article

• 64

  WischmeierW. H.SmithD. D. (1978). Predicting Rainfall Erosion Losses: A Guide to Conservation Planning (No. 537). Washington, DC: U.S. Department of Agriculture.

  • Google Scholar
  • View reference in article

• 65

  WooM. K. (2012). Permafrost Hydrology. Springer Science & Business Media.

  • Google Scholar
  • View reference in article

• 66

  YilmazV. (2021). Automated ground filtering of LiDAR and UAS point clouds with metaheuristics. Opt. & Laser Technol.138, 106890. 10.1016/j.optlastec.2020.106890

  • CrossRef
  • Google Scholar
  • View reference in article

• 67

  YoungJ. M.FarquharsonL.LuoJ.NesterovaN.Van der SluijsJ.KokeljS. V. (2026). Permafrost mass wasting in ice‐rich landscapes: recent advances (2013 to 2024) on mechanisms, dynamics and impacts. Permafr. Periglac. Process.37 (2), 284–304. 10.1002/ppp.70015

  • CrossRef
  • Google Scholar
  • View reference in article

• 68

  ZevenbergenL. W.ThorneC. R. (1987). Quantitative analysis of land surface topography. Earth Surface Processes Landforms12 (1), 47–56. 10.1002/esp.3290120107

  • CrossRef
  • Google Scholar
  • View reference in article

• 69

  ZhangK.ChenS. C.WhitmanD.ShyuM. L.YanJ.ZhangC. (2003). A progressive morphological filter for removing nonground measurements from airborne LIDAR data. IEEE Transactions Geoscience Remote Sensing41 (4), 872–882. 10.1109/TGRS.2003.810682

  • CrossRef
  • Google Scholar
  • View reference in article

• 70

  ZhangW.QiJ.WanP.WangH.XieD.WangX.et al (2016). An easy-to-use airborne LiDAR data filtering method based on cloth simulation. Remote Sensing8 (6), 501. 10.3390/rs8060501

  • CrossRef
  • Google Scholar
  • View reference in article

• 71

  ZlotnikV. A.HarpD. R.JafarovE. E.AboltC. J. (2020). A model of ice wedge polygon drainage in changing Arctic terrain. Water12 (12), 3376. 10.3390/w12123376

  • CrossRef
  • Google Scholar
  • View reference in article