Morphometric Analysis of Groundwater Icings: Intercomparison of Estimation Techniques

: Groundwater icings, typical features of permafrost hydrology, are indicative of hydrothermal interactions between surface and ground waters, and permafrost. Their main morphological parameters, i.e., icing area and volume, are generally estimated with low accuracy. Only scarce field observational data on icing volume and seasonal development exist to date. Our study evaluates and compares performance of several widely used techniques of icing morphometric estimation, based on field data, collected on a giant Icing #2 in the Samokit River basin, southern Yakutia. Groundwater icing area was estimated by: (a) staking, (b) unmanned aerial vehicle (UAV) surveys, and (c) satellite imagery analysis. Icing #2 area in late February was between 1.38·10 6 m 2 and 1.68·10 6 m 2 , icing volume, between 1.73·10 6 m 3 and 4.20·10 6 m 3 , depending on the technique used. Staking is the least accurate, but also the only direct technique, which is hence used as a baseline tool in our study. Staking-based assessment of icing morphometry is the most conservative, while UAV-based estimates of icing area are higher by 14% to 17%, and of icing volume, by 74% to 142%, compared to staking. The latter appears, in our case, to be the least accurate method, although a direct one. It requires a sufficient number of staking points and transects, which should be set up to represent all icing zones, i.e., channel branches and alluvial islands. Photogrammetry based on UAV surveys has numerous advantages, i.e., higher precision of a per pixel icing volume calculation, based on an ice-free valley bottom digital surface model (DSM), and potential reusability of a resulting DSM. However, positioning precision suffers from the overlay of multiple flyovers required because of battery replacements, and, in our case, an insufficient number of ground control points. Satellite imagery along with B.L. Sokolov’s empirical approach were used to estimate the annual maximum icing area and volume, and the empirical estimates tend to converge to satellite-based values. Finally, all thing being equal, UAV-based photogrammetry shows higher precision in estimating the icing morphometrical parameters.

Numerous direct and indirect methods have been developed for observing or estimating icing area and volume. Historically, staking is by far the most widespread technique for icing volume monitoring (see detailed descriptions of this and other methods in the Methods section). Geophysical surveys aimed at reconstructing the ice bottom topography allow estimation of ice thickness and volume [30][31][32]. Recently, photogrammetry based on unmanned aerial vehicle (UAV) surveys has emerged as a novel data source for digital surface models (DSM), allowing fast, inexpensive, and relatively accurate estimation of surface deformations, e.g., landslides, icing, thermokarst ground subsidence, and other processes [28,33,34].
Remote sensing of permafrost environment allows more detailed indirect studies of icing morphology and dynamics [5,24,35,36]. Satellite imagery interpretation allowed Yoshikawa et al. [5] to track icing distribution, groundwater temperature and discharge, water residence time, and the influence of icing on winter river runoff in the Brooks Range region, northern Alaska, USA. The interferometry synthetic aperture radar (InSAR) is a perspective technique for estimation of landform morphology changes [37,38]. However, open-source SAR mission products, i.e., Sentinel-1 data acquired in C-band, have vertical resolution more than 2 m [39], which is comparable to an average icing thickness; besides, the accuracy of InSAR products can degrade owing to temporal decorrelation effect [40].
Numerous empirical methods have been developed for converting icing areas to volumes. In Russia, the approaches developed by B.L. Sokolov [13] are widely recognized. Sokolov calculated morphometric parameters of icings in various environments and climates of the USSR, and introduced several empirical techniques for evaluating maximum icing area and volume [13].
This paper aims to compare several techniques for icing area and volume estimation, based on both published and observed data from the Samokit River valley in southern Yakutia, Siberian Russia.

