Image-Based Monitoring and Predictions of Terrain Feature Variations Under Accelerated Warming Conditions in Kaffiøyra Area, Svalbard

1. Introduction

Global climate change has been a prominent concern since the mid-20th century, impacting regions variably. Polar areas, especially the Arctic, face more severe effects. Research highlights rapid Arctic warming, known as Arctic amplification, driven by sea ice loss, heat flux changes, and cloud cover variations [1]. Arctic warming has been faster four times than the globe since 1979, based on a study using several observational datasets to analyze near-surface air temperatures over the Arctic regions [2].

Research in the Arctic has historically been challenging due to its inaccessibility and extreme weather conditions. However, remote sensing technologies have revolutionized our ability to study this region. Remote sensing provides crucial data that enhances our understanding of climate change impacts in the Arctic. Various vital indicators extracted from remote sensing monitoring for the Arctic, including surface temperature, sea ice concentration, snow cover, surface albedo, phytoplankton concentrations, and atmospheric parameters such as cloud cover, temperature profile, ozone concentration, and aerosol levels, have been derived [3]. Additionally, studies have identified key observational indicators such as air temperature, permafrost thaw, snow cover, sea-land ice dynamics, and tundra ecosystem changes, which are critical for monitoring the Arctic's response to climate change [4].

When assessing the terrain feature changes, various methods are employed based on the research scope, categorized into two main groups: direct and indirect. Direct methods, including field and Global Positioning System (GPS) surveys, ice core drilling, and seismic approaches, are commonly applied for in-depth studies in specific areas. On the other hand, indirect methods, such as remote sensing, aerial surveys, Unmanned Aerial Vehicle (UAV) imaging, or numerical modeling, are typically utilized for research covering larger geographical areas.

Over recent decades, image-based approaches have advanced significantly, becoming crucial for monitoring global environmental changes. These methods provide valuable data for mapping, environmental monitoring, and long-term assessments with high spatial and temporal accuracy [5,6,7]. Its applications in geosciences are relatively diverse, including but not limited to shoreline changes, glacier retreats, vegetation cover, temperature variations, subsidence patterns, precipitation trends, and humidity fluctuations [8,9,10,11,12,13,14,15]. The growing prominence of image-based methods underscores their pivotal role in enhancing our understanding of dynamic Earth processes and climate-related transformations.

Svalbard, a Norwegian archipelago located in the Arctic, exemplifies climate change's effects in the Arctic. Covering an area of 62,248 km², Svalbard is notable for its extensive glaciation; glaciers cover 36,591 km² (59%) of its area [16,17]. Longyearbyen's administrative center is often cited as being on the frontline of climate change, emblematic of the drastic environmental shifts occurring in the world's northernmost town [18]. Over recent decades, Svalbard has experienced the fastest temperature increases within Europe, a phenomenon primarily attributed to the influence of the warm Atlantic Water current that significantly impacts local climate conditions by altering sea-ice coverage and water temperatures in the region's fjords [19,20].

Previous research on the glaciers in the Kaffiøyra region has contributed significantly to the understanding of glacier retreat, mass balance, and surface dynamics, employing various methods such as using historical and archival maps, field and GPS measurements, and more recently, remote sensing and UAV data [21,22]. These studies documented glacier area reductions and mass loss, offering valuable insights into the physical transformations over time. The Atlas of Changes in the Glaciers of Kaffiøyra [17] provided a comprehensive view of glacier transformations from 1909 to 2020, emphasizing glacier area and retreat. Further studies, such as those conducted by Sobota et al. [23,24] based on specific glaciers such as Waldemarbreen and Irenebreen, explored snow accumulation and melt rates but similarly lacked long-term climate correlations. Research tasks focused on ground temperature [25] and active layer thickness [26,27] have been vital in monitoring climate variables, though primarily focused on short-term observations and limited to summer periods.

The overarching stability of Svalbard's shoreline has been documented due to initial characteristics such as a rocky basement and ice-covered resistance to significant shoreline changes [28,29,30]. The localized studies reveal considerable variations along the West Coast, attributed to factors such as the changes in sediment supply, alterations in nearshore waves due to storms, and sea-level changes [31,32,33,34,35,36,37]. The Kaffiøyra shoreline and surrounding environment have been extensively studied, particularly by the Nicolaus Copernicus University research team, starting their first expedition in 1975. Their work has focused on various aspects of glaciology, geomorphology, and permafrost, including the impact of climate change on glaciers and shoreline dynamics. Many valuable studies have been focused on glacier dynamics and permafrost. The regional-scale and long-term variations of different terrain features in response to the temperature variations in the Kaffiøyra area were limited. The investigations on the outwash areas in the Kaffiøyra region have primarily focused on geomorphology and sediment transport, especially around the Waldemar River. Weckwerth et al. [38,39] detailed the dynamics of confined outwash fans shaped by glacial meltwater. Their studies emphasize how meltwater influences sediment deposition and landscape formation but lack long-term analysis of outwash evolution under changing climate conditions.

Based on previous investigations, it is clear that numerous valuable studies have examined individual and local-scale terrain features in Kaffiøyra. However, the observations are limited in characterizing the influences of temperature on the regional scale and long-term variations of the critical terrain features such as glaciers, shorelines, and outwash areas in the Kaffiøyra area. The gap presents an opportunity for a unified approach in conducting image processing to assess long-term variations of glacier termini, shoreline evolution, and outwash areas. Specifically, the study addresses the issues using a multi-decadal Landsat dataset from 1985 to 2023, enabling consistent, long-term analysis. By correlating these terrain feature changes with temperature data, this study aims to provide a more detailed understanding of how climate factors influence terrain feature dynamics, contributing crucial insights into local, regional, and global climate processes.

2. Materials and Methods

2.1. Study Area

