Long-term mass-balance monitoring and evolution of ice in caves through structure from motion–multi-view stereo and ground-penetrating radar techniques

Introduction

Ice caves are defined as permanent (>2 years) ice deposits in caves or lava tubes (Perşoiu and Onac, 2012) and for this reason are recognized by most authors as an evidence of permafrost (Colucci and Guglielmin, 2019; Holmlund et al., 2005; Luetscher, 2013; Luetscher et al., 2005). Although cave climate is rather complex and allows conservation or increase of ice even in environments with mean annual air temperature (MAAT) well outside the typical values of periglacial environments (e.g. Holmlund et al., 2005; Colucci et al., 2016), it has been recently observed that anthropogenic climate warming is also drastically affecting this portion of the cryosphere (e.g. Kern and Perşoiu, 2013; Colucci and Guglielmin, 2019). Ice bodies in caves are suffering significant mass losses worldwide, and it is therefore urgent to study these sources of paleoenvironmental information that will soon disappear (e.g. Kern et al., 2009; Colucci et al., 2017; Sancho et al., 2018).

Long-term mass balance measurements are rather common on glaciers both in mountains and polar environments (e.g. Marcer et al., 2017), and in the last decades, new methodologies such as LiDAR, terrestrial/aerial photogrammetry and pseudo 3D ground penetrating radar improved the reliability of the surveys (Berenguer et al., 2014; Forte et al., 2014; Mertes et al., 2017; Piermattei et al., 2016; Santin et al., 2019). This is not the case for permanent ice deposits in caves where only scattered and short-term mass balance records exist (Bertozzi et al., 2019; Obleitner and Spötl, 2011; Perşoiu, 2018). Moreover, in all such cases, mass balance measures refer to specific discrete locations as for instance in the Eisriesenvelt ice cave (Austrian Alps) where an ultrasonic ranger continuously monitored ice thickness changes (Obleitner and Spötl, 2011), or in the Leupa ice cave (SE Alps), where surface topography changes have been monitored by using 2 fixed benchmarks since 2011 (Colucci et al., 2016; Colucci and Guglielmin, 2019).

SfM-MVS is not a stand-alone technique but rather a workflow combining different algorithms coming mainly from computer vision, dense image matching and partially from photogrammetry (Anderson et al., 2019; Smith et al., 2016). The used algorithms differ from traditional photogrammetry because camera positions, orientation and the entire scene geometry are retrieved simultaneously without any georeferencing (Smith et al., 2016). However, the acquisition does require references for scaling such as ground control points (GCPs) or even the camera positions to improve the overall accuracy of the final output, making it a reliable base for further spatial analysis (e.g. Sanz-Ablanedo et al., 2018; Anderson et al., 2019). Nevertheless, GPCs or camera positions are needed only for scaling purposes while a relative DEM can be produced without. The proposed methodology can yield a resolution comparable with those achieved by using the terrestrial laser scanner (TLS) in the creation of a digital terrain model (DTM; Piermattei et al., 2016). Moreover, the approach reduces the instrumental cost making it easier to improve temporal resolution by repeating surveys, which is fundamental in multi-temporal analysis such as long-term mass balances. Another advantage is the applicability across a range of scales, making it possible to focus on small details (e.g. bédière, sediments and cryoconites) or the whole cave area without altering the processing flow (e.g. the approach is reliable from 10⁻²–10⁶ m²; Smith and Vericat, 2015; O’Connor et al., 2017). The easy data acquisition and the lightweight instrumentation are even more useful in particularly remote and not-easily-accessible environments such as ice deposits within caves. To our best knowledge, the SfM-MVS method has never been used in ice caves for detecting multi-annual ice surface topography variation and mass balance calculation. There are very few examples exploiting laser scanning data, although limited to short periods of acquisition (e.g. Šupinský et al., 2019) and just considering topographic variations as they do not include GPR datasets.

GPR is a widespread technique in several glaciology-related surveys; thanks to the dielectric properties of frozen materials and ice which both exhibit low intrinsic attenuation and, therefore, the EM waves can reach depths conspicuously higher than in common geologic materials (e.g. Woodward and Burke, 2007; Navarro and Eisen, 2009). Other key aspects that are making GPR a widespread solution for glaciological surveys are the relatively fast acquisition times and the high lateral and vertical resolutions, which are in the order of some centimetres, at least for the shallower applications. Recent application of GPR techniques on ice bodies in the Julian Alps area are numerous and usually related to glaciers (e.g. Forte et al., 2014; Colucci et al., 2015; Del Gobbo et al., 2016) while there are only a few examples around the world of GPR in underground ice bodies (e.g. Hausmann and Behm, 2011; Colucci et al., 2016; Lende et al., 2016).