Study Site
Field studies were conducted in the Samokit River valley, about 40 km west of Neryungri city, Southern Yakutia ( Figure 1). The Samokit River drains the north-eastern flank of the Zverev's Ridge and crosses, in its lower section, the Chulman Plateau inheriting the linear shape of its valley from the local tectonic faultline, associated with high geothermal fluxes [27]. The study area lies in a zone of discontinuous permafrost, which occurs in valleys, on slopes, and on mountain tops above 1400 m above sea level (a.s.l.). Maximum permafrost thickness exceeds 400 m on mountain tops and gradually decreases to 150-200 m in river valleys [41]. The permafrost temperature varies from -1 to 0°C across the region, and reaches -4 °C in the mountainous areas. Taliks occur on flat watersheds between 800 and 1200 m a.s.l., and south-oriented slopes [41]. Open taliks are usually hydraulically connected to water-bearing tectonic fissures in valley bottoms, including the Samokit River valley [42], serving as pathways of groundwater discharge at the valley bottom. Sub-permafrost groundwater recharge, transit, and discharge are mainly controlled by regional geologic structures, tectonics, geomorphology, and geocryology. A detailed assessment of subpermafrost groundwater storage of the Samokit aquifer was performed by the South Yakutian Geological Enterprise (SYGE) between 1978 and 1985 [43]. Thirty exploration and production water wells were drilled in the valley. Hydrological, hydrogeological, and geothermic studies were conducted, and the geological structure was studied by geoelectrical survey [43].
Icing-sourcing open sub-channel talik zones and sub-permafrost aquifers belong to local groundwater flow, according to Toth [44], corresponding to free water exchange regime, after Fotiev [42], in relation to the Samokit River valley. Open sub-channel talik aquifers exist within the alluvial fill of the valley, about 10 m thick, and are fed by supra-permafrost groundwater and, predominantly, by riverine water. Its low total dissolved solids (TDS) content, about 0.07 g/L, matches that of the Samokit River runoff. The aquifer-free table corresponds to the river stage in summer and rises during winter under the piston effect of seasonal freezing.
Sub-permafrost groundwater is found at depths below 10 m within the valley, below 70 m on the slopes, and below 200 m in the watershed, as revealed by many boreholes drilled in the 1980s, and acts as confined aquifers with hydraulic head over 8 m at the valley bottom. Sub-permafrost groundwater has low TDS content, varying from 0.05 to 0.15 g/L, its table responds quickly to precipitation events. Groundwater temperature Tw in boreholes varies from 0.5 °C to 4 °C, depending on the season.
The Samokit River is a typical mountainous river, featuring a steep valley slope of around 8 m per kilometer length and high flow velocity. The river length is 33 km and the basin area is 233 km 2 [43]. Its discharge under winter low-flow conditions (February-March) is about 0.4 to 0.5 m 3 /s and is sourced by groundwater. The valley width varies from 400 to 800 m. The floodplain width generally does not exceed first tens of meters, but increases to 800 m within three large icing fields, present at the valley bottom ( Figure 1).
Numerous icings appear annually in the Samokit River valley, merging by late winter into three large icing fields ( Figure 1). Their area varies from 0.2 to 1.6·10 6 m 2 , and volume, from 0.5 to 3.8·10 6 m 2 [43]. Icings originate from both open channel talik and sub-permafrost water [8]. Ice accumulation in the valley begins with the establishment of steady negative air temperatures [3], which in this region occurs around early November [13]. On average, the ice accumulation period lasts for 180 to 200 days, until early May, when the transition to positive air temperatures is observed and icings start to lose area and volume through thermal decay. Normally, Samokit River icings melt out completely by the end of June or the beginning of July, but Icing #2 meltdown may be retarded until late August, and small ice patches may persist throughout summer. The contribution of icings to the Samokit River annual runoff was estimated by SYGE around 40% [15,42,43]. Thus, groundwater and its seasonal redistribution play an important role in local hydrology.
This study site was selected because its groundwater icings are typical and representative for the region, and also because several previous studies were performed in this valley, describing icing development and its contribution to streamflow [13,43]. Icing #2 was chosen from three major icings at the Samokit River valley bottom for the detailed field studies. Similar icing areas are known in Russia (e.g., Verkhoyansk Range, Sayan Mountains, and others), Alaska (Brooks Range), and Canada (British Mountains), though at different elevations [3,5,6,12,14,15].
Several field-based techniques were employed during our study to obtain the icing morphometrical parameters, which allows for their intercomparison.

