Assessment of Morphology Changes of the End Moraine of the Werenskiold Glacier (SW Spitsbergen) Using Active and Passive Remote Sensing Techniques

: Wedel Jarlsberg Land in Svalbard is a region with a varied periglacial landscape. In the mountains and in the valleys, the climate is polar with permafrost. During the summer, the near-surface ground layer thaws. The Werenskiold Glacier, together with its end moraine, are located in the central part of this area. The rate of morphological changes observed within the moraine varies in time and space, and depends on the environmental conditions. This study investigates four periods of archival aerial photogrammetry measurements (1936, 1960, 1990, and 2011) performed for the end moraine of the glacier. The long-term analysis was also based on a direct GNSS RTK survey from 2015. Over a period of almost 80 years, more than 14 million m 3 of rock and ice material disappeared from the end moraine of the glacier (an average of approximately 200 thousand m 3 /year). Analyses of the dynamic surface changes over one year (2015) were performed with the use of synthetic aperture radar interferometry (InSAR). The time interval between images was in this case 12 days and covered (simultaneously in each scene) the entire investigated area. In this case, the analysis demonstrated that over a period of only 4 months, the moraine lost 200 thousand m 3 of material (approximately two thousand m 3 /day), which is equivalent to the entire annual mass loss of the moraine.

Global warming significantly influences the transformations of ice-covered areas in Spitsbergen, which is the largest island of the Svalbard archipelago. In this case, the situation affects not only glaciers, which melt together with an increase in the average annual temperature in the region [25,26], but also glacial moraines [27,28].
The methods employed to observe these phenomena are primarily based on geodetic measurements. Frequently encountered harsh terrain and/or weather conditions, as well as progress observed in the field of measurement instruments and techniques, are the reasons for the increasing popularity of the so-called remote sensing techniques. They can be classified into two groups: passive and active methods. The first group includes photogrammetry, which is based on the development and analysis of numerical terrain models from photographs. These can be taken from various levels (on-ground, low, and high altitude aerial or satellite images) and with the use of various types of cameras. The type of a photogrammetry mission is adjusted to the observed phenomenon, its scale, and the required/expected accuracies. Its significant disadvantage, however, lies in the need to ensure proper lighting conditions (data acquisition is possible only during the day) and weather conditions, which-if imperfect-can negatively influence or even completely compromise a measurement session [29][30][31][32].
The second group of methods has been dominated by SAR interferometry. SAR gives a high accuracy of displacement determination, but it refers to a larger averaged area. Photogrammetry allows for a satisfactory DEM accuracy at the level of 1-2 m, but the time intervals between the flights determine the generalization of the registered deformation process. This active system makes use of the so-called radar images (scenes), and subsequently, based on the calculated signal phase shift for image sequences, it allows the identification of the values of ground surface deformations/displacements; in particular, of vertical displacements. An image of phase differences thus obtained is referred to as an interferogram. Depending on the defined research problem, the characteristics of the research area and the expected accuracy (also determined by the estimated value of change in a unit of time), different approaches of classic radar interferometry can be used, such as the PSInSAR, CRInSAR, or SBAS methods [33][34][35][36][37][38][39][40][41][42][43][44][45]. The advantages of this system include the possibility to use both commercial and free data obtained remotely at relatively dense time intervals, which means that no need exists to interfere directly with the natural environment.
The main goal of this study is to provide a description and an interpretation of changes observed on the surface of the end moraine of the Werenskiold Glacier (Figure 1), which is located in the southwest part of the Spitsbergen [27]. The time span for the analyzed morphological changes was divided into two periods determined by the measurement technique used as the basis for this study. The first, long-term period covers photogrammetry measurements first performed in 1936 and continued until 2011. These were supplemented with the direct survey from 2015. As a result, it was possible to determine changes in the shape of the moraine surface (as well as the changes of the volume of the moraine itself) over a period of 80 years.

Materials and Methods
The long-term changes were analyzed on the basis of numerical terrain models developed from aerial photographs taken in 1936,1960,1990, and 2011 for the Polar Norsk Institute in Oslo (cf. Figure S2). These photographs are the only datasets available free of charge for the analyzed area. Low frequency of the photogrammetry missions in the area In turn, the InSAR technology was used to demonstrate the values of changes over 1 year (and more specifically-for the summer period). It showed (as will be further discussed below) that these changes occur at a fast pace and are impossible to be traced using traditional survey methods (e.g., GNSS). Importantly, photogrammetry is prone to external limitations (such as atmospheric precipitation, clouds, wind), and a GNSS survey is time-dependent: it is prolonged and requires the involvement of a significant number of people. For the above reason, it was impossible to select only one universal method to describe the problem addressed by the authors of this article.