The aims of this paper are to (i) quantify volume and surface changes of permanent ice deposits in selected caves to calculate seasonal and annual mass balance from SfM-MVS point clouds; (ii) estimate the total volume of ice present in each deposit through GPR survey data; (iii) develop, test and validate a repeatable and relatively fast workflow to monitor the evolution of the underground cryosphere over time. This would in turn (iv) allow a more accurate estimation of water resource stored in the underground cryosphere of high elevated karstic areas, which represents a potential secondary source of water supply often underestimated.

Study area

Vasto and Leupa ice caves are both located on the northern side of the Mt. Canin–Kanin massif (2587 m; Julian Alps, N-E Italy, 46°21′ N, 13°26′ E) (Figure 1) and contain a well-documented permanent ice deposit (Bertozzi et al., 2019; Colucci et al., 2014, 2017; Colucci and Guglielmin, 2019). The two caves are characterized by different cave air dynamics (CAD) and ice accumulation mechanisms. The lithologies of the Canin massif are mainly limestone and dolostone from the upper Triassic. Particularly, stratified and faulted Dachstein limestone is highly karstified and permeable with a reported number of caves close to 1000 (including the Slovenian side of the massif; Colucci et al., 2016). On the opposite, dolostones (i.e., Dolomia Principale formation) exhibits secondary porosity, and to some extent, might be considered as a basement for groundwater drainage given the presence of a high permeability contrast (Turk et al., 2015). From the geomorphological point of view, the area is characterized by the combination of glacial and karstic phenomena typical of high mountain karst terrain where glaciokarst is widespread (Žebre and Stepišnik, 2015). Mean annual precipitation in the Julian Alps and Canin massif is the highest in the Alps (Isotta et al., 2014); the 1981–2010 annual average was around 3300 mm (Colucci and Guglielmin, 2015; Crespi et al., 2018). Precipitation influences mean winter snow accumulation, which at an altitude of 1830 m a.s.l. equates to about 7.0 m (Colucci and Guglielmin, 2015). The 1981–2010 mean annual air temperature was 1.1 ± 0.6°C.OPEN IN VIEWER

The main entrance to Vasto ice cave (VIC) is located on the north side of the Canin (a.k.a. Kanin) massif at 2162 m a.s.l. The area close to the cave entrance hosts several small glaciers, glacierets and ice patches representing some of the lowest evidence of glaciation in the European Alps (Colucci and Zebre, 2016; Colucci et al., 2021). The presence of steep rock walls and a deep karstic gorge opening to the north favour wind drifted snow accumulation right next to the large entrance of the cave (15 x 4 m approximately). A perennial firn ice cone partially filled by small avalanches or debris coming from the surrounding walls is located here. Even though the VIC main entrance is quite large compared to other typical ice caves, during very wet winters (e.g. 2008–2009 and 2013–2014), the snow can almost entirely block the entrance. The presence of liquid water inside the cave is a frequent phenomenon in summer. Although this water is almost always drained downwards, there is a time lag that generally coincides with the beginning of the ablation season, where a temporary surface pond formation is observed (on rare occasions, e.g. 2020 season, the pond was observed till the end of August). This accumulated water, which is mainly due to melting snow and rain, can exceed several decimetres of depth and it is usually abruptly drained at some point between late June and late July through lateral conduits in contact with the hosting rock (Colucci et al., 2016). The primary ice accumulation mechanism inside the VIC is represented by snow metamorphism (i.e. firnification) combined to a lesser extent with the refreezing of water dripping from the cave ceiling. CAD is rather complicated due to the strong interactions with the external environment, having such a large main opening combined with a lateral chimney open or closed depending on the amount of snow. In detail, when the chimney is filled by snow, the CAD results in a static ice cave behaviour. During summer and autumn, when the snow partially melts or during very dry seasons when there is not enough snow to block the chimney, the circulation is dynamically active and driven by density flows (Bertozzi et al., 2019; Colucci et al., 2016).