Materials and Methods
Icing #2 in the Samokit River valley was subject to detailed studies since autumn 2016. In September 2016, the staking grid was installed in the valley bottom, resembling previously used grid set up by the SYGE staff in 1970s and 1980s ( Figure 2). Then, the first UAV survey was performed in order to produce a baseline DSM of the ice-free valley bottom. Repeat observations on staking grid and UAV surveys were performed on 25 February 2017 and 2018. The fieldwork was conducted at the end of February due to complicated logistics in spring, around the dates of maximum icing extent and volume. In addition, satellite imagery was collected and used for estimating icing areas and volumes by using Sokolov's empirical approach [13], as described in Section 3.4.

Staking
Staking technique used in the present study is based on published guidelines [45,46] and is similar to the one used previously by SYGE (Figure 2, a) in the 1970s and 1980s, which allows direct comparison of the results obtained. Icing stakes are wooden rods with length from 2 to 5 m, dug vertically in the ground deep enough to ensure their stability, and labelled at 0.5 m intervals from the ground surface. When the stake is partially covered by icing, ice thickness is measured by a ruler from the closest non-covered label, and snow thickness is also recorded.
A staking network was set up in the Samokit River valley, within the Icing #2 limits, in September 2016. First, a reconnaissance survey was undertaken to plan the staking grid layout preliminarily, where the distance between stake transects were not exceeding one-sixth of the icing length, and between stakes within each transect, a quarter of the transect length, according to guideline [45]. At this stage, cryogenic landforms, i.e., frost mounds, and points of groundwater discharge in the valley bottom were mapped. Second, a staking grid was installed according to this preliminary design. Four staking transects were set up, with 400 to 500 m between them. Along each transect, five stakes were installed at 100 m spacing (Figure 2, b). Each stake was 4 m long, and labelled each 0.5 m. Third, a geodetic survey of the valley bottom was conducted during ice-free period, on 24 September 2016, to serve as a reference surface for further ice volume calculations, and several stakes were Global Positioning System (GPS) surveyed as control points (see Section 3.2 for details). In 2017 and 2018, ice thickness observations on the staking grid were performed on 25 February. If the icing surface was covered with snow, its thickness was also recorded. Based on ice thickness data from each stake, the icing volume was calculated according to guidelines [45], in two steps. First, the icing cross-section area under each transect was calculated ( Figure 2, b, transects II, III, V, VI). For each partial segment within a transect, segment width equaled spacing between the stakes, or, for the end segments, between the stake and the icing boundary ( Figure 2, c). Segment depth is average ice thickness between the neighboring stakes, calculated as arithmetical mean, or, for end sections, as half of the observed ice thickness at the outer stake. The same approach was then applied to calculate segment areas under individual longitudinal profiles, connecting stakes from different transects ( Figure 2, profiles a to e). From these data, partial ice volume was calculated for each block limited by individual transect and profile segments, and the sum of all partial volumes finally gives total icing volume.