Materials and Methods
The long-term changes were analyzed on the basis of numerical terrain models developed from aerial photographs taken in 1936,1960,1990, and 2011 for the Polar Norsk Institute in Oslo (cf. Figure S2). These photographs are the only datasets available free of charge for the analyzed area. Low frequency of the photogrammetry missions in the area has been mainly the result of difficult weather conditions. Depending on the date of the photographs, they are either black and white or in color. The photographs are orthogonal, taken with a camera installed under the plane fuselage. The photographs from 1936 are an exception in that they are oblique and taken from the board of the plane (Figure 2a).  (Table S1).
In the first step, the input data (Table 1) in the form of images were used to develop the Digital Surface Model (DSM) for each analyzed year. The Digital Terrain Model (DTM) was based on data in the form of a point cloud obtained by processing photogrammetric data in the Agisoft computing environment-with the use of Agisoft PhotoScan Professional tools. The height field surface algorithm with the closest neighbor interpolation served to prepare a Triangulated Irregular Network (TIN) from the point cloud. All Digital Terrain Models (DTM) were fitted to the Universal Transverse Mercator (UTM) system, in the 33X zone, on the basis of 7 characteristic points identified in each photograph in each of the measurement periods (the location of these points is presented in Figure S3). The accuracy of the models was within 1-2 m in the case of horizontal components and within 2-3 m for the vertical component.  The analyses were also based on the DTM developed from direct in situ surveys of the surface of the end moraine performed in 2015 ( Figure 3). The surveys followed the GNSS method and used the reference stations existing in Spitsbergen: NYA1 in Ny-Alesund and ASTRO in the vicinity of Polish Polar Station in Hornsund. It is worth mentioning here that the length of the analyzed area of the end moraine is approximately 3 km, the distance between the measurement points and the closest reference point in Hornsund is approximately 12.1 km, and the distance to the farthest point is 14 km. However, the distance to the second reference point (Ny-Alesund) is above 200 km. Corrections from the reference points were obtained after the measurement, at intervals of 5" and 30" from the closest station, and of 30" from the farthest station. The data was postprocessed and adjusted to the two reference points. The GNSS measurement has the greatest accuracywith allowance for the corrections from the reference points, and the location accuracy for each of the measurement points was at a level of below ten centimeters (for the horizontal components) and above ten centimeters (for the vertical components). The lowest accuracy obtained for a single point was characterized by horizontal accuracy of ±8 cm and vertical accuracy of ±11 cm as well (this level of accuracy was achieved in the post-processing mode). The highest achieved accuracy of points measured in the GNSS RTK method was ±1.5 cm for the horizontal components and ±2.5 for the vertical component. The survey was performed by 4 people and lasted for 6 weeks, providing a total of 20,721 points. pairs in common images, using Atmospheric Phase Screening (APS) for each scene [57]. All of the InSAR results spatially relate to the first image of August 02, 2015. The implementation of multiple time series reduces the atmospheric effects owing to the assumed temporarily non-correlated tropospheric effects [58,59].
Based on the interferograms, the authors estimated the displacement time series by following the Small Baseline Subset (SBAS) method [39]. The phases were reversed with the use of the cost function based on the L1 standard, in relation to the unwrapping errors [60]. The results for the analysis of the InSAR time series for the data set correspond to single-dimensional (1D) displacement speeds along the LOS (Line Of Sight; details of the analysis are provided in Table S2). On the map, the obtained pixel corresponds to an area of 40 × 40 m.  The short-term changes were analyzed with the use of synthetic aperture radar interferometry (InSAR). The scenes used in the study were acquired from the Sentinel-1A satellite, Single Look Complex (SLC), Path 14, Ascending, Polarization HH, Frequency band C (λ: 5.55 cm), LOS (orientation/incidence angle: 69.5 • /37.3 • ). It allows the terrain surface changes to be traced at practically any time intervals. The only limitation is due to time periods in which the radar images were taken and the availability of the source materials in the European Space Agency (ESA) service or in its mirror services [48][49][50].
The InSAR calculations were performed according to an algorithm, with the use of the GMTSAR algorithm [51]. We co-registered and multi-looked single-look complex (SLC) images using a range/azimuth multi-looking factor of 2.3 × 13.9 m pixel resolution, providing a ground resolution of approximately 40 × 40 m. Simplified locations for the radar scenes are provided in the supplementary material (cf. Figure S2).
Due to a significant diversity of the investigated processes and to the substantial change rate in the end moraine, the interferograms were generated at a maximum base timeline of 24 days, which allowed data coherence and a possibility to determine phase differences [39]. The interferences occur when the displacement speed exceeds one fourth of the wavelength during the time interval. In the case of the prepared interferograms, this is 1.39 cm, i.e., a period of 6-24 days. Spatial base lines were limited to 150 m (details to be found in Figure S1). This study was based on the scenes from 2015. The details  . Table S2). However, the presentation includes only seasonal displacements (June-November). The data series was ended in November, as the winter snowfall causes the images to be decorrelated [52]. The noise level was reduced in all interferograms by employing a spatially adaptive, coherencedependent Goldstein filter [53,54]. Strongly decorrelated interferograms were removed and not involved in the calculations of the deformations. The impact of the atmospheric layers was mitigated by fitting the linear relationship between the residual phase and the topography [55], using the Digital Elevation Model (DEM) obtained from the GNSS surveys of 2015. Pixels with noise were removed by employing a coherence filter (coherence above 0.3-0.48 for 48% of the interferograms). The interferograms were unwrapped with the use of the Statistical-Cost, Network-Flow Algorithm for Phase Unwrapping (SNAPHU) [56]. The interferograms were corrected by averaging all of the pairs in common images, using Atmospheric Phase Screening (APS) for each scene [57]. All of the InSAR results spatially relate to the first image of August 02, 2015. The implementation of multiple time series reduces the atmospheric effects owing to the assumed temporarily non-correlated tropospheric effects [58,59].
Based on the interferograms, the authors estimated the displacement time series by following the Small Baseline Subset (SBAS) method [39]. The phases were reversed with the use of the cost function based on the L1 standard, in relation to the unwrapping errors [60]. The results for the analysis of the InSAR time series for the data set correspond to singledimensional (1D) displacement speeds along the LOS (Line Of Sight; details of the analysis are provided in Table S2). On the map, the obtained pixel corresponds to an area of 40 × 40 m.