The study area, Kaffiøyra, is located in the southwest region of Svalbard. The research conducted in this area will focus on three main types of terrain features to understand the impact of accelerated warming conditions on landform variations. Firstly, the study extracts and calculates the rate of shoreline changes to both land and glacier shorelines. Secondly, we plan to monitor and quantify the changes in the area and rate of seven glaciers, including the tidewater glacier - Aavatsmarkbreen, and six land glaciers: Waldemarbreen, Irenebreen, Elisebreen, Eivindbreen, Andreasbreen, and Oliverbreen. Lastly, we will assess the changes in the three outwash areas, which are zones where sediment is transported to the sea through the glacier's melting process during the summer seasons (Figure 1).

The geology of the Kaffiøyra plain is defined by two dominant formations. The western section comprises a tectonic trench separating Oscar II Land from Prins Karls Forland, while the eastern part is characterized by a crystalline bedrock folded and formed during the Caledonian orogeny. Over time, these early formations were overlaid by Quaternary deposits, consisting primarily of marine and glaciofluvial sediments such as marine pebbles, gravels, and sands. The surface of the region is now mainly covered by these Quaternary deposits. Marine sediments dominate the coastal areas, contributing to complex sedimentary dynamics significantly affecting shoreline stability and glacial retreat [38,40,41,42].

The river systems (outwash) on the Kaffiøyra plain are seasonal, primarily supplied by water melted from snow, glaciers, and permafrost, contributing significantly to sediment transport during the summer months [41]. The river systems primarily consist of networks that drain individual glaciers, forming distinct fluvial networks. These systems feature dense and intricate braided-distributary channels, particularly in the outwash fans and plains where sediment is deposited by glacial meltwater [38,43]. Precipitation in Kaffiøyra during the summer (July 21–August 31) is 40.7 mm on average but shows significant variations between years. For example, the wettest season recorded 122.5 mm (1997), while the driest was only 8.5 mm (2010). Interestingly, dry summers are more common in Kaffiøyra, especially in the later years of the study period (1975–2010) [44]. These seasonal and climatic variations play a significant role in the hydrological and geological dynamics of the region, influencing both sediment transport and the long-term retreat of glaciers.

2.2. Data for the Site

Landsat images and GIS techniques extracted key terrain features, including shorelines, glaciers, and outwash areas. Landsat's long-term archive makes it an ideal selection for monitoring surface changes because of the consistent data sources and resolution. Other satellites, such as Sentinel-2, might offer higher resolution, but the duration of available data is limited. The Landsat multi-decadal coverage provides a 30-meter resolution that could balance the detail and coverage, which is crucial in Arctic regions with frequent cloud cover. Additionally, the free, well-calibrated data offered by Landsat is cost-effective, making it a practical choice for large-scale, long-term environmental studies [45].

The study utilized multitemporal remote sensing data of the Landsat Thematic Mapper (TM), Landsat Enhanced Thematic Mapper (ETM+), and Landsat Operational Land Imager (OLI) from 1985 until 2023 (see Table 1). The data are available for free at the archive of the U.S. Geological Survey. The dataset spans from 1985 to 2023, including Landsat 4-5 TM, 7 ETM+, 8 OLI/TIRS, using the WGS_1984_UTM_Zone_33N reference system and targeting annual coverage. However, due to various challenges, such as the Landsat 7 scan line errors and cloud cover, usable images were obtained for 29 out of the 38 years (Table 1).

2.3. Methods

The primary methodology for this research involves integrating remote sensing and GIS techniques to analyze Landsat images. Table 2 lists the main steps and the associated toolboxes to analyze the observed images. The following is a brief description of the processes for our analyses.

• Data collection:

The study area is in the polar region, typically covered with snow and ice from October to May next year. The Landsat images taken in this period cannot be used for the study. Therefore, the summertime images taken from June to September are used to extract the terrain features of shorelines, glaciers, and outwash areas in different years.

The image quality is critical. Thus, image selection will consider factors such as cloud cover, scanline error in satellite images, and shadow, which might cause difficulty extracting objects. From 1985 until 2023 (38 years), the study has extracted 29/38 shorelines, 28/38 glacier termini, and 21/38 outwash based on the available images.

• Image processing:

Based on the selected images, there are several processing steps, including geometric correction to align and position the image accurately, image cropping for the area of interest, and color adjustment to enhance boundary visibility.

• Boundary extraction:

The refined images are then used to delineate three terrain feature boundaries, i.e., the shoreline (sea-land boundary), glacier outline (glacier-land and glacier-sea boundaries), and outwash outline (outwash-land boundary). The Landsat Toolbox is employed to assist in distinguishing between different objects (e.g., vegetation, water, glaciers, land, lake, and others) based on their spectral reflections, categorizing them into classes (often ten classes). These are then simplified into two primary groups, the interest group and the remaining group, and then the boundary is obtained from those two groups [46]. Manual digitization supplements the automated processes in cases where clouds or shadows cover the boundaries. For more detailed processes, one can refer to Table 2.

The shoreline positions can be referenced from any consistent linear feature, such as the dry/wet line, the vegetation line, the low-water line, or the high-water line. The shoreline can be digitized from various sources, such as satellite imagery, digital orthophotos, and historical coastal survey maps). This study detected the shoreline based on consistent position and source, the dry/wet line, and the source of satellite imagery. Each shoreline vector demonstrates a specific position in time and must be assigned a date in the shoreline attribute table.

• Changes analysis

Once the boundaries are established, they are analyzed over time. The areas of glaciers and outwash are calculated to determine how they have varied over the years. The shoreline data is processed through the Digital Shoreline Analysis System (DSAS) [47,48]. It is an automated method for establishing measurement locations, performing rate calculations, and providing statistical data necessary to assess the robustness of the rates and distances. In the study area, both land-shorelines and glacier-shorelines exist, which may raise concerns about DSAS effectiveness. However, the DSAS manual states that it can be used for any boundary-change problem with clearly identifiable positions over time, such as glacier limits or riverbanks [48].