Unmanned Aerial Vehicle (UAV) Surveys and Digital Surface Model (DSM) Creation
Aerial surveys were performed within the Icing #2 limits on 24 September 2016, 25 February 2017, and 25 February 2018. A Phantom 3 Advanced UAV (DJI, China) was used throughout this study. The average flyover height was 150 m, resulting in spatial resolution of 7 cm/pixel. Photogrammetric processing of the UAV survey imagery was done in Agisoft Photoscan Professional software v.1.4.5, and included survey alignment, alignment control using ground reference points, and dense point cloud creation. A digital surface model (DSM) and an orthomosaic with spatial resolution of 1.3 m/pixel were created for each survey date ( Figure 3). The ice-free DSM was further used as a reference, or subtraction, surface in icing volume calculations within the calculation boundary, i.e., within the actual icing areal extent on the survey date, from photogrammetric analysis. The spatial and vertical references of the UAV survey were controlled by using a reference point network established with two Ashteck ProMark 2 Global Navigation Satellite System (GNSS) receivers, working in static mode. The base station was installed on a Russian State Geodetic Network mark, and the other 12 ground control points were relatively evenly distributed within the icing ( Figure 3). In addition, several control points at the staking grid were surveyed, so that they could be further reused as common points for multitemporal surveys and subsequent analysis. Based on these common points, we checked that our future DSM was adequate in the long term, because other points were established on either a ground surface or an icing surface.
The point observation period varied from 20 to 40 minutes; ground control data were processed in GNSS Solutions v.3.80.8 software according to the manufacturer's guidelines [47]. DSM planar coordinate errors ranged from 2 to 51 mm (root mean square deviation of 0.005), with vertical reference errors varying from 3 to 100 mm (root mean square deviation of 0.012). Additionally, for processing of GNSS data, we used actual regional continuous GPS stations by the International GNSS Service (IGS) [48] for adjustment of triangulation. The data of actual regional continuous GPS stations were received from an open source for the period of our survey (25 February 2017 and 2018) [49]. The use of continuously operating reference station (CORS) data allowed excluding the potential bias from frost heave and other sources of long-term instability of the geodetic mark.
Icing boundary delineation, important for correct icing volume estimation, was performed manually for each winter survey date. Differentiation criterion between the icing and the surrounding terrain was a color of the orthophoto pixel, falling into three main gradations: (a) areas without icing, mainly shown in white shades (snow) or sometimes with dark shades corresponding to vegetation, (b) areas with current groundwater outpouring and icing accumulation are shown in yellowish-blue shades, and (c) areas having intermediate color between white and yellowish-blue. These latter shades correspond to areas where the icing accumulation was active in the last period, but had ceased since, and the freshly accumulated ice was covered with either falling or blowing snow. The icing surface in such areas is a slurry-like mixture of outpouring groundwater and snow, producing an intermediate color on the orthophoto. This color-based approach to icing boundary delineation is by far more precise as compared to linear approximation approach used in staking, thus the precision of the icing volume estimates is also expected to be significantly higher. A GPS survey might produce even more accurate delineation, but is more labor-intense and time-consuming, and is not expected to result in proportional increase in volume estimation accuracy. Thus, we expect that our approach to icing delineation yields sufficient accuracy for the Samokit River icing volume analysis.
Quantum GIS software v.3.4.8 was used to obtain the cumulative ice thickness distribution within the Icing #2 limits (Figure 4) from the ice-free DSM as reference, and orthomosaics for each survey date in 2017 and 2018 ( Figure 3). The icing boundary was traced for each year (see Figure 3), and total ice accumulation was calculated on a per pixel basis using this boundary as a mask.

Satellite Imagery Treatment
European Space Agency (ESA) Sentinel-2 optical imagery from and Sentinel-1 synthetic aperture radar (SAR) images were used in the present study, mostly because of open access to data, high image quality and acquisition frequency. ESA user products were acquired through the Copernicus Open Access Hub web interface [50]. Image treatment and analysis was done using ESA SNAP v. 6.0.10 and QGIS v. 3.4.8.
Icings are best detected via remote sensing when (a) they show the most contrast with the surrounding terrain, and (b) their area and volume are close to the annual maximum [46], which in the studied region is observed around mid-May [13]. At this time, ice bodies at the valley bottoms differ greatly from the surrounding areas, where the snow have already melted away. In early June, icings rapidly lose area and volume through thermal decay; this process can also be traced and quantified, but the annual maximum values of these morphometric parameters are the most demonstrative for long-term monitoring.
Satellite imagery was acquired for two temporal windows on each year: late February, i.e., close to field survey dates, and mid-May, the period of maximum ice accumulation. In the first case, the icing area was calculated from remotely sensed data to be used in comparative analysis against other techniques (Figure 5a-d), while in the second case, the reliability of Sokolov's equations was assessed (Figure 6a-f). In winter 2017, the closest optical image was from 6 March, and the scene was completely covered by clouds; while in 2018, the acquisition was close to field survey date and cloudfree. Spring images were acquired on 18 May 2017 ( Figure 6b) and 18 May 2018 (Figure 6a), under almost cloud-free conditions. Icing boundary delineation was performed using Sentinel-2 false color infrared (FCIR) band combination (RGB = 8-4-3), which produces highly contrasting images for ice surface and snow-free terrain (Figure 5a,b; Figure 6a,b). Also, FCIR combination allows color-based icing surface zonation as influenced by overflowing water (blue and green shades), or intact ice (white-blueish shades), snow over the icing or ground surface (pure white) and bare ground (red). The configuration of channels cut into the icing surface is also clearly detectable on these images.
The performance of the normalized-difference snow index (NDSI) in icing boundary extraction was also tested for mid-May images, since it was shown to be useful in ice cover studies [51]. NDSI was calculated using Sentinel-2 short-wave infrared (Band 11) and green (Band 3) channels with wavelengths 1.610 µm and 0.560 µm, respectively [52]. Icing surface shows NDSI values ranging from 0.80 to 0.99, falling close to a typical NDSI range for snow, and as such cannot be differentiated from snow-covered terrain [53]. ESA Sentinel-1 SAR imagery (Figure 5c,d; Figure 6e,f) was used in extracting icing areal extent and detecting groundwater outpouring over the icing surface. The acquisitions closest to the field survey dates (late February) and to available Sentinel-2 images (mid-May) were downloaded and post-processed. These images were acquired in the interferometric wide swath mode (incidence angle range 20°-46°) with VV (vertical-to-vertical) and VH (vertical-to-horizontal) polarizations [50]. The spatial resolution of the raw SAR data is 5 m in range and 20 m in azimuth, over a 250 km wide swath. SAR post-processing was undertaken in three steps [38]: radiometric calibration, multitemporal speckle filtering, and range-Doppler terrain correction. Finally, both VV and VH data were transformed from unitless intensity to backscatter reflectance, dB. Previously, the threshold backscatter values for RADARSAT and European Remote-Sensing Satellites (ERS-1, 2), used in icing boundary tracking, were discussed [5]. In southern Yakutia, the corresponding threshold value was -13 dB, and overflow water at the icing surface ranged below -18 dB.