Leupa ice cave (LIC) entrance is located at 2285 m a.s.l. on the north side of Mt. Leupa (2402 m a.s.l.). The cave is characterized by the presence of two or more openings at different altitudes and hosts an approximately 3 m thick sub-horizontal ice body. According to Bertozzi et al. (2019), the LIC’s CAD can be classified as dynamically active. In the middle of the cave, a noticeable chimney drains water dripping which, during extreme precipitation events, can bring large amounts of water into the cave producing intense melting and the formation of deep bédière (glacial stream), as observed in November 2014 (Colucci et al., 2016). Inside the cave, air inflows even through the debris resulting in the formation of a dome-shaped cavity underneath the ice body that is observed to be getting wider with time. Progressive degradation of permafrost induced by recent climate change has been monitored over the years (Colucci and Guglielmin, 2019). Presently, permafrost has almost completely disappeared from the host rock, and density-driven air flows are the dominant control of the cave’s thermal regime. For these reasons, the ice body has been in a state of recession since 2014.

Methods

The proposed workflow is divided into three phases (Figure 2) and aims to be replicable in any underground cryosphere environment with only minor constraints. The first phase consists of photographs and GPR data acquisition, while the second one involves SfM-MVS and GPR data processing. The third phase consists of a comparison of the topography and volume variations between different entities preceded by some point cloud processing. This includes the coupling between SfM-MVS and GPR data outputs.OPEN IN VIEWER

The 12 SfM acquisitions performed within the VIC were acquired on a monthly basis from July to October in the summer seasons of both 2017 and 2018 (8 surveys) while they were made only at the beginning and at the end of the ablation season in 2019 and 2020 (4 surveys). Once the methodology proved to be effective in the more straightforward scenario of the VIC, it was also tested in the LIC, with acquisitions in July and November 2020 (2 surveys). In addition to the absence of light from outside, the LIC is characterized by a much more complex and irregular ice deposit that challenges the SfM-MVS capabilities.

All the photographs were taken with a full-frame (FF) sensor (36 x 24 mm) DSLR camera (Nikon D800 E) with 36.2 megapixels of resolution, equipped with 35 mm fixed lens (63.5° field of view on a FF sensor). All the files were recorded in a lossless Raw format (.nef) to take advantage of increased dynamic range during post-processing, given the complex lighting and shadow conditions that arise when using artificial lighting in caves. Except for a few frames, all photos were taken with a fairly closed aperture (>f/5.6), high ISO and shutter speeds of at least equal to 1/160 s to prevent possible micro-blurring. To be able to shoot in dark environments with these settings, it was necessary to use artificial lighting. This was successfully obtained with both a photographic flashlight (e.g. Nikon SpeedLight SB-700) and a 100W battery-powered flooding light that proved to be the better alternative, due to its wider and more uniform light diffusion capability. The use of these artificial light sources combined with the full-frame sensor low light performance allowed all photos to be acquired by hand, significantly reducing the acquisition time.

The ground control points (GCPs) referencing surveys were made once per cave using a lightweight instrument already widespread in cave surveys, namely, the Disto-X (Leica Camera AG). It is an electronic surveying device for cavers consisting of a distance metre with a built-in replacement board that extends the functionality of the instrument by a three-axis electronic compass/clinometer and a Bluetooth connection. The three-axis compass allows measurements of directions with arbitrary orientations of the device without degradation of the precision (0.002 m with a 0.5–100 m range). Each georeferenced point cloud was used as a source of coordinates. In practice, the distance between some easily distinguishable features in the rock walls of the georeferenced point clouds were then used as scale bars to optimize and scale the other point clouds, based on the recognition of the same physical features.

The GPR acquisitions were made with a ProEx GPR equipment (Malå Geoscience) connected with 250 MHz and 500 MHz shielded antennas as a function of the different objectives of the performed surveys and of the expected ice thicknesses (e.g. Colucci et al., 2016). The instrument was triggered by an electro-mechanical odometer allowing a constant trace interval in the range from 0.02 up to 0.15 m. Inside the VIC (Figure 1(c)), 11 profiles were acquired, while inside the LIC (Figure 1(d)), 7 profiles were acquired.