Results
The quality of the developed terrain models varies. The oblique photographs from 1936 allowed the model of the terrain to be precise with respect to altitudes-the fitted photographs advantageously result in a dense point cloud. The effect is satisfactory owing to the use of seven photographs taken from different locations. The pixel resolution, however, is on the order of more than ten meters, and as a result affects the quality of the model. The locations of the photographs are available at [46]. The 1960 model is acceptable, as it shows the differences of both the altitudes and the terrain forms. The point cloud from this model is sufficiently diversified-it is based only on three photographs showing the end moraine of the glacier. The black and white photographs have a very good light exposure and a resolution in the order of more than ten meters, which increases the value of the model. The 1990 model has the worst quality. The point cloud is blurred and "flat", and does not allow the values of the surface changes to be fully identified. The model is based on three photographs. The images have a very small depth of field-the photographs are in color, but the colors do not reflect the altitude differences. The pixel resolution is in the order of more than ten meters. The 2011 model is probably the worst, despite the quality of the input data used (Figure 4). The reason for this lies in heavy clouds on the day of the measurement, which led to the development of pits and spikes in the model.
The results of the GNSS direct survey are the most accurate and therefore allow the most precise DTM model. This precision results from the possibility to allow for corrections in identifying positions in 3D space provided by the reference points (GNSS permanent stations). The obtained location precision for each of the measurement points is in the order of several centimeters (for the horizontal components) and above ten centimeters (for the vertical component).
The SAR-based analysis of the surface shape of the end moraine indicates very dynamic changes. The surface of the moraine is undergoing constant modifications, with unexpectedly occurring sink holes and scarps, as well as uplifts and banks in their vicinity. The accuracy of the calculated displacements in the quasi-vertical (LOS) plane is ±3 cm (the analysis did not involve changes in the horizontal plane), and the displacements observed on the surface of the end moraine are greater than the measurement errors, suggesting that the employed method was adequate. cloud is blurred and "flat", and does not allow the values of the surface changes to be fully identified. The model is based on three photographs. The images have a very small depth of field-the photographs are in color, but the colors do not reflect the altitude differences. The pixel resolution is in the order of more than ten meters. The 2011 model is probably the worst, despite the quality of the input data used (Figure 4). The reason for this lies in heavy clouds on the day of the measurement, which led to the development of pits and spikes in the model. The results of the GNSS direct survey are the most accurate and therefore allow the most precise DTM model. This precision results from the possibility to allow for

Long-Term Analysis
The DTMs prepared on the basis of the photogrammetric data were subjected to a differential analysis for the successive measurement periods. The obtained results allow both the evaluation of the morphological changes, i.e., the changes of the geometry of the end moraine ( Figure 5) and the calculation of the volume differences for the entire structure ( Table 2).