Empirical Icing Volume Estimation
The first studies of the Samokit River icing morphometry and seasonal evolution were performed on Icing #1 by Sokolov in the mid-1960s [13]. His observations, employing staking technique, allowed documentation of temporal evolution of the icing, and based on these data, Multitemporal analysis of groundwater icing development, performed by Sokolov, led to two major conclusions. First, icing temporal evolution shows similar patterns regardless of geographical location. Icing area and volume at a given moment of time, expressed in relative values as percentage of annual maximum values, depend on the ice accumulation duration (Figure 7). This approach allows a first approximation of the annual maximum icing area and volume, using the values observed in the field at a given icing accumulation duration. In our case, by 25 February 2017 and 2018 the icing area was about 90%, average ice thickness, about 73%, and icing volume, about 64% of the annual maximum. These ratios may be further used to recalculate the annual maximum values of given icing morphology parameters. Second, icing volume depends almost linearly on icing area. Sokolov obtained a power law empirical relation, linking annual maximum icing area F, m 2 , and volume W, m 3 : W = aF b (1) where a = 0.96 and b = 1.094 [13]. This approach to calculations of icing volume at the end of winter is widely used in regional hydrological studies [25,54]. In the present study, Equation 1 was used to estimate annual maximum Icing #2 volume from satellite-derived areal extent of this icing in May 2017 and 2018 (see Figure 6a,b).

Ice Thickness Distribution
Ice thickness distribution in the Samokit River valley in winter 2017 and 2018 was estimated using two methods: staking (Table 1) and photogrammetry based on UAV surveys (see Figure 4). Icing #2 thickness is distributed unevenly within the icing field, with two maximums in the upper and the lower part of the icing, and only insignificant ice accumulation in the middle section is observed. The maximum ice thickness observed by staking in 2017 was up to 2.1 m in the upper part (transect II), and up to 3.0 m in the lower part (transect VI); in 2018, the same icing sections showed higher accumulation, up to 2.8 and 4.1 m, respectively (see Table 1).
Photogrammetry-based results (see Figure 4) were controlled to fit the observation data from the staking grid, so that during equalization, ice thickness values from the stalking grid nodes were used in DSM processing. Finally, the DSM product accuracy degrades with distance from the stakes, because of low vertical positioning accuracy of the UAV unit. Nonetheless, these spatially distributed data have enormously higher resolution, as they are based on per pixel analysis of orthomosaics. In general, maximum ice thickness observed with UAV-based photogrammetry is higher than obtained by staking, by 2 to 3 m. Maximum ice thickness reached 6.5 m (Figure 4) in the upper part of the valley, which is comparable with previously published values, up to 7 m [13].
Ice thickness distribution is patchy, with clearly detectable increased ice accumulation areas in the upper and lower parts of the valley. In the upper valley section, this area has circular shape and visibly correlates to a concave ice accumulation basin at the valley bottom, where locations of groundwater discharge points vary from year to year. A marked breakline is visible in the middle section of the valley (Figure 4a), probably reflecting differences in local hydrogeology and interaction with seasonal ground-freezing conditions. Numerous seasonal frost mounds were observed in the upper section of the valley, around transect VI (see Figure 2b). Frost mounds are oval-shaped, with linear dimensions about 5 × 10 m. In both observation years, their location within the valley varied only slightly. At the onset of melting period, 1 m deep voids were observed under the ice at these locations. The development of frost mounds evidences higher cryotic hydraulic head in the subsurface compartment (sub-icing talik), and hence higher seasonal freezing depth. This may further lead to increased groundwater discharge in the lower section of the valley, downstream the breakline, thus promoting ice accumulation.