The required inputs to calculate the shoreline change through DSAS include shorelines (Figure 2 - green, yellow, blue, red lines) and baseline (Figure 2 - black lines). The baselines in this study were created manually, with several iterations needed, due to the irregularity of the shoreline in the Kaffiøyra plain. While the goal was to create transects perpendicular to the shorelines, the complexity of feature boundaries meant some transects were not perfectly aligned. However, every effort was made to ensure the transects fit as closely as possible. For the simpler shoreline profile in other zones, the baseline creation was more straightforward, with parallel alignment to the shoreline being sufficient. The baseline serves as the starting point for all transects cast by DSAS. It can be drawn anywhere relative to the shoreline vectors: onshore, offshore, and midshore among the shorelines. The DSAS uses a measurement baseline method [49] to calculate distance and rate-of-change statistics for a time series of shorelines. The transects will be cast perpendicular from baseline to shorelines (Figure 2 - gray line). The user can set the smoothness of the transect, distance from the baseline, and transect interval (Table A3).

After DSAS creates a transect cast from baseline to shorelines, the system will provide statistical analyses to calculate shoreline changes such as net shoreline movement (NSM), shoreline change envelope (SCE) (Figure 2), end point rate (EPR), linear regression rate (LRR) (Table 3). These parameters are calculated for each transect.

Notice that the NSM measures the greatest distance between the earliest and the latest shoreline, and its values can be negative or positive depending on whether it is eroded or accreted. The SCE measures the greatest distance between all shorelines without counting on the years, and the values are always positive. The EPR is calculated by NSM divided by time, and it is simple to calculate with only two shorelines sufficient, but the value will miss the trend or magnitude if you have more than two shorelines [50,51]. The LRR is estimated by fitting a least-squares regression to all intersections on each transect. Therefore, the LRR will provide a more comprehensive view of the trend over many years.

2.4. Uncertainty and Data Accuracy

Several factors contribute to uncertainty in the change analysis, mainly when using Landsat imagery. This study predominantly utilizes Landsat L1TP (Level-1 Terrain Precision-Corrected Product) and L2SP (Level-2 Surface Reflectance Product), which are terrain-corrected and orthorectified for high positional accuracy. These corrections are essential for high-latitude regions like Svalbard, where distortions from Earth's curvature and atmospheric conditions can be significant. The georeferencing error for these products is estimated at 12–15 meters, ensuring reliable spatial alignment throughout the dataset [52].

In Svalbard, the tidal range is small, typically 1–2 meters, with no significant tidal flats, limiting tidal influence on shoreline displacement. Studies show that climate-related factors such as glacial retreat, permafrost thawing, and increased wave energy due to longer ice-free seasons are the primary drivers of coastal erosion [28,53]. Additionally, the loss of sea ice has led to more decisive wave action and storm surges, further reducing the role of tides in shoreline changes [54]. Based on these studies, we assume that the tides have a negligible effect on shoreline changes in the Kaffiøyra region. Therefore, we might receive 15 meters of uncertainty in shoreline delineation due to the 30-meter resolution of Landsat imagery (half a pixel) [55,56].

We calculate total shoreline position uncertainty [57] by combining georeferencing and spatial resolution errors using the formula:

$$U=\sqrt{{\left(GeoreferencingError\right)}^2+{\left(SpatialResolutionError\right)}^2}$$

(1)

Assuming 15 meters for each error, the total $U$ is about 21.2 meters. This uncertainty is consistently applied throughout the analysis to ensure confidence in shoreline changes. Additionally, using the 95% confidence interval, incorporated into rate-of-change statistics [47], could provide a robust statistical measure of uncertainty and ensure that the derived shoreline change rates are interpreted with appropriate consideration of potential errors and enhance the reliability of the results.

3. Results

3.1. Analyze Shoreline Changes by DSAS

Figure A1 shows an example of the 2014, 2015, 2022, and 2023 boundary extraction results. Once the image extraction is complete, all the shorelines used as input for DSAS are overlapped (Figure 3a). The baselines are created along the shorelines. A baseline is good only if the transect starts from the baseline and must intersect all shorelines we have. Due to the different shoreline geometry, this study divided the baseline into five segments, as the white lines shown in Figure 3a. The transects are cast perpendicular from the baseline to all shorelines. In different baselines, we can set various parameters for the transects, such as the interval between two transects or the length of the transect (Figure 3b).

An overview of Figure 3a, the overlap shoreline map shows that the shoreline of the study area has erosion and deposition but is insignificant. Except for two regions with glacier terminus bordering the sea, the shoreline is dramatically eroded.

The study considers two main parameters (i.e., NSM and LRR) to interpret the results. The results relevant to the remaining parameters (i.e., SCE and EPR) will be presented in the appendix (see Figure A2 and Figure A3 for the calculated SCE and EPR). Figure A4 shows the specific values for each transect. Figure 4 shows the results of the calculated LRRs on different transects. The values in Figure 4 are shown with the representative color ramps. The dark red transects show values between minimum to -3 m/yr, meaning that the shoreline is experiencing intense erosion. However, the dark blue color represents a strong accretion for an area.

Figure 5 further zooms into the local areas to show the temporal behavior of the shoreline variations. Table 4 summarizes each zone's LRR and NSM values and the associated percentage of erosion and accretion for each zone. The shorelines in zones 1, 3, and 5 appear relatively stable. Zone 1 is situated in the northern part of the study area and has a shoreline length of approximately 4.5 km. Here, 255 transects were created, revealing annual erosion and accretion rates ranging between -2.2 ±0.62 m/yr and 2.3±0.73 m/yr, with an average shoreline change of -0.1 m/yr. In this zone, erosion is slightly more prevalent than accretion, affecting 56.9% of transects, compared to 43.1% showing accretion.

