Arctic coastal hazard assessment considering permafrost thaw subsidence, coastal erosion, and flooding

The thawing of permafrost in the Arctic has led to an increase in coastal land loss, flooding, and ground subsidence, seriously threatening civil infrastructure and coastal communities. However, a lack of tools for synthetic hazard assessment of the Arctic coast has hindered effective response measures. We developed a holistic framework, the Arctic Coastal Hazard Index (ACHI), to assess the vulnerability of Arctic coasts to permafrost thawing, coastal erosion, and coastal flooding. We quantified the coastal permafrost thaw potential (PTP) through regional assessment of thaw subsidence using ground settlement index. The calculations of the ground settlement index involve utilizing projections of permafrost conditions, including future regional mean annual ground temperature, active layer thickness, and talik thickness. The predicted thaw subsidence was validated through a comparison with observed long-term subsidence data. The ACHI incorporates the PTP into seven physical and ecological variables for coastal hazard assessment: shoreline type, habitat, relief, wind exposure, wave exposure, surge potential, and sea-level rise. The coastal hazard assessment was conducted for each 1 km2 coastline of North Slope Borough, Alaska in the 2060s under the Representative Concentration Pathway 4.5 and 8.5 forcing scenarios. The areas that are prone to coastal hazards were identified by mapping the distribution pattern of the ACHI. The calculated coastal hazards potential was subjected to validation by comparing it with the observed and historical long-term coastal erosion mean rates. This framework for Arctic coastal assessment may assist policy and decision-making for adaptation, mitigation strategies, and civil infrastructure planning.