The entire SfM-MVS processing chain was carried out using Agisoft Metashape suite, following the state-of-the-art steps that make up the increasingly popular applications of this methodology in the geosciences (e.g. Carrivick, 2016; Piermattei et al., 2016; Smith et al., 2016; Verma and Bourke, 2019), with notable methodological differences only in the georeferencing phase (also depending by the scale of interest). After feature extraction from each frame and keypoint correspondence evaluation pair by pair, a sparse point cloud reproducing the acquired scene was generated. Given our datasets consisting of 100–400 photos maximum, the standard Metashape’s 40,000 key points and 4000 tie point limit for every frame provided a qualitatively adequate output. According to Alfio et al. (2020) as well as to our preliminary tests, there is no further increase in quality and accuracy from directly using a compression-less format such as TIFF (Tagged Image File Format) instead of high-quality jpgs. All the sparse point clouds generated after this first step of the processing lack an absolute scale, and therefore need referencing with the introduction of external spatial information coming from a surveyed network of targets (GCPs). The GCPs can also serve as check-points for error estimations of the methodology (e.g. Sanz-Ablanedo et al., 2018). The sparse point clouds were therefore optimized by re-executing the SfM and then implemented through MVS algorithms (e.g. Furukawa et al., 2010). The obtained output were dense point clouds composed of hundreds of millions of points. This process was streamlined in very complex scenes (e.g. LIC) by working on different sections of the cave independently, and then aligning and merging them together in a later step. The dense point clouds were at this point ready for further analyses that would lead to the actual multi-temporal comparison. It may be useful, even if this is not among the aims of this workflow, to produce a polygonal mesh from the available clouds that can be lately texturized, thus providing a photorealistic 3D model of the caves. Another derivable product, especially useful for surface morphologies analysis, is the DTM of the permanent ice deposits.

GPR datasets were processed using a standard workflow, which included drift removal (zero-time correction), geometrical spreading correction, bandpass filtering, background removal, 2D migration and depth conversion. The primary purpose of the processing was to increase the signal/noise ratio to improve imaging thus making the data more easily interpretable and informative. The EM speed considered in the entire processing flow was 17 cm ns⁻¹, that is, the typical value for pure ice, while for the cavity within the ice in the LIC, we set a velocity equal to 30 cm ns⁻¹ (air velocity). The actual speed of EM waves through the ice was verified by analysing the diffraction hyperbolas, finding non-significant variations from the reference value previously chosen. Horizon interpretation was done manually exploiting the software Petrel (Schlumberger) originally developed for reflection seismic interpretation. Dedicated attributes analyses were used to better constraint the horizon picking (Zhao et al., 2016).

The third and last phase, including the actual comparison and data integration, was entirely made with an open-source point cloud processing software (CloudCompare®). The initial steps were focused on noise reduction through a statistical outlier removal (SOR) filter that calculates for every single point an average distance from the near neighbours (NN) number of our choice and discards the points at a greater distance. Each cloud was then subsampled with a minimum space between two points of 0.002 m to streamline subsequent calculations without significant loss of spatial resolution. As the multi-temporal point clouds were very different due to ice melting inside the caves, they needed to be finely aligned through point picking in the stable portions of the caves (rock walls). Once those 10–15 notable feature points were used for the alignment, quality control was executed performing a cloud-to-cloud distance (C2C; Girardeau-Montaut et al., 2005) involving the stable rock walls only. The same C2C approach could be applied to the entire scene for change detection and measurements as it is fast and can generate X-Y-Z direction scalar fields simultaneously. However, as evidenced by Lague et al. (2013), the C2C measured distances are too sensitive to cloud roughness, outliers and point spacing (Figures S2, S3). The alternative used was the multiscale model to model cloud comparison (M3C2; Lague et al., 2013), an algorithm that works fully in 3D, does not require meshing or gridding of the point cloud and does not compute a difference when there are missing data. The M3C2 estimates a surface normal vector (N) for each core-point of the reference cloud and then compares the two clouds along N; using the local point roughness, size and co-registration error to estimate a local confidence interval (Lague et al., 2013). The computation of geometric features based on the spatial information of 3D points in a local neighbourhood (e.g. Weinmann et al., 2013; Hackel et al., 2016) can be useful for automatically selecting and isolating a certain area of the cave in which to compute the M3C2 comparison. If the permanent ice deposit is sub-horizontal, as in the VIC and LIC, it is possible to use the verticality attribute (Weinmann et al., 2013) to filter the entire point cloud while retaining only a certain range of values.