Zone 3 is centrally located within the study area, spanning roughly 21 km of shoreline. This zone is well-documented with 1,058 transects, demonstrating erosion and accretion rates varying from -1.1 ±0.76 to 3.7 ±1.62 m/yr and an average LRR of 0.5 m/yr. Accretion is the dominant process here, as indicated by 76.7% of transects, with only 23.3% showing erosion.

Zone 5 is located in the southern part of the study area, with a shoreline length of about 5.5 km. The number of transects is 305. The erosion and accretion rates here vary from -1.7 ±0.87 to 1.7 ±0.71 m/yr, with an average deposition rate of 0.2 m/yr. Accretion predominates in this zone, affecting 77.7% of transects, while the remaining 22.3% exhibit erosion.

In summary, the results of LRRs indicate that the shoreline on the glacier outwash plain is relatively stable compared to the LRRs of the tidewater glaciers. These zones are the boundary between the outwash plain and the sea. The variation of LRRs in these zones highly depends on the geometry of the coastal line and the sources of sediments carried by water from melted glaciers and snow. In addition, ocean currents, wind, sea waves, and floating ice could be the driving forces to control the overall behavior of the coastal line. Detailed dynamics of the variation along the coastal line might need further modeling and the associated on-site observations to quantify the erosion and deposition processes.

In contrast to the three zones with relatively stable shorelines, zones 2 and 4, which host tidewater glacier termini, have experienced significant changes, as shown in Figure 5 and Table 4. Zone 2, which encompasses the Aavatsmarkbreen tidewater glacier terminus shoreline, lasts about 8.3 km. A total of 122 transects were created here, showing a wide range of erosion or retreat rates from -64.7 ±9.37 to 9.2 ±2.92 m/yr, with an average retreat rate of -25.4 m/yr. Erosion is the dominant process, evident in 70% of the transects, particularly concentrated in the central terminus, and 30% of transects do indicate accretion.

Zone 4 is the Dalhbreen tidewater glacier terminus, with the shoreline length of approximately 3 km. Here, 152 transects were created, and the erosion and accretion rates ranged from -39.0 ±4.97 to 0.7±0.78 m/yr. The average shoreline change rate is significantly erosive at -27.5 m/yr, with erosion evident in 96.7% of transects.

There are also relatively high LRR areas near the corner of the glacier valleys. The coastal weather, ocean currents, or floating ice from melted glaciers and sea ice might influence these erosion areas. The retreating glaciers might also lead to barren zones that appear downstream (shown with light or dark blue colors) along the coastal lines and beside the glacier-sea boundaries. We observe the relatively high deposition areas along the barren zones, where the drifts accumulate, or the rocky outcrops that are covered by glaciers.

The Aavatsmarkbreen LRR indicates significant retreat, consistent with the 360 to over 540 meters of retreat reported [58] for the period 2000–2006 and the ongoing retreat documented in the Atlas of Changes in the Glaciers of Kaffiøyra [17]. Similarly, the Dalhbreen LRR aligns with the 310 meters of retreat recorded [58] from 1986 – 2005, with continued retreat also confirmed by the Atlas [17]. These comparisons underscore the accelerating glacier retreat and highlight the results of this study that match the previous ones.

In Figure 6 and Table 4, the NSM illustrates the distances between the youngest (2023) and oldest (1985) shorelines along each transect. It is observed that the values in zones 1, 3, and 5 vary from -56 m to +215 m. The negative values are marked by light red to dark red, meaning shoreline erosion has occurred. The longest eroded distance is -56 m to the landward in zone 5 (dark red). Meanwhile, the NSM value is positive at the transects marked by light blue to dark blue, which means the shoreline is accreting. The greatest distance is 215 m to the seaward in zone 5 (dark blue). The NSM values range from -2,762 m to + 426 m between the youngest (2023) and the oldest (1985) shorelines in zone 2 and zone 4. In zone 2, the longest transect with a value of -2,762 meters to landward, indicating that in 38 years, the Aavastmarkbreen has retreated the farthest distance of 2,762 meters (dark red). In zone 4, the glacier has retreated to the farthest distance of 1,654 meters (dark red).

3.2. The Glacier Terminus Areas, Retreat Rates, and the Correlation to Temperature Variations

3.2.1. Area Changes in Glacier Termini

The glacier terminus is the toe or snout of a glacier, and monitoring its changes is vital for assessing the retreat or growth of the glacier, in other words, the impact of climate change. The research utilized Landsat imagery to delineate the areas of seven glaciers within the study region, which comprised one tidewater glacier (Aavatmarkbreen) and six on-land glaciers (i.e., Waldemarbreen, Irenebreen, Elisebreen, Eivindbreen, Andreasbreen, and Oliverbreen), as shown in Figure 7. The extraction period spanned from 1985 to 2023, although data for some years were unavailable due to image errors or cloud cover, resulting in successful extraction for 28 of the 38 years. The terminus area changes were calculated by creating a reference line called the glacier baseline, shown in Figure 7, across the glaciers (yellow line). The baseline for calculating glacier terminus area changes was determined to ensure accuracy and consistency. The baseline was positioned behind all recorded glacier termini and tended to be perpendicular to the glacier's flow to capture advances and retreats. The baseline remains fixed across all years, allowing for consistent measurements of terminus positions.

Figure 8 and Table A1 illustrate a clear retreat trend among all seven glaciers. Results show that the area of the Aavatsmarkbreen glacier terminus decreased from 11.05 km² in 1985 to 3.63 km² by 2023 based on the defined baseline. This reduction of 7.42 km² translates to a loss of 67.1% of its total area, leaving only 32.9% of its original size. Aavatsmarkbreen experienced the most significant decline in the area. The on-land glaciers experienced varying degrees of area loss: Waldemarbreen at 45.2%, Irenebreen at 38.7%, Elisebreen at 35%, Eivindbreen at 26%, Andreasbreen at 44.5%, and Oliverbreen at 56.3%, as shown in the percent of remaining area chart in Figure 8. Note that the reference areas of the original glacier in the calculations are those obtained in 1985.

