Ice Roughness Estimation via Remotely Piloted Aircraft and Photogrammetry

Structure-from-Motion Photogrammetry conducted with images obtained via Remotely Piloted Aircraft (RPA) has revolutionized the field of land surface monitoring. RPA-Photogrammetry can quickly and easily capture a full 3D representation of a study area. The result of this process is a high-definition Digital Elevation Model (DEM) representing the land surface of a given study area. It is particularly useful in applications where land surface data collection would otherwise be expensive or dangerous. The :::: The monitoring of fluvial ice covers can be time-intensive, dangerous, and costly, if detailed data are re5 quired. Fluvial ice roughness is a sensitive parameter in hydraulic models and is incredibly difficult to measure directly using traditional field methods. This research hypothesized that the surface roughness of a newly-frozen fluvial ice cover is :: Ice :::::: covers :: on : a ::::: river :::::: surface ::::: causes ::::::::: resistance :: to ::::: water :::: flow, ::::: which :::::::: increases :::::::: upstream ::::: water ::::: levels. ::: Ice :::: with :: a ::::: higher :::::: degree :: of ::::::::: roughness :::: cause :::::::: increased :::: flow ::::::::: resistance ::: and :::::::: therefore :::: even :::::: higher :::::::: upstream :::: water :::::: levels. :::::: Aerial :::::: images ::::::: collected ::: via ::::::::: Remotely :::::: Piloted :::::: Aircraft :::::: (RPA) :::: were ::::::::: processed :::: with :::::::: Structure :::: from ::::::: Motion :::::::::::::: photogrammetry :: to ::::: create :: a :::::: Digital :::::::: Elevation :::::: Model :::::: (DEM) :::: and 10 ::: then :::::::: produce :::::::::: quantitative :::::::::::: measurements :: of ::::::: surface ::: ice ::::::::: roughness. ::::::: Images ::: and ::::::: surface ::: ice :::::::: roughness :::::: values ::::: were :::::::: collected ::: over :::: two ::::: years :: on ::: the ::::::: Dauphin ::::: River :: in ::::::::: Manitoba, ::::::: Canada. :: It ::: was ::::::::::: hypothesized :::: that :::::: surface ::: ice ::::::::: roughness ::::: would :: be : indicative of subsurface roughness. The :: ice ::::::::: roughness. :::: This : hypothesis was tested through a comparison of ice roughness determined through the statistical analysis of RPA-photogrammetry DEMs to ice roughness values :: by ::::::::: comparing :::::::::::: RPA-measured ::::::: surface :: ice ::::::::: roughness ::::: values :: to ::::: those predicted by the Nezhikhovskiy equation. The Nezhikhovskiy equation is a widely used empirical 15 method for estimating ice roughness based on observed , :::::::: wherein ::::::::: subsurface ::: ice :::::::: roughness :: is ::::::::::: proportional :: to : ice thickness. Hydraulic and topographic data were collected over two years of field research on the Dauphin River in Manitoba, Canada. Various statistical metrics were used to represent the roughness ::::: height : of the DEMs. Strong trends were identified in the comparison of ice cover roughness values determined through RPA-photogrammetry and those calculated via :::::::::::: RPA-measured ::: ice :::::: surface :::::::: roughness ::: to ::::::::: subsurface ::: ice :::::::: roughness :::::: values :::::::: predicted ::: by the Nezhikhovskiy equation, as well as with :::::::::: comparisons 20 :: to ice thickness. The :::::: standard :::::::: deviation :::: and inter-quartile range of observed roughness heights was :::::::: roughness :::::: heights ::::: were determined to be the most representative roughness metric. The maximum peak value performed better in some cases, but the fact that this metric would be heavily influenced by outliers led to it being rejected as a representative metric. Three distinct forms of surface ice roughness were noted: rough, smooth, and ridged. Statistical :::::::: statistical :::::: metrics :::: and :::::: several : properties of the DEMs of fluvial ice covers were calculated ::: and :::::::: observed. No DEMs were found to be normally distributed. k-means 25 clustering analysis was used to group sampled data into two categories, which were interpreted as rough and smooth ice . The inter-quartile range of the smooth and rough categories were found to be 0.01 0.05 meters and 0.07 0.12 meters,