Two different approaches were tested out to compute the volume differences (Table S1). The first one used cloud gridding (0.002 grid step) to calculate volumes in 2.5D (i.e. difference of DTMs; DoDs), interpolating average relative heights for missing data. The computed global volume difference equals the sum of all the volumes of the elementary parallelepipeds (grid step*grid step*difference of height). The second approach involved another software, 3DFlow Zephyr, with volume with projection algorithms. This method meshes every point cloud using Delaunay triangulation and then calculates the volume between the entity (i.e. ceiling) and a chosen plane under it (i.e. ground). The resulting values of this second approach were used just as feedback for the DoDs results, which seemed to be more robust to small scale variations between point clouds.

In the final phase of the workflow, the interpreted GPR horizons were interpolated to generate a 2.5D grid (Kriging gridding method) of the bedrock depth or the ice thickness throughout the scene. The lack of absolute georeferencing was solved by using visual references in the point clouds set as marker points during GPR data collection. The overall error detected by the positioning was almost entirely due to the variation of the scene from the GPR (2012) and SfM (2017–2020) acquisitions. The total ice volumes were lastly computed in 2.5D with a DoDs approach.

Results

The entire set of dense point clouds generated through SfM-MVS (Figure 3; Figures 6(a), (b); Table S2) provided accurate and reliable data for multi-temporal comparisons. The only scene that was not successfully reconstructed was the VIC in July 2020 (Figure 3, VIC11), where the presence of a surface pond made the methodology ineffective. As expected, the first surveys carried out inside the VIC generated dense point clouds with more empty areas due to the lower number of camera stations throughout the scene. The total reconstructed area (including rock walls) was 1645 m² for the VIC (excluding the lateral chimney) and 961 m² for the LIC (including part of the main chimney). It should be highlighted that the central areas of VIC and LIC, where the photographic coverage is greater, were more comprehensively resolved and less noisy than the areas at the margins of the caves. These external areas were included in the acquisitions to get additional data but were not fundamental for the purposes of this work.OPEN IN VIEWER

As far as GPR data, picked horizons (i.e. ice bottom for both VIC and LIC and top of the air-filled cavity for LIC) were exported in relative X, Y, Z coordinates (Figure S3) and combined with the dense point clouds generated through SfM-MVS. In particular, we exported such a data for the portions of all profiles in which the horizons were interpreted in every 8 or 15 cm (along track) depending on the original spatial sampling interval of profiles.

With regard to the VIC, it was possible to reconstruct both seasonal and multi-year variations of the permanent ice deposit (Table 1). The presence of snow within the cavity did not follow any particular trend while an increasing volume of debris, mainly derived from rock wall collapse, was observed over the years. The use of M3C2 algorithms made it possible to highlight the distribution and the intensity of the surface variations within the VIC. The largest observed values (slightly > 1.20 m) were all due to snow entering the cave from the outlet ice patch and the lateral chimney. A constant area of 175 m², fully including the permanent ice deposit, was used for all the seasonal volume calculations (Table 1).

OPEN IN VIEWER

Table 1.

DoDs computed volume differences inside the VIC between different seasonal surveys. The difference between ice volume loss and total volume loss is entirely due to snow accumulation/melting inside the cave. The area considered for the estimates is constant and equal to 175 m².

Year	Period	Ice volume loss [m3]	Total volume loss [m3]
2017	July–August	−64.05	−100.36
2017	August–September	−5.67	−5.67
2017	September–October	−1.68	−1.68
2017	July–October	−71.41	−107.71
2018	July–August	−14.50	−47.85
2018	August–September	−13.41	−41.63
2018	September–October	−1.95	−8.68
2018	July–October	−29.85	−98.17
2019	July–October	−24.80	−64.96

To separate volume differences related only to ice melting phenomena from the total volume difference, several corrections to consider snow and debris variations were applied to the multi-year volume measurements (Table 2). The reconstructed 2012 surface was approximated to a plane, considering the almost flat geometry of the deposit at the time, with almost no observed debris on the ice surface. All the annual volume differences were carried out on the last acquired survey of each ablation season (Figure 4; VIC4, VIC8, VIC10 and VIC12). A major benefit of the M3C2 algorithm for this application was the possibility of displaying both positive and negative mass variations from season to season (Table 2).

OPEN IN VIEWER

Table 2.