The study shows an average land glacier area loss of 41%, while the average remaining area is 59% by 2023. Such results are close to 43% of area reduction reported in previous studies [17,22], observed from the late 19th century to 2015. Despite the different periods, the percentage similarity suggests a consistent long-term trend of glacier retreat for land glaciers.

3.2.2. Aavatsmarkbreen – Tidewater Glacier Surge Event

The long-term time series data are valuable for assessing the dynamic of glaciers in specific areas. The available images and resolution have enabled the identification of historical events for the important terrain features in the study area. Significant fluctuations in the area of Aavatsmarkbreen were observed in 2015 (Figure 8), with an expansion of 2.5 km² from 2013, marking a substantial deviation from the typical retreat pattern. The Aavatsmarkbreen surge events in 2014 and 2015 have been previously documented, with glacier surges being a well-researched phenomenon in Kaffiøyra [59,60]. These surges are driven by complex interactions of thermal, hydrological, and mechanical processes, wherein surface meltwater drains into the glacier base, causing basal lubrication that reduces friction and increases glacier velocity. During the 2014–2015 surge, Aavatsmarkbreen expanded from 5.57 km² in 2013 to 6.52 km² in 2014, peaking at 8.07 km² in 2015 (Figure A5). On average, before the surge, the glacier retreated at 0.24 km²/yr between 1985 and 2013. During the surge period, the area gained at a rate of 1.25 km²/yr before resuming retreat at a rate of 0.53 km²/yr from 2015 to 2023 (Figure A6). These surge events are consistent with the documented behavior of Aavatsmarkbreen and reflect a known pattern of periodic glacier advances followed by retreat.

Smaller increases were observed in both tidewater and land glaciers during 1987, 1990, 2001, 2014, and 2016. Upon investigating these periods (Figure 7), it was evident that these increases were not the result of surge events, as the other glaciers continued their steady retreat. The anomalies in these glaciers primarily stem from changes in the lateral areas. Since the satellite images were captured during the summer, snow cover in the lateral regions may have temporarily increased the measured glacier area from June to September. However, the glacier termini remained consistently clear, ensuring accurate tracking of the main glacier boundaries. These minor anomalies, therefore, reflect a limitation of the method, as seasonal snow cover in the lateral areas can influence the apparent glacier extent.

3.2.3 Retreat Rate in Glacier Termini

In addition to assessing changes in areas, evaluating the rate of glacier retreat is crucial. The study generated centerlines based on the glacier outlines in 1985 (see Figure 7, marked with white lines). One centerline (AA') and two additional reference lines (BB' and CC') will be established for the tidewater glacier Aavatsmarkbreen. For the land glaciers, one centerline is defined for each glacier (i.e., Waldemarbreen, Irenebreen, Elisebreen, Eivindbreen, Andreasbreen, and Oliverbreen). These centerlines will serve as the basis for tracking retreats over time. The study will employ a linear regression model to analyze all intersection points where the centerlines meet the glacier outlines at different times. By fitting these data points, we can derive a rate of retreat for each glacier at a specific time.

In long-term retreat, Table 5 and Figure A7 show that the tidewater glacier Aavatsmarkbreen exhibits the most rapid retreat along its centerline, labeled AA', with a rate of 62.93 m/yr. Note that the reference lines BB' and CC' show similar patterns but with relatively slow retreat rates. Such slow variation of retreat along the reference lines might be influenced by the land of the coastal line and the glacier's geometry. Among the land glaciers, Andreasbreen shows the highest retreat rate at 43.74 m/yr. The retreat rates for the other land glaciers decrease as follows: Irenebreen at 22.96 m/yr, Elisebreen at 19.48 m/yr, Oliverbreen at 18.99 m/yr, Waldemarbreen at 9.24 m/yr, and Eivindbreen at 7.03 m/yr. The average retreat rate of Aavatsmarkbreen is 43.74 m/yr, while the average rate for land glaciers is 20.24 m/yr. It is worth mentioning that the retreat rate of Elisebreen at 19.48 m/yr in the study is close to the average value of 18 m/yr proposed by Sobota et al.(2016) [22]. The results of Sobota et al.(2016) were based on the observations from 1909 to 2015. Furthermore, the retreat rate of Waldemarbreen recorded by Sobota et al.(2016) was 8 m/yr, which is also close to the observation in the study. The verification shows that the image-based approach could be consistent with the ground measurements.

In short-term retreat, Table 5 shows that Aavatsmarkbreen continues to exhibit rapid retreat along AA', with an average of 83.97 m/yr over the 10-year recent boundaries, rising sharply to 210.48 m/yr in 2019-2023. The BB' and CC' show relatively lower rates, with BB' increasing from 33.49 m/yr (10-year) to 75.13 m/yr and CC' shifting from -3.16 m/yr to 165.57 m/yr, likely influenced by coastal land and glacial geometry. Among land glaciers, Andreasbreen has the highest short-term retreat, at 98.94 m/yr (10-year) and 101.02 m/yr in 2019-2023. Other land glaciers show lower rates, with Irenebreen, Elisebreen, Oliverbreen, Waldemarbreen, and Eivindbreen showing varying increases, reflecting a recent acceleration in glacier retreat. These findings align with long-term trends.

Figure 10 illustrates how the distance between measurements taken over two-year intervals varies along each centerline of a glacier. The boxplot chart helps to assess the overall behavior of different glaciers and whether a glacier's retreat is accelerating or decelerating along specific transects. Significant surge events dramatically impacted its measurements in the case of Aavatsmarkbreen (see the results of AA', BB', and CC'). To accurately assess the glacier's retreat behavior, we excluded data from the years affected by surges (2014-2019). This exclusion is vital because surges characterized by rapid and periodic advances might distort long-term trends and are generally not observed from land glaciers in the study. By focusing on non-surge periods, we aim to provide a more explicit analysis of the gradual changes indicative of the glacier's response to climatic conditions, allowing for more consistent comparisons with retreat patterns of land glaciers.