This relationship is based on measurements conducted on rivers in Russia several decades ago and it has served well as an estimation tool for engineering applications. Using more complex data, Equation 1 was adapted by Beltaos (2013) for use in the estimation of the roughness of newly-formed ice jams, resulting in Equation 3.
The value given for c n = 0.095 has been determined to be representative for ice jams. Additionally, the inclusion of the 160 hydraulic radius R[m] accounts for the fact that the roughness elements of ice jams are often of such magnitude as to increase relative roughness to the point where it has significant impact on the hydraulic radius. This relationship is only valid for newly-formed ice jams. Immediately : ; ::::::::::: immediately after formation, the ice is subject to shear forces from the water flowing underneath, which slowly smoothens the sub-surface of the ice cover (Ashton, 1986).

175
RPAs were even used by Alfredsen et al. (2018) in the mapping of river ice in Norway.
There are several methods of producing a DEM from a set of overlapping close-range aerial images. Currently, the most widely used and time-effective method is that of Structure from Motion (SfM) photogrammetry (Matthews, Fraser). For the purposes of this research, SfM photogrammetry was the sole method employed. An important consideration when using the RPA-photogrammetry method for topographic data collection is the impact of doming errors on the final product (James) .

180
Doming errors are most prominent when all images are taken from a parallel axis (Eltner and Schneider, 2015) . In the case of RPA-photogrammetry, this is when the camera angle is set to 0°tilt, however, some distortion is also caused by the shape of the camera lens. Most advanced software packages used to produce DEMs from photogrammetry data include a self-calibration process that develops a model of the distortion caused by the lens of the camera. Eltner et al. (2016) makes the distinction between local surface quality and more systematic errors such as doming, relating these two categories to the precision and 185 accuracy of the DEM, respectively.
After targets were placed, their centres were surveyed using a Leica Viva GS14 ® survey-grade Real-Time Kinematic (RTK)

205
Global Navigation Satellite System (GNSS) base-and-rover system, which is typically observed to have an in-field reported horizontal error of ≈ 2 cm, and a vertical error of ≈ 3 cm. The Canadian Geodetic Vertical Datum of 2013 (CGVD2013) geoid was used in the recording of all surveyed elevations. Localization was assessed using a Manitoba Infrastructure (MI) benchmark located near DRLL03, and verified using the Natural Resources of Canada Canadian Spatial Reference System Precise Point Positioning (CSRS-PPP) service. Further benchmarks were established using the CSRS-PPP and"leap-frogging" 210 to further benchmarks. In the 2019-2020 season, some RPA flights were completed without targets, to allow for more flights to be completed during the field visit. A comparison between the representative metrics calculated from a DEM with and without geo-rectification is presented in Section 4.1. and 230 m across the river. Ideally, geo-rectification targets used in RPA-photogrammetry would be evenly distributed across the study area. In this research, targets could only be placed on the left bank of the study area, and none could be placed directly on floating ice, due to safety concerns.
During the 2019-2020 field season, the RPA mission planning application Pix4Dcapture ® was used to plan and automate RPA flights over study areas. This greatly reduced the required flight time, and produced similar, if not better photo coverage.

Field Accuracy Tests
There was a need to quantify the impact of the ground control targets being clustered :::::: required :::::: ground ::::::: control ::::::: grouping : on the left bank of the study area. The field methods described in Section 3.1.1 were repeated at River's Edge Nursery in La Barriere, Manitoba. A fully dry land study area of equivalent size to typical study areas flown at the Dauphin River was delineated, and 15 targets were distributed. The targets were conceptually grouped into three areas: typical, middle, and end. The "typical" 255 group represented a target distribution that was generally produced during field work at the Dauphin River sites. The "middle" and "end" target groups were supplemental, which would be added or subtracted from the photogrammetry analysis to test their respective impacts on DEM accuracy. The distribution of targets in the study area is represented in Figure 4. Finally, after the RPA flight was conducted, 10 independent and unmarked locations were captured by RTK-GNSS survey as a check for accuracy in subsequent data analysis. The :::::::::::::: photogrammetry ::::::::: processing :::::::: software ::::::: selected ::: for :::: this ::::: study :::: was : PhotoScan Professional ® from AgiSoft LLCuses the SfM method of photogrammetry, which has been widely used in geosciences in recent years (Westoby) , as well at in this study. Gini et al. (2013) compared their custom research-grade photogrammetry algorithms to results obtained from Pix4UAV 265 Desktop ® and PhotoScan Professional ® . Their findings suggested that these commercial packages performed similarly to their software, with PhotoScan Professional ® performing somewhat better than Pix4UAV Desktop ® . PhotoScan Professional ® is also considered to be a relatively fully-featured and complex (Eltner and Schneider, 2015) tool as compared to other options, and includes an automated process for estimating and correcting doming errors.
Images were imported into PhotoScan Professional ® , and aligned to create a sparse point cloud of tie points. Where targets 270 were used, they were identified in all images containing them, and their coordinates were imported to provide geo-rectification of the resultant point cloud. A dense point cloud was then generated, followed by a DEM. An example point cloud consisting of ≈ 15 million points and corresponding DEM are shown in Figures 5 and 6.