Long-Term Analysis
The DTMs prepared on the basis of the photogrammetric data were subjected to a differential analysis for the successive measurement periods. The obtained results allow both the evaluation of the morphological changes, i.e., the changes of the geometry of the end moraine ( Figure 5) and the calculation of the volume differences for the entire structure (Table 2).  The greatest volume change occurred between 1960 and 1990. The outflow of the material from the end moraine was at 8.5 million m 3 , with an inflow of 4 million m 3 . In the subsequent period, a greater part of the moraine material was accumulated-more than 4.5 million m 3 . However, this result may not be reliable, and a significant measurement error is possible. The observed precision of the DTM allows a conclusion that the 1990 model has a limited accuracy. Its bias is indicated by the significant changes with respect to both the previous and the following models. The model is flat, which is well visible in the longitudinal cross-section of the overlapped models from individual measurement periods.
The longitudinal profile is drawn along the entire length of the analyzed object and across each of its surfaces ( Figure 6). The graph shows the greatest changes in the northern part of the moraine; it is also visible in the photogrammetric photos attached in the study. The cause is most likely the central moraine, which in the 1940s reached the end moraine. Then the central moraine was completely melted, retreating more than one kilometer along with the glacier, as seen in the photo from 2011. The central moraine is also made of an ice core, and as it melted, it released a large amount of water, which washed out the northern part of the end moraine. The process is visible on the longitudinal profile. In addition, the attached photo from 1936 clearly shows that in the northern part of the moraine there was the highest terrain point, which currently is in the middle part of the moraine.  The morphological changes of the moraine vary. The greatest changes, and this being even despite the limited accuracy of the spatial models, were observed in the northern part of the research area, where the terrain subsidence reaches~33 m-in its central part the subsidence is~6 m, and in the southern part it is~18 m. Over the period of 80 years, a rise visible in the northern part of the moraine was degraded. Additionally, in its northern part, the end moraine was approached by the medial moraine, which is now completely washed, as visible in the photogrammetric images ( Figure 2). The total change of the moraine volume between 1936 and 2015, as calculated by comparing the DTMs, reached almost 14 million m 3 . The washing of the moraine in its northern and southern parts may result from the path of the glacial river bed. Earlier, the bed extended over the northern end of the moraine, whereas presently no glacial river is found in the northern part of the moraine and the bed is completely dry. The glacial river currently extends in the southern part of the moraine and most probably washes it. The end moraine of the glacier is built of the ice core [27,61], and the flowing river has a significant impact on the thawing of the ice, resulting in the degradation of the moraine manifested in its changing geometrical shape.

Short-Term Analysis
The measurement results and the models constructed on their basis were used to prepare differential analyses. The scene of 2 August, 2015 is the reference image, and each successive image is the difference between the shape of the terrain surface and the reference scene. The deformation models were named in accordance with the image dates. The displacement values are symbolized with the pixel color, on a scale from dark blue in areas of most significant subsidences to dark brown in areas of the greatest rise. Stable areas are represented in white. The LOS displacement values are shown in millimeters (Figure 7).
The above allowed a conclusion that the most significant changes on the surface of the moraine occurred in August and September 2015 (Table 3), which may be related to the atmospheric conditions. The turn of September was a period of the most intensive precipitation over the three analyzed months, as measured in the Hornsund station [25,26], which is located at a distance of 12 km in a straight line from the research area. The least significant changes were identified in September, possibly due to stable weather conditions observed in the area over that period. The total change of the material volume in the moraine reaches 250 thousand m 3 . The weather variations were not recorded in the research area, and the 12 km distance in a polar region separated by a mountain range may be a source of error and lead to the misinterpretation of the obtained values.   The results of the longitudinal profile analysis allow a conclusion that the entire moraine is subject to a phenomenon of general subsidence, albeit at a varied rate depending on the part of the moraine. The speed of the changes may result from different reasons, such as exposure, thickness of the moraine material cover, the presence of rockfalls, and flows or permanent water reservoirs (e.g., lakes) and waterways in the surface of the moraine. The displacement diagrams for two selected points located within the profile (Figure 8) indicate that the changes in these points reach a value of above 10 cm.

Accuracy of the Results
The accuracy of the results depends largely on the employed measurement method. For research purposes, the accuracy analysis was simplified and limited only to the order of magnitude for the values obtained in the measurements. The main results of this research are the calculated vertical changes of the end moraine surface of the analyzed glacier. The altitude determined on the basis of the photogrammetric images from 1936, Additionally, an analysis was performed of the obtained results, i.e., the DTMs from satellite images, constructed along the longitudinal profile drawn between the southern and the northern part of the moraine (Figure 8).