The boxplot for centerline AA' reveals a wide dispersion of data points, indicating substantial year-to-year variability in glacier retreat, suggesting that the retreat pattern along this transect is highly irregular and variable. In contrast, reference lines BB' and CC' show progressively minor variability, indicating a more consistent retreat pattern, particularly at CC', where the behavior is the most stable among the tidewater glaciers. Because the two reference lines are close to the coastal land, the geometry of the glacier or glacier flow patterns could constrain the variation of these reference lines. In addition, the coastal lands could be essential in binding two sides of the glacier front.

When we compare tidewater glaciers to land glacier centerlines, the land glaciers generally exhibit a more stable retreat. The result is evident from their boxplots, which show a narrow range and data clustered close to the median. Waldemarbreen, Eivindbreen, and Oliverbreen demonstrate the typical stable retreat patterns among the land glaciers. Conversely, Andreasbreen exhibits the most irregular retreat behavior, standing out for its inconsistency within the group of land glaciers.

The observable retreat of glaciers nearly the past four decades indicates climate change's impact on the Arctic. Various factors, including changes in precipitation, the albedo effect, ocean temperature increases, shifts in atmospheric circulation, deposition of black carbon, thawing of permafrost, increased iceberg calving, and the loss of sea ice, all contribute to glacier retreat. However, the principal driver among these is the temperature change [61,62]. Temperature changes might directly and indirectly interact with and impact other factors. Thus, the following section will discuss the temperature changes within the study area.

3.2.4. Temperature Increasing Trend

In this study, we considered Ny-Ålesund temperature data due to the discontinuous records from Kaffiøyra. The Ny-Ålesund monitoring station (78.9230N, 11.9331E), located about 25 km north of Kaffiøyra, provides a detailed temperature record from 1975 to the present, capturing monthly variations and serving as a critical resource for understanding regional climate dynamics. Research has demonstrated a high correlation between temperatures at Ny-Ålesund and Kaffiøyra [63]. Therefore, using Ny-Ålesund temperature data in this study could be acceptable. The data is accessible via the NASA Goddard Institute for Space Studies (GISS).

The annual mean temperature chart in Figure 11 illustrates a gradual upward trend from 1975 to the present, with an approximate average yearly increase of 0.09 ºC. This data underscores a consistent warming pattern over the years.

Analyzing the monthly temperature records from the Ny-Ålesund station, alongside Landsat imagery, reveals that the study area experiences summer during June, July, and August. During these months, surface temperatures rise above 0 ºC, and snow cover is typically absent. Figure 12 displays the mean summer temperatures spanning from 1975 to 2023 and indicates a warming trend. The data show an average annual increase in summer temperatures of 0.04 ºC. The average summer temperature in Ny-Ålesund closely matches that of Kaffiøyra, with a difference of approximately 0.04 ºC compared to 0.03 ºC. This slight difference may be due to the summer period defined in the study of Kejna and Sobota (2019) [63], which considers summer from July 21 to August 31. Notice that the variation behavior of annual mean temperatures might not be consistent with the summertime mean temperature. An obvious example is that in 2020, we found a historically high mean summer temperature, but the annual mean temperature in 2020 was relatively low compared to previous years. In 2015 and 2016, the high annual mean temperature came with the high mean summer temperature.

These charts reveal the warming temperature trends in the study area. The following part of the study will quantify the correlation between increasing temperatures and the glacier's remaining area by calculating the correlation coefficient of the datasets. This effort seeks to elucidate the specific impact of temperature increases on glacier retreat, offering a clearer understanding of the dynamics at play in temperature-induced changes to glacier area variations.

3.2.5. Correlation Between Glacier Terminus Area and Summer Time Mean Temperature

Before establishing the correlation between temperature increases and the reduction in the glacier terminus area, the study initially normalized each variable through z-score normalization [64]. This process was essential to ensure fair and accurate comparisons between the two variables. Z-score normalization adjusts the data to have a mean of 0 and a standard deviation of 1, thereby allowing for a direct comparison of different datasets on a like-for-like basis, regardless of their original scales or units of measurement.

This preparatory step is crucial for producing reliable results and analyzing the relationship between temperature variations and the remaining area.

$$z={\displaystyle \frac{x-\mu }{\sigma }}$$

(2)

Where $z$ is the z-score, $x$is the value of the data point, $\mu$represents the mean of the dataset, and $\sigma$is the standard deviation of the dataset.

In this analysis, we utilize the Pearson correlation coefficient $\left(r\right)$ [65] to measure the linear relationship between independent variables (temperature) and dependent variables (remaining glacier area). The formula to calculate the Pearson correlation coefficient is:

$$r={\displaystyle \frac{\sum \left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}{\sqrt{{\sum \left(x_i-\overline{x}\right)}^2{\sum \left(y_i-\overline{y}\right)}^2}}}$$

(3)

where the notations $x_i,y_i$ are the individual sample points indexed with $i$, and $\overline{x},\overline{y}$ are the means of $x,y$ samples, respectively. Figure 13 shows the scatterplot of the correlation between summertime mean temperature and the glacier's remaining areas after using the z-score normalization. The study proposes a linear regression equation to evaluate the relationship between two variables.

The findings indicate that all seven glaciers exhibit a negative correlation with temperature, with values ranging between -0.65 and -0.73 (see Table 6). These correlation coefficients are significantly negative (approaching -1), suggesting a relatively strong inverse relationship. The remaining glacier terminus areas tend to decrease as temperatures increase.

3.2.6. Prediction of Glacier Area Changes