Icing Areal Extent and Volume
Icing #2 area and volume were calculated using different data sources and at different points in time: late February and mid-May (Table 2), as explained in detail in the Methods section. Most values correspond to late February surveys, but several estimations of annual maximum values, derived using Sokolov's empirical approaches (see Section 3.4) and satellite imagery analysis, are also given for comparison.

Comparison of Methods: Icing Area
Icing #2 areal extent, observed by different methods, is comparable (see Table 2); still, there is substantial difference between values obtained from direct observations (staking) and indirect assessment (photogrammetry and satellite imagery). Taking staking as a baseline approach, relative deviation was between 14% (2017) and 17% (2018) for photogrammetry, and between 17% (2018) and 24% (2017) for satellite data. Inverse interannual difference in satellite-based estimates, where Icing #2 area was larger in 2018 than in 2017, may reflect the uncertainty in SAR-based icing boundary delineation, which is clearly visible upon direct comparison (see Figure 5d, f). We also suppose that the staking-based icing area is underestimated, because this method only allows linear interpolation of the icing boundary between neighboring transects and does not account for real icing geometry.
Maximum areal extent of the Icing #2 was also assessed from satellite imagery (see Table 2). Only slight increase was observed, compared to late February values, which means that in the second half of winter, ice accumulation mostly occurs at the icing surface, and not along its boundaries. By late February, the Icing #2 reaches about 93% to 94% of its maximum linear dimensions, which is close and insignificantly higher than 90%, previously published by Sokolov [13]. This may indicate either the difference in methods, by which these values were observed, or a slight shift in icing areal growth rate throughout winter. In the latter case, climate may be a key driver of changes in icing evolution, through variations in seasonal freezing depth and water redistribution in the subsurface. If this is the case, then Sokolov's findings and equations (see Section 3.4) are to be revised; however, there are still insufficient data for definitive conclusions on the matter.

Comparison of Methods: Icing Volume
Staking, again, is the only method providing direct and most conservative estimates of the Icing #2 volume on observation dates in late February (Table 2). Photogrammetry-based estimates are significantly higher, by 142% in 2017 and by 74% in 2018. These estimates then were corrected for snow cover, since UAV-based DSM surface corresponds essentially to snow cover surface, and not the icing surface. Snow covered about 72% and 51% of icing areal extent, with an average snow cover thickness about 16 cm and 10 cm in February 2017 and 2018, respectively. This corresponds to an excess icing volume of 0.18 × 10 6 m 3 (2017) and 0.08 × 10 6 m 3 (2018), or less then 3% of the total icing volume. Finally, snow-corrected Icing #2 volume equals 4.02 × 10 6 m 3 (2017) and 3.81 × 10 6 m 3 (2018), or by 132% (2017) and by 71% (2018) higher, than obtained by staking.
This difference originates, in our opinion, from the inability of our staking grid setup to capture variations of ice thickness within the icing. First, the number of staking profiles is insufficient and should be at least doubled to properly account for icing bed topography and variable ice accumulation across the Samokit River valley bottom. Second, our staking data appear to be negatively biased. The highest ice accumulation corresponds to the sections of the anabranching Samokit River channel, but our stakes were mostly installed on alluvial islands, separating numerous River channels, to secure the stakes' stability. The Samokit River depth locally reaches from 1.5 to 2.0 m, hence our staking data underestimates the ice accumulation in the channel zone by roughly this order of values.
Low spatial-positioning quality of the Phantom 3 UAV could also affect the resulting icing volume estimates, but it is yet impossible to adequately quantify this effect.
Satellite-based mid-May Icing #2 area was used to calculate its volume using Equation 1, and photogrammetry-based data from late February (see Table 2) were recalculated from mid-winter to maximum annual values using the ratio given in Sokolov's publications (see Section 3.4 and Figure  7). These two sets of values are coherent, although the difference is much higher for 2018 than for 2017. Considering the technical difference in these approaches, we judge the degree of convergence quite satisfactory, and overall supportive of the Sokolov's approaches.
Technical progress during recent decades lead to the development of new techniques for icing morphometry observations. These techniques are nonetheless associated with various sources of uncertainty; quantifying these sources, even at a local scale, is an important task.