Accuracy Testing
The impact of placing all control points on only one bank :: an ::::::: extreme ::: end ::: of ::: the ::::: study :::: area : was tested through a detailed 275 trial on an open field. Groups of targets were used as input to the photogrammetry software and the resultant DEMs were  compared to the 10 independent survey points. The DEM generated using all available targets was assumed to be the most correct representation of the land surface, against which all other target groupings were compared. A maximum acceptable vertical difference of 0.03 m between the DEMs and the independent survey points was adopted. This value was chosen based on the typical error observed in the data gathered by the RTK GNSS base-and-rover system. This system was the limiting 280 factor for accuracy in this study since it was the tool which informs the absolute spatial position of all field equipment. The following target scenarios were tested: "all points" utilizing all ground control targets, "typical points" using all the targets identified in the "typical" subset, "three points" using a subset of three targets from the "typical" subset, and "two points" using a subset of two targets from the "typical" subset. In the two and three point tests, the most spatially distributed targets within the "typical" subset were selected. An additional test was required to determine if systemic errors were introduced in 285 DEMs generated without the use of geo-rectification targets. The data :::: This ::::::: scenario :: is ::::::: referred :: to ::: as ::: the :::: "no :::::: points" ::::: case. :::: Data collected at site DRLL06 on 2019-11-13 was prepared with and without the inclusion of control point data. A maximum acceptable percent error of 5% was adopted to evaluate the results of this test.

Roughness Characterization
To avoid unwanted influences in the surface slope and texture, a 50 m 2 sub-sample from the center of the river was taken, 290 which excluded all overbank objects and sections of the ice cover that were near to the bank. Additionally, a three-dimensional plane-of-best fit Linear Model (LM) was found for each sample, and then subtracted from the surface data. The goal of this was to normalize each data set, setting the average surface elevation to 0, and removing the river slope from the sample data. Gadelmawla et al. (2002) noted that the average surface elevation is the most commonly used, and most sensible reference standard from which to assess roughness height. By shifting the elevation data down to a base elevation of 0, and removing 295 unwanted patterns, each data point was transformed from an elevation to a roughness height.
A two-dimensional Fast Fourier Transform (FFT) was then applied to each sub-sample, with the goal of filtering the input data and removing other surface trends beyond those addressed with the plane-of-best-fit. The combination of the LM and FFT adjustment and filtering process will be referred to as LMFFT. Though an analysis of dominant frequencies it was found that the lowest frequencies (< 1 m −1 ) had the largest amplitudes, while the highest frequency signals (> 5 m −1 ) had the 300 lowest amplitudes. A typical distribution of amplitudes is shown in Figure 7 a. These results are arranged such that the highest frequency waves are found in the center of the matrix, while the lowest frequencies are found in the corners of the matrix.