In this section, the study will utilize a linear regression equation to predict the changes in the glacier terminus areas, drawing on the proven linear correlation between temperature increases and glacier retreat. The prediction, covering the period from 2020 to 2039, will incorporate temperature projections derived from the Shared Socioeconomic Pathways (SSPs) of the Svalbard archipelago. These SSP scenarios comprehensively account for how socioeconomic factors influence greenhouse gas emissions and the strategies for adaptation. In addition, the study uses future temperature prediction data from the Ny-Ålesund station (Figure 12 – marked with green rectangular symbols) to make the parallel comparison for the same period. This comparison aims to clarify how glacier areas might evolve in response to regional SSP-derived temperatures and local temperature trends at the Ny-Ålesund station.

Temperature data for the SSP scenarios were sourced from 2020 - 2039, as detailed on the Climate Change Knowledge Portal. The result involved collecting temperature data for each year's summer months (June, July, August) within this period. The process entailed identifying the minimum, maximum, and average temperature values across these months for 2020 to 2039. Subsequently, the averages of these minimum (Min.), maximum (Max.), and average (Avg.) monthly temperatures were calculated to serve as the input data. In the study, the temperature data for Ny-Ålesund from 2020 to 2039 was analyzed using the linear equation provided based on the station's summer data (Figure 12) to determine the minimum, maximum, and average temperature values for this period.

Figure 14 shows the future temperature predictions and the remaining glacier areas for 2020–2039. In Figure 14, we have observed that the local summer temperature at the Ny-Ålesund station exhibits a much higher temperature increase group than that obtained from the SSPs of the Svalbard archipelago. The predicted glacier area changes could be considerably different between SSPs and the local observations at the Ny-Ålesund station. Results presented in Table 7 show that the temperature projections under SSP scenarios vary between 0.44 and 4.32 ºC. Consequently, the predicted glacier area remains significantly large, with estimated ranges varying from 8.6 to 16.25 km². This projection contrasts sharply with the image-based measured remaining area of Aavatsmarkbreen from 1985 to 2023. Based on our observations, the remaining area of Aavatsmarkbreen is 3.63 km² by 2023 (see Figure 8). The temperature data from the Ny-Ålesund station yields results that align more closely with the actual changes observed at the Aavatsmarkbreen glacier (see Figure 14 and Table 7). However, all the predictions considerably underestimate the high variations of glacier remaining areas. The example of Aavatsmarkbreen in the study suggests that the impact of climate change, particularly the effect of rising temperatures, on glacier retreat is more severe at the local scale. The findings highlight the crucial role of localized temperature measurements in accurately assessing climate change impacts on glacier dynamics. Establishing additional monitoring stations will be vital for enhancing the precision of future forecasts and developing a more comprehensive understanding of climate change's local effects.

3.3. Variations of Outwash Areas

The retreat of glaciers reduces the extent of the glacier terminus and significantly impacts the surrounding landforms. The outwash is one of the main terrain features that might be developed with the land glacier retreat. The outwash areas in the study area are identified and shown in Figure 1. Specifically, Outwash 1 is located downstream of Waldemarbreen, while Ourwash 2 might receive melting water and sediments from Irenebreen and Elisebreen. Outwash 3 is located downstream of Eivindbreen, Andreasbreen, and Oliverbreen. There is no clear interface to differentiate the sources of water and sediment depositions from Oliverbreen, Eivindbreen, and Andreasbreen.

Figure 15 and Table A2 illustrate the changes in the areas of three large outwash areas over the past 38 years. Based on the original sizes of the outwash areas estimated in 1985, Outwash 1 experienced a decrease of approximately 13.1% in area, while Outwash 2 and 3 showed increases of 6.7% and 9.6%, respectively. There has been a notable decrease in Outwash 1 since 2019. A significant drop in the Outwash 1 area appeared in 2020 when the highest mean summer temperature was obtained (See Figure 12). We could not find a specific correlation between the Waldemarbreen remaining area and the Outwash 1 area changes. However, the Outwash 1 area reduction could be induced by the disconnected ice source of Waldemarbreen based on the image shown in Figure 1.

Outwash 2 and 3 show notable increases since 2015. The rise in outwash areas might be relevant to the relatively high temperature near the Ny-Ålesund. In 2015, the summer and annual mean temperatures reached historically high values compared to the data before 2015. Again, no clear evidence exists to correlate the glacier's remaining areas with the outwash areas. Unlike the conditions in Outwash 1, the increase in summer and annual temperatures are the possible factors that drive the rise in outwash areas for Outwash 2 and 3 if the upstream glaciers have sufficient ice supplied to their sources.

Water sources and sediments from glaciers or snow could influence the changes in outwash areas. As the glacier retreats, the increased meltwater output may expand outwash plains where the water and sediments are deposited. However, the time series data in Figure 15 show a slight increase in outwash area changes from 1985 to 2023. These changes are consistent with the dynamics of glacier retreat; as the glacier terminus recedes due to increased melting, the volume of meltwater feeding the outwash areas changes.

4. Discussion

The analysis reveals a significant increase in long-term shoreline change rates across different coastal zones, categorized into land-shorelines (Zones 1, 3, 5) and glacier-shorelines (Zones 2, 4). Sediment dynamics and hydrodynamic conditions predominantly influence the land-shorelines, whereas climatic factors primarily affect glacier-shorelines. Prins Karls Forland Island is a noteworthy feature within this landscape, located approximately 15 km from Kaffiøyra (Figure 1). This island is a natural buffer, shielding the nearby land-shorelines from the Atlantic's hydrodynamics. This protection is particularly evident in Zone 3, where we observe substantial accretion. Here, the average LRR over a 21 km stretch is approximately 0.5 m/year, and the NSM indicates an advance of about 56 m towards the sea over almost four decades. This advancement can be attributed to the island's shielding effect and the abundant sediment supplied from the outwash areas in the area, as evidenced by the accretion transect distribution detailed in Figure 6 (zone 3).