Yearly volume differences inside the VIC from SfM-MVS surveys (2017–2020 period) and estimated volume loss from the 2012 GPR interpreted ice depth. The surface area considered is not constant due to the presence of snow or debris and rocks in the scene. The total volume difference includes also the snow melting and it is therefore influenced by the amount of snow entering inside the VIC.

Period	Ice volume loss [m3]	Total volume difference [m3]
2012–2017	−45.38	−45.38
2017–2018	−10.95	+36.55
2018–2019	−25.37	−27.21
2019–2020	−52.30	−65.33

OPEN IN VIEWER

Figure 4.

(a) Planar view of the RGB dense point clouds at the end of each ablation season (2017–2020) for the VIC; (b) seasonal M3C2 differences: 2017–2018, 2018–2019 and 2019–2020; (c) for each of the M3C2 related histograms showing the distribution of the computed values in 256 classes from −1.0 m to 1.0 m. The point clouds are always compared starting from the oldest to the most recent one. Negative and positive M3C2 values indicate a mass loss and a mass gain in the considered period, respectively. Please notice that mass losses are not always associated with ice variations but can be related also to the input/output of snow or debris.

Computing the verticality on point clouds enhanced the visualization of geomorphological features and water drainage dynamics involved in the VIC’s deposit during the ablation season (Figure 5). The bédière carried the water to the same area between the ice and the rock that probably served as a drain for the surficial lake that was observed above VIC’s surface. The presence of debris was much more significant towards the end of the season when there was generally only old ice inside the cave. The melting of all the new ice formed in each winter is always observed before the end of the ablation season.OPEN IN VIEWER

Despite the artificial lighting of the scene, the entire LIC was also accurately reconstructed (Figure 6 and Figure 7). The much shorter time series compared to the VIC did not allow a continuous long-term observation; nevertheless, data were retrieved from 2012 GPR surveys. In contrast to the VIC, where volume increases were also observed due to the presence of snow, only volume losses during the surveyed interval of time were observed. On the surface of the ice deposit, the highest melting occurred below the main chimney and above the inner cavity. In the latter, the rapid growth of a hole was observed (maximum diameter ≈ 2.5 m; Figure 6). The volumetric losses between July and November 2020 were significant, averaging about 10–15 cm in thickness with peaks of 30 cm. There were some areas further into the cave where the ablation was lower. Here, it was not clear from our reconstruction if the point cloud was reliable enough in terms of accuracy (Figure 6(c)). Given the complicated air circulation already observed within the LIC (Bertozzi et al., 2019), it was also interesting to focus on ice melting phenomena acting in the inner cavity (Figure 7). In this area, melt rates seemed more heterogeneous. From the M3C2 of the ceiling only, the differences in thicknesses between July and November ranged between a few centimetres at the entrance to more than 15 cm in some of the internal portions of the cave (Figure 7(a)).OPEN IN VIEWER

OPEN IN VIEWER

Figure 7.

M3C2 (July–November 2020) of inner air-filled cave ceiling inside the LIC (a) and a 3D view from July 2020 textured mesh (b). Holes at markers 1 and 3 are artifacts due to errors in the SfM reconstruction. The central hole is real and appeared in the interval between the two considered surveys.

The analysis and interpretation of the GPR profiles inside the VIC revealed a maximum ice thickness of about 8.5 m (Figure 8(b)). On all profiles, the presence of debris over the surface with centimetric to decimetric diameter could be seen, as confirmed by the morphology observed in the subsequent years. Below this debris-rich ice layer, the ice appeared cleaner with the rare diffraction hyperbolas fitting perfectly with the typical values of pure ice (17 cm ns⁻¹). The total ice volume present at the time of the GPR surveys was estimated as 987 m³. The LIC ice deposit was thinner, reaching a maximum of 4.5 m (Figure 8(d)). The bedrock’s shape was easily reconstructed as an irregular funnel with a slope of about 25°. Near the LIC entrance, the presence of debris in the ice was more significant, although in the 2020 surveys the presence of decimetric or centimetric debris increased. The 2012 reconstructed inner cavity ice ceiling (Figure 9) made it possible also to calculate the air volume changes inside the inner cavity. The total air volume in 2012 was 22.43 m³, while in July 2020, it was 48.01 m³ reaching 54.19 m³ in November 2020. The bedrock surface differences observed between the 2012 GPR and the 2020 SfM (Figure 9) could be partly due to geophysical errors (e.g. related to the velocity used for depth conversion) and/or to the profile positioning; however, they can be also related to the fact that debris within the ice fell from the melting vault zone into the area below, thus producing a local elevation of the ice cave bottom (Figure 9, difference between orange and light blue levels). The 2012–2020 volume difference came from the estimated relative height of the 2012 ice deposit. Overall, the LIC bedrock was less dipping and more regular compared to the one in the VIC (Figure 8). The total ice volume present at the time of the GPR surveys in the LIC was estimated to be 365 m³.OPEN IN VIEWER