RPA Performance
The RPA chosen for this study , the DJI Phantom 4 Professional ® , performed well during all field visits in various weather and cloud conditions. In extreme cold (-20°C and below) the RPA performed all functions well, although the battery life was reduced by approximately 50%. Additionally, it : It : was found that the RPA had to be powered on in a warm area, such as the 370 heated cab of the field vehicle. Once powered on, it could then be placed outside for normal operations.
RPA-Photogrammetry performed very well across a variety of scenarios in the land-based field accuracy tests. The worst observed :::: DEM : accuracy was found in the "two points" scenario with an average vertical difference of 6.30 m calculated between the DEM and the 10 independently surveyed test points. The "three points", "typical points", and "all points" scenarios all had the same average vertical difference of 0.03 m. These three scenarios were all within the acceptable vertical difference 375 criteria of 0.03 m identified in Section 3.2.2. The maximum error of the DEM was determined to be the sum of the expected vertical accuracy of the RTK-GNSS unit, 0.03 m, and the observed vertical difference difference between the DEMs and surveyed locations, 0.03 m, for a total expected error of 0.06 m. The results of the test ::: "no :::::: points" :::: test :::: case, :::::: which :::::: sought to determine if a lack of geo-rectification targets introduces systemic errors into the DEM data : , are presented in Table 2. prior years lower flows were noted. This resulted in notably thinner ice covers, more thermal ice growth, and less extensive ice Table 2. Percent and absolute difference between various statistical metrics computed from geo-rectified and non geo-rectified ::: ("no ::::::: points") DEMs produced from data collected and site DRLL06 on 2019-11-13

Dauphin River Results
Max. Min.

Statistical Properties of Ice Roughness Height Distributions
Several different forms of ice roughness were observed in DEMs produced using RPA-photogrammetry. Figure 8 illustrates three different ice roughness forms ::::::: observed :: at :: a ::::: single :::: site, and their appearance in cross section along the indicated transects.
Cluster analysis was conducted using k-means clustering. The optimal number of clusters was found to be two using average silhouette analysis. This was enforced by observations taken at the time of sample collection. The samples were broadly 410 separated into two groups corresponding to the visual extent of "rough" or "smooth" ice, as defined earlier in this section. All observations at site DRLL06 were found to belong to cluster "1", while observations at site DRLL03b, DRLL05, DRLL08, and DRLL08a belonged to cluster "2". The within cluster sum of squares was found to be 5.44 and 11.96 for clusters 1 and 2 respectively.
This limitation is well known within photogrammetry (Harwin, Hamshaw,Lane) , but since there were no trees in the areas of interest of our study areas, this limitation was not consequential to the study.
As described in Eltner et al. (2016) systemic errors causing deficiencies in DEM accuracy differ from local-scale errors 425 causing deficiencies in DEM precision. Systemic errors include those incurred by improper sampling technique and by limitations in the analysis. These errors were largely assumed to have been handled to the extent that is possible in this research by the automated processes in AgiSoft PhotoScan Professional ® . Eltner and Schneider (2015) found that AgiSoft PhotoScan Professional ® also performed well in reproducing the texture of complex natural surfaces. Direct comparisons could not be made in this research between the naturally occurring ice surfaces and the RPA-photogrammetry reproductions. However, the 430 magnitude of such results as the maximum ice peak value matched visual observations and field notes. The accuracy test performed at the La Barrier field site confirmed that with appropriate ground control points this method could accurately reproduce snow-covered land surfaces. It also showed that the method could precisely reproduce features of the same order of magnitude as the chunks of ice expected to be measured on the Dauphin River.
Some of these differences could be explained by differences in light conditions between the RPA flights, indeed this may have given rise to some of the differences observed at site DRLL05 over the same time period. This was corroborated by the RPA flight conducted at DRLL06 on 2019-11-26, under similar light conditions to the flight on the 13 th , which showed very similar general and peak distributions. Similar general and peak distribution data were noted at site DRLL08a between the two RPA 470 flights conducted on the 13 th and 2019-11-26. The raw distributions had some noted differences in composition, although they were similar in height and width.
Using the results of the k-means clustering analysis, a breakdown of ice roughness cluster characteristics was generated. Table ?? illustrates the site grouping under the cluster analysis results, with maximum, mean, and minimum values reported within each cluster. The two best performing roughness metrics in the above subsections, IQR and 84 th percentile of peaks are 475 included, as well as the strongest defining metrics in the principal components analysis, IQR and kurtosis are included in this table. This table may serve as a defining guide for further categorization of ice surface types. Since cluster "1" corresponded to ice surfaces with higher IQR and 84 th percentile of peaks values, this group was interpreted as the "rough" ice group.
This analysis could be extended to examine shear walls, and determine maximum ice elevation of a freeze-up jam, after the fact.
Competing interests. The authors declare that they have no conflict of interest.