Conversely, the glacier-shorelines in Zones 2 and 4, known as Aavatsmarkbreen and Dalhbreen, experience a significant retreat, influenced mainly by climatic conditions, highlighting that tidewater glaciers are distinctly more vulnerable under the same climate conditions as land glaciers. Tidewater glacier terminus has retreated by 67.1% of the area, substantially more than the 40.9% average retreat observed inland glaciers. Additionally, the average retreat rate (i.e., from AA', BB', and CC') for the tidewater glacier (43.74 m/yr) is roughly twice that of land glaciers (20.24 m/yr). Recent data also show accelerated retreat rates. The Aavatsmarkbreen AA' rose from 83.97 m/yr over the 10-year boundary to 210.48 m/yr in 2019-2023, with BB' and CC' showing similar but smaller increases. Land glaciers follow this trend, with Elisebreen increasing from 8.10 m/yr to 60.69 m/yr and other glaciers like Oliverbreen and Waldemarbreen showing notable rises. This widespread acceleration suggests glaciers respond more intensely to recent climate changes, signaling a rapid shift in ice loss dynamics.

Based on previous reports, surface events in the study area occurred twice, first in 1987 and then again in 2014 and 2015, exclusively affecting the tidewater glaciers. While these events led to advances toward the sea and increased area at the glacier terminus, they should be viewed as a consequence of climate change. The development of crevasse systems during these periods likely facilitated meltwater flow to the glacier base. Following these surge events, the area loss rate of the glacier terminus more than doubled, indicating increased instability in the glacier structure post-surge, which made the terminus more prone to collapse and accelerated area loss. When examining retreat behavior via the boxplot chart (Figure 10), it is evident that the tidewater glacier exhibits more significant variability than land glaciers. The distance between two years is often significant and spans a wide range of values, reflecting this variability. In contrast, land glaciers demonstrate more stable retreat patterns; the distance between two years is smaller and distributed close to the median. This analysis supports that the tidewater glacier is more susceptible to rapid changes and potentially more hazardous in consistent climate conditions.

Temperature analyses further underscore this vulnerability, showing a strong negative correlation (between -0.65 to -0.76) between summertime mean temperatures and glacier remaining areas, signifying that rising temperatures are a primary driver of glacier retreat. A predictive model based on linear regression and temperature data from Ny-Ålesund closely aligns with observed changes more accurately than regional climate models, which fail to capture the current status of the Aavatsmarkbreen and deviate from observed data. This discrepancy emphasizes the need for improved monitoring systems that can provide more accurate local data to understand better the effects of climate change on local and regional scales.

Furthermore, outwash plains 2 and 3 show a growth trend, though the increase is less than 10%. These outwash areas are critical for shoreline evolution in Zone 3, serving as significant sediment sources. The gradual expansion of these outwash areas aligns with the steady melting of ice and snow during the summer, providing a consistent supply of sediment that supports sustained shoreline accretion. In the study, Outwash 1 area shows a decreased behavior triggered from 2019 to 2020. The image-based assessment could give an overview of the area changes, but the correlation between the outwash and the glacier is unclear. More detailed, accurate, on-site monitoring and estimations of glacier and outwash evolution might improve understanding of the interaction mechanism between glacier and downstream outwash.

Lastly, the study recognizes several limitations. While the primary focus is on nearly 40 years of long-term observation, the 30-meter resolution of Landsat imagery poses challenges in capturing small-scale landform changes. Additionally, the manual creation of baselines and transects introduces a degree of subjectivity, particularly in complex areas. The reliance on summer imagery further limits the consideration of seasonal variability in glacier and shoreline dynamics. Future research could address these limitations by utilizing higher-resolution remote sensing imagery or UAVs for more detailed analysis of small-scale changes. Incorporating numerical modeling could also enhance understanding of seasonal variability.

5. Conclusion

This study has comprehensively explored the impacts of climate change on various landforms from 1985 to 2023, utilizing cost-effective remote sensing and GIS techniques to monitor spatial and temporal changes. Key findings include the relative stability of the land-shoreline, with average changes ranging from -2.2 to +3.7 m/yr. Considerable changes were observed in the tidewater glacier terminus under the climatic conditions in the past 4 decades. The results show that the Aavatsmarkbreen exhibited the greatest loss in the terminus area, with a reduction of 67.1%, in contrast to an average area loss of 40.9% among land glaciers. Regarding the retreat rates, Aavatsmarkbreen also displayed the most rapid retreat behavior, averaging 62.93 m/yr, compared to 43.74 m/yr for land glaciers. In recent years, retreat rates have increased dramatically. For example, Aavatsmarkbreen AA' increased from 83.97 m/yr (10-year recent boundaries) to 210.48 m/yr (2019-2023). This sharp increase underscores the intensified impact of climate change on glacial dynamics, particularly on tidewater glaciers, which are retreating at accelerating rates due to heightened sensitivity to warming temperatures.

Additionally, temperature analyses from the Ny-Ålesund station revealed a strong negative correlation (ranging from -0.65 to -0.72) with the remaining glacier areas, indicating a direct impact of rising temperatures induced glacier retreat. Furthermore, the study employed these correlations alongside regional climate projections from Svalbard's SSPs and linear temperature trends from Ny-Ålesund to predict future declines in Aavatsmarkbreen. These projections suggest that regional climate models may fail to capture local scale variations, potentially leading to underestimations of climate impact.

As for the outwash plains, there is limited data availability (21 out of 38 Landsat image sets). It was noted that two outwash areas, Outwash 2 and 3, exhibited moderate increases in size, though less than 10%.

Overall, this study underscores the responses of Arctic terrain features to climate change and highlights the need for finely scaled models to predict future environmental transformations accurately. It overviews long-term landform behavior under accelerated climate change, forming a foundation for future hydrodynamic and glacier numerical modeling research. These models will further explore the interactions between climate dynamics, shoreline evolution, and glacier retreat, helping to fill the gap between images, predict future changes, and inform strategies for managing the impacts of climate change in Arctic regions.