Accuracy of the Results
The accuracy of the results depends largely on the employed measurement method. For research purposes, the accuracy analysis was simplified and limited only to the order of magnitude for the values obtained in the measurements. The main results of this research are the calculated vertical changes of the end moraine surface of the analyzed glacier. The altitude determined on the basis of the photogrammetric images from 1936, 1960, and 2011 remains accurate to 5 m and depends on the size of the pixel and the corresponding terrain unit. The photogrammetry measurement from 1990 is less precise-the altitude is calculated with an accuracy of approximately 10 m. The accuracy of the direct survey from 2015 is in the order of 10 cm and depends on the number of satellites observed during the measurement and on the corrections obtained from the reference stations. If the changes of the surface shape are considered as a difference between two models, then the accuracy of displacement calculations in one point (pixel) is a value of ±7 m in photogrammetrybased methods. For the entire, very diversified, research area (the end moraine), the error value of the calculated difference between the surface models reaches as much as ±300 thousand m 3 .
In the case of InSAR imaging, this error is in the order of ±30 mm for the calculations of pseudo-vertical (LOS) displacements in which the surface changes are considered as a difference between successive images. The on-surface location of points for which the displacements were calculated depends on their precise identification. The volume calculation error in this case is of ±500 m 3 .

Conclusions
The long-term measurements indicate a tendency and a direction of changes which have occurred on the surface of the end moraine of the Werenskiold Glacier over the last 85 years. The moraine is changing, significantly losing its volume and collapsing inwards without changes on the surface, on which it is located. The most significant changes are observed in its northern part, where the hill located northwards over the medial moraine has undergone an almost complete degradation. The second region affected by the changes is the southern part of the moraine, where large-surface sink holes occur, forming a very diversified landscape.
The short-term changes in the shape of the moraine surface result, among other things, from rockfalls, creep of the erratic material, and favorable ground conditions. This is particularly well visible in the northern and southern parts, where the steep moraine slopes accelerate the process of loose material falling towards the moraine foot. The morphological changes are also influenced by the melting of the ice deposited in the moraine; the so-called ice core of the moraine, which is clearly visible in the form of numerous outflows. The analysis of SAR images enabled the detection of locations in which this process is the most intensive. This would be impossible with the use of traditional measurement methods such as photogrammetry, which is expensive, or direct surveys, which are expensive and timeconsuming. Unfortunately, the pixel resolution in InSAR is significantly smaller than in the above traditional methods. During traditional GPS-based measurements over the period of 6 weeks, the moraine changed its volume in a manner both noticeable and measurable with the use of SAR technology. Its deformations reached 10 cm, as illustrated in Figure 8. Unfortunately, the assumed method of GNSS surveys and result processing did not allow a similar standard accuracy level for the GNSS RTK technique during the 2015 measurement session. Photogrammetric images could allow a precise identification of locations most prone to erosion, but the cost related to taking and processing the photographs is very high. In addition, atmospheric conditions frequently prevent photogrammetry flights.
The analysis of both the short-term and the long-term changes allows the identification of the degradation processes in the structure of the moraine. The most significant changes were observed in the northern and southern parts of the moraine, which is most probably related to the changing path of the glacial river bed. The northern part does not have a river, which completely dried in the erosion process. The changes which can be observed in the aerial photographs significantly influence the washing and melting of the ice core in the moraine, especially during the summer period. The InSAR analyses confirm both the locations of the changes and the areas of their greatest dynamics.
The long-term changes are identified with a limited accuracy which depends only on the quality, in this case on the resolution of the aerial photographs, or on the accuracy of direct surveys. The analysis of the short-term changes, illustrated with the use of LOS images, which precisely indicate the areas of deformations, allows the identification of the locations in which the shape of the end moraine surface is most prone to changes. The size of the pixel is not a limitation in this case, as the emphasis is on the location of the area potentially affected by erosion.
Being instantaneous measurement methods, aerial photographs and SAR images enable a quick identification of current surface shape, with the entire moraine represented at the same time. Direct surveys require a long time to be performed (between 5 and 6 weeks). They are also an expensive solution, as is photogrammetry. The SAR technology does not involve any costs (in the case of the Sentinel 1 satellite) and at present is available at 6-day intervals.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/rs13112134/s1. Table S1: Photogrammetry parameters, Table S2: InSAR processing methods and parameters, Figure S1: Baseline plot of interferograms generated and selected for dataset, Figure S2: Location of the pictures and measurements, Figure S3: Location of photogrammetry control points.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy restrictions.