OPEN IN VIEWER

Figure 9.

Profile P6 within the LIC (see Figure 8(d) for the position) with interpreted bedrock and inner cavity ice ceiling from GPR surveys (2012) are compared with the bedrock and the ice surface from SfM-MVS survey (July 2020). See text for details.

During the 8 years (2012–2020) between GPR and SfM-MVS last surveys, the VIC lost 13.6% of its estimated ice volume, which equates to 134.00 m³. In the same period, the LIC lost 49.4% of its ice volume, which equates to 180.10 m³ of ice.

Discussion

Terrestrial SfM-MVS in underground environments has great potential, but its reliability is highly dependent on the quality of the photos, their overlapping and the geomorphological setting of the surveyed area. Another critical point is the presence of ice, whose surface is often lacking in easily distinguishable features and rich in reflections that are inconvenient for keypoint correspondence algorithms, especially when artificial lighting is used. Therefore, the errors derived from our measurements should only serve as a possible reference example. The ice in VIC and LIC is rich in heterogeneously sized debris and is therefore not as uniform as in other caves. This characteristic combined with the proximity to the target and the abundance of overlapping photos made it possible to achieve a qualitatively remarkable result. Although not recommended, we found no negative feedback in the combined use of vertical and horizontal photos as a source of points for SfM-MVS.

Given the scarcity of studies on the subject and the fact that these focus on discrete locations measured by stake reading or ultrasonic range sounding, it was difficult to compare our results with broader findings on ice losses rates in karstic cavities around the world, although they all indicate ice reserves in caves are diminishing (e.g. Rachlewicz and Szczucinski, 2004; Obleitner and Spötl, 2011; Perşoiu and Pazdur, 2011; Schöner et al., 2011). Continuous measurements could still be a good solution for seasonal monitoring where a network of sensors is available, whereas SfM-MVS is much better suited to multi-annual monitoring as it takes into consideration the entire ice deposit with a significant increase in the overall accuracy of the volumes involved.

To our knowledge, no positive multi-year mass balances were observed in other locations in recent decades inside ice caves located in temperate areas, but the rate of ice loss could be expected to be spatially very variable. The interplay between different CAD, ice accumulation mechanisms and external climate is rather complex and site-dependent. In an ice cave in Ciemniaku (Tatra Mountains, Poland), significant changes comparable to those we measured in the VIC and LIC were observed (−15–30 cm y⁻¹ of ice, for example,. Rachlewicz and Szczucinski, 2004). Other larger and much developed ice caves have not experienced this magnitude of loss, such as Eisenriensvelt cave in Austria (−3.5 cm y⁻¹, Obleitner and Spötl, 2011) and Scărioara ice cave in Romania (−1.5 cm y⁻¹, Perşoiu and Pazdur, 2011). According to Colucci and Guglielmin (2019), the LIC was similarly stable (less than −5 cm y⁻¹), until permafrost degradation, combined with a series of severe precipitation events, occurred in 2014 and abruptly accelerated the ice loss rate (up to −10/−15 cm y⁻¹).

The SfM-MVS errors (E_SfM) were equal to 0.01 ± 0.01 m (Table S2) in the targeted areas of reconstructed dense point clouds. This was achieved by excluding a few GCPs from the SfM-MVS processing and using them as check-points. In addition to these, two metric scale bars were used in each cave. The difference between the distance of these check-points in the reconstructed model and the Disto-X topographic survey represents the error due to SfM-MVS reconstruction. This could be partially attributed to the hand-held Disto-X measurements. A second and less significant error comes from the alignment of different point clouds through manual point picking in CloudCompare® (E_reg). This equates to 0.005 ± 0.005 m and has been measured performing C2C distance algorithms on the stable portions of the scenes after their alignment. The total error (E_tot) of any point forming the aligned point clouds can be estimated, according to Collins et al. (2009) and Šupinský et al. (2019), with a deterministic error analysis where