Icing Volume Data: Sources of Uncertainty
Staking-based observations are direct, in the sense of field ice thickness measurements, but these results are highly uncertain. First, ice thickness values are interpolated between adjacent stakes and neighboring transects, neglecting potential influence of valley bottom topography, especially at the end sections of transects. Ice accumulation is unevenly distributed across the icing surface, and for a given year, high accumulation zones may or may not coincide with the staking grid points or transects, so that the calculation accuracy may vary correspondingly. Second, linear interpolation of the icing boundary is used between the same end sections of each transect, while remotely observed icing boundary is highly curved. Third, as is the case of our study, most stakes were set up at alluvial islands and on the floodplain, since it is almost impossible to ensure the staking rod stability installed in the channel. The highest ice accumulation, however, is observed within the deepest parts of the channel. As a result, staking technique tends to underestimate the icing parameters.
The photogrammetric approach based on UAV surveys and inferred DSMs appears to yield sufficiently precise estimations of icing area and volume. However, our field experience showcased several important complications associated with this technique, potentially affecting its accuracy. Low negative air temperatures significantly reduced battery life of the UAV, and large aerial extent of the surveyed area in the Samokit River valley, about 2.1 km 2 , required several, usually three to four, flight sessions, interrupted by battery replacement. In its turn, battery replacement was associated with significant positioning bias, up to hundreds of meters, between the flight sessions. This complicated subsequent data binding during photogrammetric processing. Moreover, an adequate number and distribution of ground control points (GCPs) is essential in minimizing this positioning bias. For example, for the surveyed area, our 12 GCPs were not enough to significantly improve the lateral and vertical accuracy of the DJI Phantom 3 UAV data in post-processing. The DEM equalization procedure is effective within first tens of meters around the GCPs, while in our case the distance between them was of the order of first hundreds of meters.
DSMs, inferred from the UAV surveys are not void of uncertainties associated with either trees and shrubs in the valley bottom or snow accumulation over the icing body. The former can be neglected because of scarce presence of vegetation in the valley bottom; the latter can significantly influence the estimation accuracy.
Falling and/or blowing snow covers the icing surface in winter, and can either be melted by outpouring groundwater and incorporated into the icing volume, or stay intact, in which case it directly affects the icing volume estimations using photogrammetry. According to our field observations, snow accumulation does not exceed 30 cm at the icing surface, and can reach 50 to 60 cm in areas adjacent to the icing boundary. Assuming the average Icing #2 area about 1.61 × 10 6 m 2 , an average snow cover thickness of 16 cm, and its uniform distribution over the entire icing surface, then the associated snow volume would be 0.24 × 10 6 m 3 , or approximately 6% of the icing volume. Normally, the extent of snow cover reaches 72% of the total icing area. Thus, the associated error is generally <4%-5%.
Satellite images, used for icing morphometry estimation, also have their limitations. First, clouds on visible spectrum images may partly or totally cover the area of interest. Second, both overrunning surface water and snow cover complicate icing boundary delineation. NDSI mapping is of limited utility, since its values for snow cover are close to that of the icing surface.
Cryogenic processes transform the icing surface and body; multiple frost mounds were observed at the Icing #2 surface in 2017 and 2018. Frost mounds, as we observed in the field, are associated with voids in the icing body, and the presence of such voids may influence the volume estimate accuracy. The unit volume of one void can reach 50 m 3 , so their total volumetric effect is in the order of 250 to 500 m 3 , or less than 0.01% of the total icing volume, and in our case can be neglected.
Every spring during the freshest period, the snow meltwater erodes both the icing surface and the channel bed. In the latter case, if only one DSM is used in the volumetric analysis, uncorrected for these deformations during repeat ice-free surveys. Our visual observations in 2016-2018 indicate no significant changes in the topography of the valley bottom and river channel. By the end of the spring freshet, we observed grass cover and shrub disturbance, and accumulation of sand alluvium sediments up to 10 m long, 1 to 2 m wide, and 0.3 to 0.5 m thick, with total volume not exceeding 10 m 3 , several orders less than the icing volume.
Sokolov's equations, i.e., Equation 1, can be applied when the maximum areal extent of icing is known, but can yield large errors. First, these equations were derived from field data of winter 1964/1965, and may give erroneous results for other weather conditions and physiographical settings. Second, icing volume estimates can be imprecise; since an increase in ice accumulation does not necessarily result in larger areal extent, undermining the efficiency of Sokolov's Equation 1.