Introduction
The Arctic is experiencing significant changes and warming up to four times faster than the rate at lower latitudes (Rantanen et al 2022). An increase in air temperature leads to permafrost thaw and ground ice melt, which can lead to talik formation, thermokarst development, and associated ponding. Taliks are defined as perennially unfrozen zones above or within the permafrost (Ferrians et al 1969). In coastal areas, permafrost thaw, talik formation, and thermokarst processes can contribute to coastal erosion, as thawed sediments are more easily eroded than frozen ones.
Ocean wave dynamics also play an important role in coastal erosion in the Arctic. During calm conditions, waves can thermally erode frozen bluffs, whereas extreme storm surges and surface ocean waves mechanically increase the vulnerability of the Arctic coast to erosion, leading to flooding and enhanced ground thawing in flooded areas (Irrgang et al 2022). The decline of sea-ice extent and duration increases the exposure to storms, enhancing Arctic wave-and storm-driven erosion. Projection suggests that sea ice will continue to decline as stormy conditions overlap more frequently with the open water sea (Barnhart et al 2016, Crawford et al 2021. Consequently, Arctic coastal hazards are expected to increase over the coming decades. Predicting the natural hazards of Arctic coasts for civil infrastructure planning and understanding society's capacity to adapt and transform are crucial for effectively preparing for an uncertain future. Several assessment methods have been developed to evaluate and visualize the hazard potential associated with permafrost thaw subsidence, coastal erosion, and flooding. In thaw subsidence assessment, settlement (Nelson et al 2001) and bearing capacity (Streletskiy et al 2012) indices are the two primary approaches. Other methods related to subsidence assessment include the flow-diagrambased risk zonation index (Daanen et al 2011) and the analytical hierarchy process-based index (Hong et al 2014). For the evaluation of coastal erosion and flooding, Gornitz (1994) developed a Coastal Vulnerability Index that considers six variables: geomorphology, coastal slope, rate of relative sea-level rise, shoreline erosion, mean tide range, and mean wave height. Jaskólski et al (2018) applied this Coastal Vulnerability Index to assess high Arctic coastal vulnerability specifically in Longyearbyen, Svalbard. Arkema et al (2013) calculated a coastal hazard index (HI) for Alaska coastlines, aiming to identify areas at risk of inundation and erosion. The coastal HI integrates seven variables: habitat characteristics, shoreline type, coastal relief, sea-level rise, wind, wave, and surge potential. Manson et al (2019) employed a coastal sensitivity index to evaluate the sensitivity of Canada's marine coasts to climate change impacts. Casas-Prat and Wang (2020) computed pan-Arctic coastal hazards, focusing on wave characteristics and using the Inundation Index and Erosion Index (Jiménez et al 2012). The Inundation Index and Erosion Index comprise three variables: wave height, peak wave period, and storm duration.
To date, the thaw subsidence indices do not include zones of talik formation where previously frozen ground is now thawed year-round. This exclusion might lead to an underestimation of thaw subsidence in the analysis (Wang et al 2023). The development of assessment methods for evaluating combined Arctic coastal hazards is relatively rare when compared to their widespread utilization in low-latitude coastal regions. Previous efforts have been made by developing or applying indicator-based methods (Gornitz et al 1994, Ford andSmith 2004) or an empirical model (Casas-Prat andWang et al 2020, Nielsen et al 2022) to assess Arctic coastal hazard or vulnerability. However, permafrost thaw was not explicitly considered in these frameworks. No attempt was made to integrate permafrost thaw, coastal erosion, and coastal flooding into a synthetic hazard assessment of the Arctic coast. To fill this knowledge gap, this study conducted an integrated analysis of the coastal plains of Alaska's northernmost borough, North Slope Borough.
The objectives of this research are to (1) present a regional-scale permafrost thaw subsidence assessment for North Slope Borough, Alaska, (2) develop a holistic framework for Arctic coastal hazards assessment that considers permafrost thaw subsidence, coastal erosion, and flooding, and (3) identify coastal areas of North Slope Borough that are prone to the three types of hazards.

Study area
This research studies the coastal regions of the North Slope Borough in Alaska. Figure 1 depicts the location and four coastal communities. They are Utqiaġvik, Wainwright, Point Lay, and Kaktovik. Continuous permafrost (90%-100% of areas underlain by perennially frozen ground) underlies the North Slope Borough with a thickness ranging from 200 m to more than 600 m (Osterkamp and Payne 1981). The permafrost has moderate to high ground ice content, with some parts of the coastal lowlands reaching up to approximately 80% ground ice by volume (Kanevskiy et al 2013). Owing to the presence of ice-rich permafrost and low-lying tundra terrain, the landscapes in these villages are vulnerable to widespread thaw subsidence, coastal erosion, and coastal flooding.

Data
In our assessment of thaw subsidence hazard, we utilized the pre-existing ground ice map of northern Alaska developed by Jorgenson et al (2014). This map provides information on the distribution of massive ice volume within the upper five meters of permafrost. For reference, the ground ice map can be found in figure S1 in the supplementary materials of this study.
The input variables that were utilized for the assessment of coastal hazard, including habitat and shoreline type, were derived from the ShoreZone database (NOAA ShoreZone). Specific details regarding the assignment of exposure ranks for different shoreline types and habitat classes can be found in supplementary tables S1 and S2, respectively. The spatial distribution of shore-type exposure rankings is visualized in figure S2, while the visualized habitat classes are depicted in supplementary figure S3. To determine the variable of coastal relief, we utilized an Arctic digital elevation model obtained from the Polar Geospatial Center (Porter et al 2018). For the exposure factors related to wind, wave, and surge potential, we utilized a compiled dataset derived from the WAVEWATCH III model spanning eight years (2008-2016) (NCP 2018).
To determine the sea level projections, we compiled data of mean projections of sea-level rise along the U.S. coastlines from the National Oceanic and Atmospheric Administration, as presented by Sweet et al (2022). These sea-level rise scenarios are related to the air temperature scenarios in this study, estimated based on the fifth phase of Coupled Model Intercomparison Project (CMICP5) models (Sweet et al 2022). For more detailed information regarding the classification and representation of the sea-level rise scenarios, please refer to figure S4.

General methodological framework
We integrate the index-based permafrost thaw potential (PTP) as an input variable into the Arctic coastal hazard assessment framework, as shown in figure 2. The conceptual framework consists of three steps: (1) assessment of permafrost conditions; (2) assessment of permafrost thaw subsidence; (3) assessment of Arctic coastal hazards.
First, in-situ field observations of mean annual ground temperature (MAGT) and active layer thickness (ALT), ground thermal properties, and environmental data serve as inputs for the Geophysical Institute's Permafrost Laboratory (GIPL)-2.0 model. MAGT, ALT and talik thickness projections from the GIPL model for the Representative Concentration Pathway (RCP) 4.5 and 8.4 are visualized in Geographical Information System (GIS)-based maps for the 2020s and the 2060s.
Second, using the input from thaw subsidence factors derived from existing ground ice maps (Jorgenson et al 2014) and the predicted MAGT, ALT, and talik thickness from the GIPL-2.0 model, the ground settlement index is quantified. The hazard indices are then used to delineate areas into five categories of PTP caused by the thawing of near-surface permafrost.
Third, the erosion and flooding factors combined with the gridded PTP from the second step constitute the input factors of the Arctic Coastal Hazard Index (ACHI). In this study, coastal hazard refers to thaw subsidence, flooding, and erosion caused by storms and permafrost thaw acting upon shorelines. Environmental risk encompasses the product of hazard occurrence and its potential societal consequences. ACHI is determined by incorporating thaw subsidence and erosion-flooding potential. In the following section, we describe each methodological step in detail.

Assessment of permafrost conditions
The GIPL-2.0 model simulates soil temperature dynamics and the depth of seasonal freezing and thawing by numerically solving 1D quasi-linear heat conductive equations with phase change . The heat equation used in the GIPL model requires spatially distributed thermal conductivity, volumetric heat capacity, volumetric unfrozen water content, and volumetric latent heat of freezing and thawing as the input parameters (figure 2).
The permafrost conditions of North Slope Borough, Alaska were simulated using the GIPL-2.0 model by the co-authors (Nicolsky et al 2017). The simulation results are visualized on the website: https://permamap.gi.alaska.edu/. This study uses existing simulation data for hazard assessment.

Assessment of permafrost thaw subsidence
We considered the role of taliks in the potential of permafrost thaw subsidence. Taliks in this study are defined by their physical state as unfrozen zones above or within permafrost. As mean annual air temperature increases, the depth of summer thaw begins to exceed the depth of winter freeze, resulting in the formation of sub-aerial taliks, as shown in figure 3. Farquharson et al (2022) used the concepts of potential freeze and potential thaw to determine the point in time when a sub-aerial supra-permafrost talik may begin to develop. The depth of potential freeze is characterized as the maximum extent of freeze reached by the end of a freezing period, assuming complete thawing of the subsurface material at the onset of freezing (Alexiades 1992). The depth of potential thaw is determined as the maximum depth of thaw attained by the conclusion of the thawing period, assuming complete subsurface freezing at the initiation of thawing (Alexiades 1992). Nelson et al (2001) developed an approach to estimate the potential thaw subsidence of permafrost. The ground settlement index is constructed using data of changes in ALT and ground ice content. The calculation is based on the following assumptions: the liquid water produced by the thawing of ground ice drains away from the affected sites and thaw settlement is proportional to the thickness of the ice lost (Nelson et al 2002). The calculated thaw settlement signifies an upper bound estimate, as the computation assumes a transformation of the thawed permafrost into a non-porous state. In permafrost regions with sub-aerial taliks, talik thickness plays a crucial role in the occurrence of potential thaw settlement. Liquid water resulting from the talik formation drains away from the impacted sites. Solely considering changes in the ALT or the seasonally thawed/frozen layer may underestimate the maximum potential thaw settlement. To address this limitation, we propose an improved approach that incorporates changes in talik thickness to accurately estimate ground subsidence. The ground settlement index (I s ) is computed using equation (1): where ∆Z ALT + TT = ∆(Z ALT + Z TT ) is the relative changes in the combined ALT and talik thickness (m), as shown in figure 3, and V ice is the volumetric ground ice content. In order to achieve a detailed delineation of hazard conditions, we categorized the distribution of I s for the 2060s in comparison to the 2020s using manual intervals. These intervals were determined based on the theoretical potential range of absolute thaw settlement (measured in meters). The five ranges for potential thaw settlement were defined as follows: 0.0-0.05 m, 0.05-0.15 m, 0.15-0.25 m, 0.25-0.35 m, and greater than 0.35 m. The upper bound of each range is inclusive. It is worth noting that the class breaks for both the RCP4.5 and RCP8.5 scenarios remained consistent to facilitate a comparison of temporal and spatial variations in hazard conditions.

Assessment of Arctic coastal hazards
We calculate the values of Arctic coastal hazards with future climate change scenarios through the coastal where R Habitats refers to the habitat type, R ShorelineType represents the geomorphology type, R Relief represents the coastal relief (i.e. topography), R Wind is the wind exposure, R Waves is the wave exposure, R SurgePotential is the surge potential, and R SLR represents sea-level rise scenarios. Because each of the seven factors has a range of 1-5 and HI is an average value of these factors. HI has a range of 1-5. A detailed explanation of each class is provided in the supplementary materials. The Arctic coastal setting is influenced by complex interaction mechanisms between thermal and mechanical drivers. The variability in Arctic Coastal Dynamics (ACD) results from coastal settings and environmental drivers. Arctic coastal settings are determined by the wave energy, coastal morphology, lithological characteristics (Wang et al 2022a(Wang et al , 2022b, and the ground thermal regime. Environmental drivers include air and water temperature, sea ice, wave climatology, storm intensity and timing, and sea-level changes (Irrgang et al 2022). When ice-rich permafrost bluffs are subjected to thawing due to MAGT increases, the volumetric land loss may be three times higher than that directly associated with coastal erosion (Lim et al 2020). We developed a composite ACHI based on equation (2), by including the permafrost thaw subsidence potential, as shown in equation (3): where R PTP represents the permafrost thaw potential (PTP). R PTP is derived from index-based permafrost thaw hazard maps from step 2 with a spatial resolution of 1 km. The complete coastline was partitioned into 10 376 segments, with a spatial resolution of 1 km × 1 km for each individual segment. Raster maps of the potential thaw settlement are used to calculate the PTP of each 1 km 2 shoreline area (i.e. the segment point). The five-class thaw subsidence hazard maps are first given a ranked number of 1, 2, 3, 4, or 5 representing low to high hazard potential. We calculate the average PTP within a radius of 3 km to decrease the uncertainty. The average PTP is then assigned to the nearest segment point to represent the shoreline PTP. The gridded shoreline PTP can be found in supplementary figure S7.

Results and discussion
3.1. Regional-scale projection of permafrost thaw subsidence potential Two GIS-based thaw subsidence hazard maps were generated in this study, comparing the 2060s to the 2020s under two climate forcing scenarios, RCP4.5 and RCP8.5 (figure 4). The two hazard maps exhibit relatively consistent geographic variation. Regions with relatively higher potential for thaw subsidence are primarily concentrated along the west coast of the Beaufort Sea, the middle portion of the North Slope Borough, and low-latitude upland areas.
Under the RCP4.5 scenario, a significant portion of the Alaska coastal plain exhibits potential thaw settlement ranges between 0 and 5 cm from the 2020s to the 2060s. Some sections of the west Beaufort Sea coast show potential thaw subsidence ranges between 5 and 15 cm during the same time frame. On the other hand, under the RCP8.5 scenario, major parts of the Alaska coastal plain display potential thaw settlement ranges between 5 and 15 cm from the 2020s to the 2060s. Additionally, certain areas along the west Beaufort Sea coast exhibit potential thaw subsidence greater than 35 cm during this period.
The spatial patterns observed for thaw subsidence in our study align with the findings of a previous assessment conducted by Jorgenson et al (2014) regarding potential thaw settlement. Jorgenson et al (2014) developed a settlement map that primarily relied on the characterization of ground ice. We propose that the spatial patterns for potential thaw subsidence are primarily influenced by the presence and distribution of ground ice, which, in turn, is closely associated with cryolithology and Quaternary depositional environments.

Validation of permafrost thaw subsidence
In order to evaluate the efficacy of the ground settlement index in estimating thaw subsidence, a comparative analysis was conducted between our calculated estimates for potential thaw settlement and the observed surface subsidence in regions of Northern Alaska where permafrost is prevalent. Figure 5(a) shows the prediction of potential thaw settlement along the Beaufort Sea coast using the ground settlement index. The calculated potential thaw settlement is for the 2060s compared to the 2020s under RCP8.5 forcing scenarios. The study region that is delineated by the dashed black line is situated in Prudhoe Bay region, positioned at the central portion of the Beaufort Sea coast. Figure 5(b) depicts the long-term subsidence rates in the study region, spanning from 1992 to 2000, examined using Interferometric Synthetic Aperture Radar (InSAR) measurements. Figure 5(b) shows observed longterm trends in surface subsidence at rates of 1-4 cm per decade, as reported by Liu et al (2010). In figure 5(c), we have rescaled the predictions of regional settlement in centimeters per decade. According to figure 5(c), over 95% of the study region (calculated based on pixel numbers) exhibits a potential settlement ranging from 1 to 4 cm per decade. The projected values and distribution for thaw subsidence rates between 2020s and 2060s align with the longterm trends observed between 1992 and 2000.

Coastal hazards evaluation for North Slope Borough, Alaska
The distributions of the ACHI are classified into quartiles to represent different levels of coastal hazard scores. The quartiles are defined as follows: high (upper 25%, ranging between 3.5 and 5), moderate (50%-75%, ranging between 3.2 and 3.5), low (25%-50%, ranging between 2.8 and 3.2), and stable (lower 25%, ranging from 1 to 2.8); the upper bound of each range is inclusive. Figure 6 shows the predicted ACHI results of coastal hazard potential of North Slope Borough, Alaska by 2060s using RCP4.5 and 8.5 climate forcing. Both hazard maps demonstrate a relatively consistent geographic distribution of ACHI. Regions with relatively higher potential for coastal hazards are primarily concentrated along the west coast of the Beaufort Sea, including the coastline of Point Barrow, Smith Bay, and Pogik Bay. Additionally, the regions between Wainwright and Peard Bay also exhibit elevated levels of potential coastal hazards.
In the RCP4.5 scenario, the coastal areas spanning from Point Barrow to Pogik Bay primarily exhibit a coastal hazard potential ranging from low to moderate levels until the 2060s. However, under the RCP8.5 scenario, a significant portion of the same coastal region is characterized by a high potential for coastal hazards. Additionally, specific areas along the Wainwright, Point Lay, and Kaktovik coastline also demonstrate a notable high potential for coastal hazards, primarily due to high ground ice contents.
Regions characterized by a relatively higher potential for coastal hazards are driven primarily by cryolithology and Quaternary geology. Specifically, the higher potential for coastal hazards observed between Smith Bay and Pogik Bay can be attributed to several factors. Notably, the increased wind, wave, and surge potential in this region, along with the exposed shoreline type characterized by ice-rich silt, contribute significantly to the heightened hazard potential. Moreover, when considering permafrost thawing, it becomes evident that areas along the west Chukchi Sea with greater hazard potential are primarily influenced by permafrost thaw and the absence of protective features within their exposed habitat class.

Validation of Arctic coastal hazard potential
To assess the efficacy of ACHI in evaluating Arctic coastal hazard, we conducted a comparative analysis between the estimated coastal hazard potential based on RCP8.5 forcing scenario and the observed long-term coastal erosion mean rates obtained from the ACD database (Lantuit et al 2020) (figure 7(a)). We also considered the historical long-term rates of shoreline change reported by Gibbs et al (2017) (figure 7(b)). As shown in figures 7(a) and (b), the historical data were reclassified into four distinct classes, namely stable or aggrading, slow erosion (0-1 m per year), moderate erosion (1-2 m per year), and rapid erosion (>2 m per year). The projected ACHI results were then categorized into quartiles to represent varying levels of coastal hazard potential: stable (0%-25%), low (25%-50%), moderate (50%-75%), and high (>75%); the upper bound of each range is inclusive.
The spatial distribution of the predicted patterns by ACHI ( figure 7(c)) demonstrated a strong agreement with both the observed and historical long-term coastal hazard data. Within the coastal region from Point Barrow to Pogik Bay, 92% of segments display a notable potential for moderate to high coastal hazards. This finding aligns with the historical data reported by Gibbs et al (2017) and Lantuit et al (2020), revealing nearly 90% and 86%, respectively.
In other regions, such as the east and west Chukchi Sea coasts, as well as the central and eastern Beaufort coast, our projected values also align with the historical data, indicating a propensity for stable to moderate coastal hazards. It is worth noting that the regions between Wainwright and Peard Bay exhibit a high coastal hazard potential, primarily attributed to the increased likelihood of thaw settlement in these areas. The comparisons were drawn across different timelines: the historical data represent shoreline movement or erosion under past conditions, whereas the ACHI is designed to forecast future shoreline conditions.

Limitations of this study
The proposed framework has several limitations: 1) Crylithology is a key factor influencing the ACHI yet ground ice mapping to date has been conducted at a relatively coarse resolution. As such, for many regions, the ground ice volumes are not well constrained. 2) Similarly, vegetation cover, organic layer thickness, and lithology are essential in the ground thermal regime. Still, field-based observations are very sparse, which limits the accuracy of numerical modeling efforts. 3) Decline in Arctic sea ice extent and increase in the open water period accelerate coastal erosion and impact biological productivity (Barnhart et al 2016). However, sea ice extent factors are not considered in the ACHI because of data deficiency. 4) The presence of barrier islands plays a role in reducing wave energy and protecting the coast from erosion (Irrgang et al 2022, Nederhoff et al 2022. However, the computation of wind and wave exposure did not consider the impact of barrier islands, which might lead to an overestimation of the exposure in protected regions. Furthermore, how barrier island morphology and movement may change with an increasingly longer open water season is very poorly understood (Moore et al 2014, Gibbs et al 2023. 5) The validation of predicted results has challenges due to the limited availability of field observations and ground-truthing data. Few long-term monitoring sites are situated along the Arctic coast.

Conclusions
This study developed the ACHI to comprehensively analyze Arctic coastal hazards by integrating thermal and mechanical drivers into a single spatial representation. To develop the ACHI, a large-scale permafrost thaw subsidence assessment is first conducted for North Slope Borough, Alaska. The assessment calculates the potential of maximum thaw settlement and determines the relative hazard ranking of PTP of those regions in 2060s compared to 2020s. The development of taliks is considered in the ground settlement index. The predicted values of thaw settlement demonstrate alignment with long-term observations of thaw subsidence in Northern Alaska. The coastal PTP of the North Slope Borough is then gridded and assigned to the coastline segments as an input factor of the ACHI. The ACHI incorporates PTP with classified geomorphology, habitat, wind, wave, surge potential, coastal relief, and sea-level rise in coastal hazard assessment. The ACHI provides a macro-picture of the North Slope Borough's heterogeneous coastal hazard potential at spatial scale of 1 km 2 under current and future climate change scenarios. The ACHI prediction shows that by the middle of this century, the areas categorized as high hazard are mainly distributed along the west Beaufort Sea coast, including Point Barrow, Smith Bay, and Pogik Bay. The predicted coastal hazard potential aligns with the observed and historical longterm coastal erosion mean rates.

Data availability statements
All data that support the findings of this study are included within the article (and any supplementary files).