E t o t = E S f m 2 + E r e g 2

E_tot was determined to be 0.01 ± 0.01 m. As the survey benchmark (GCPs) are always the same, there is no need to consider GNSS measurements or other absolute referencing errors (Collins et al., 2009). The above-mentioned errors are related to the best-coverage areas of the cave, where the overlap between the photos is high (i.e. > 9 cameras). This error estimation, based on independently measured check-points, provides only a global value and masks the spatial variability that results from photogrammetric and georeferencing effects. Unfortunately, state-of-the-art photogrammetric and georeferencing uncertainty and confidence-bounded 3D topographic change detection used for UAVs surveying (James et al., 2017) are not possible without absolute coordinates.

As far as GPR data, the most relevant error can be related to inaccurate estimations of the EM velocity, in turn involved in the time-to-depth conversion. The velocity estimation strategies are limited for common offset datasets (e.g. Forte and Pipan, 2017) which, on the other hand, are the only possible solution within ice caves due to the obvious logistical constraints making impossible multifold techniques. However, for cold ice conditions, the error is negligible, while when some free water is present, it can be in the order of 5% of the estimated depth, which is overestimated since the presence of unfrozen water decreases the electromagnetic velocity (Bradford et al., 2009). In addition to such uncertainty, there are errors related to the intrinsic resolution of GPR system mainly depending by the frequency of the antenna used and, to less extent, to the depth of the target and to the electromagnetic velocity of the surveyed materials. The major limitations of the SfM-MVS in underground environment are due to the presence of narrow passages where each acquired frame cannot provide sufficient matching with the others for a reliable 3D reconstruction. In those places, other tools such as the TLS would also be ineffective. Another situation where the SfM struggles is in the cavities entrance areas, where external light and internal darkness make the dynamic range of photo exposures too wide to recover the details during the post-processing of the photos. This aspect and a higher detail in the textured models could be achieved with a better artificial lighting set-up, ensuring uniform illumination across the scene or shooting every photo from a tripod. In both cases, the acquisition time would increase enormously.

Increasing both the global reliability and accuracy of the surveys can be achieved using a more precise data referencing system, both for SfM-MVS (e.g. TS) and GPR surveys. Besides, performing all acquisitions at the same time would allow much more precise GPR profile positioning and would also make error estimations of this SfM-GPR coupling procedure possible.

The proposed combined surveying technique offers a versatile, fast-to-acquire and relatively inexpensive approach for analysing underground cryospheric environments. The excellent resolution obtained makes this methodology ideal for investigating any variable geometry in remote underground environments, where light and space-saving instruments are mandatory. The methodology itself can be applied to other branches of glaciology where the combined SfM-MVS and GPR multi-temporal surveys would be a convenient solution to simultaneously study the surface and internal evolution of glacial bodies with ease.

Conclusions

The joint application of multi-temporal SfM and GPR inside ice caves led to the following methodological and applicative conclusions. The dense point clouds and the polygonal meshes reconstructed through SfM-MVS allowed the reproduction of the inner cave with a sub-centimetre resolution. The point cloud comparison made it possible to calculate the distribution and entity of topographic changes within the caves, thus obtaining information on mass balance, as well as on snow, debris and geomorphological settings in each scene. The GPR allowed volume estimations of the ice deposits with high precision, especially in cold ice conditions, proving to be an excellent and versatile methodology for investigating the underground cryosphere. The only limitations are related to the logistics of acquisition in such narrow spaces, often not allowing free movement within the cave, as well as problems related to the complex absolute georeferencing.

The proposed methodology constitutes a tool to a general improved quantification of ice volume reduction occurring in alpine caves under global warming scenarios. Focussing more specifically on our results for volume changes over the 2012–2020 period, it is evident how inside the Leupa ice cave where water dripping is the primary source of ice accretion, ice losses are much more significant. The Vasto ice cave, where firnification processes play a fundamental role in the formation of new ice, is apparently more resilient to multi-year ice losses and should be better investigated in future glaciological and geomorphological studies.

In synthesis, we conclude that SfM-MVS combined with GPR surveys presents a simple and robust procedure to reconstruct the geometry of all structures above and below ice deposits inside caves allowing the seasonal and multi-annual monitoring of surface and subsurface changes.

Supplemental Material