Conclusions
A giant Icing #2 was studied in the Samokit River valley, southern Yakutia, and its morphometry, i.e., icing area and volume, was assessed using several widespread techniques: staking, UAV-based photogrammetry and satellite imagery analysis. These techniques generally converge in mid-winter icing area estimates, ranging from 1.35·10 6 m 2 to 1.68·10 6 m 2 , but significantly diverge in icing volume, varying from about 2.0·10 6 m 3 to about 4.0·10 6 m 3 , obtained by staking and photogrammetry, respectively. The end-of-spring, or annual maximum icing area and volume, were calculated from satellite imagery analysis and Sokolov empirical approaches; the areal extent is just slightly higher, than in mid-winter, about 1.76·10 6 m 2 , while the volume increases to about 6.3·10 6 m 3 , averaged across methods and years.
Intercomparing widespread techniques of icing morphometry estimation highlights their advantages and drawbacks. Staking is the only direct method used in our assessment, based on in situ measurement of both ice and snow cover thickness, and is regarded as a baseline tool. Once set up, the staking grid can be further reused, providing uniform results with comparable accuracy. This method, however, tends to underestimate icing dimensions, especially icing volume. It needs a thorough pre-implementation planning, to ensure that staking grid covers all representative icing sections.
The photogrammetry technique based on UAV survey data is promising, owing to more accurate DSM-based per pixel calculations of icing volume, rather than point-based staking. The derived surface model can be further reused in other studies concerning the river valley topography, i.e., in fluvial geomorphology. The shortcomings of this approach, however, are also numerous.
High-precision UAVs have significantly higher price, and cheaper options do not necessarily possess sufficiently long autonomous flight time, and adequate image and positioning quality, especially in a vertical dimension. The survey precision further degrades when numerous overlapping flyovers are performed for a single object, forced by low battery capacity at negative air temperatures. Correction for these effects require an extensive network of ground control points, surveyed with high-precision GNSS receivers and used in the DSM equalization. The resulting DSM is also not void of inherent errors, since this is a surface model, and not a terrain model. The effects of vegetation and snow cover can lower the estimation accuracy, and require additional snow cover surveys.
A significant difference between staking-and photogrammetry-based estimates, about 2.4 times, in our case originates from insufficient number of staking transects and their distribution across the icing surface, overrepresenting alluvial islands and lacking data points from the river channels with higher ice accumulation. On the other hand, unsatisfactory DEM equalization by landmarks under photogrammetry post-processing, owing to a low number of GCPs, precludes us from claiming this technique as highly accurate.
The utility of satellite imagery, both optical and SAR, is limited by the scarcity of cloud-free scenes (Sentinel-2) and ambiguity in icing boundary delimiting (Sentinel-1). Nonetheless, satellitebased icing area estimates are close to those obtained by other techniques, and the potential of this approach, especially based on the interpretation of SAR acquisitions, merits further studies. It may significantly limit the presence in the field.
As a result, icing areal extent is estimated with higher precision than icing volume. The presented data and analysis are the first approximations, and this research will be continued. First, an UAV unit with higher positioning accuracy is expected to increase our DSM quality. Second, satellite imagery analysis will be enhanced in order to develop an adequate technique for icing surface zonation and icing boundary delineation, where SAR data would play a crucial role. The interferometric SAR approach using proprietary L-band SAR data, is also a prospective approach to test against our field data.