Assessing Drivers of Coastal Tundra Retreat at Point Hope, Alaska

Shoreline erosion over the last few decades along the Point Hope, Alaska, USA coastline has led to considerable loss of tundra behind the beach. Analysis of available remote sensing, morphology, and hydrodynamic datasets are used to inform the time scales and mechanisms of coastal land loss at Point Hope, in part through integration of these data into an analytical retreat model originally developed for coastal foredunes. At shorter time scales, these analyses indicate that steeper local beach slopes influence higher wave runup that can enhance the scale of tundra retreat at a particular section of coast. Coastal stretches that are most vulnerable to wave attack of the tundra scarp appear to shift with time related to complex spatio‐temporal variability in shoreline change rates. These local beach erosion rates primarily control the retreat of the tundra at decadal scales. However, independent of the fronting beach morphology, the largest erosion events throughout the region generally occur when there is the coincidence of both high wave energy and elevated still water levels. Additional exploratory modeling results indicate that future changes to sea level and changes in wave energy will further enhance the magnitude of tundra erosion.

In-situ measurements of environmental and morphological data are relatively sparse along the Arctic Coast (e.g., Forbes, 2011;Irrgang et al., 2022). Development of long-term and/or continuous datasets related to beach processes is hampered by the remote nature of the region, seasonally frozen conditions, limited communications networks for data relay, and a lack of established datums. Local to global remote sensing datasets from aerial photography, satellite imagery, or drone photogrammetry are thus commonly used in these analyses (e.g., Cunliffe et al., 2019;Gibbs et al., 2019), with shoreline (or bluff) retreat rates most often derived from visual proxies derived from instantaneous waterlines, wrack lines, human-made features, or local geomorphologic characteristics (Boak & Turner, 2005). Land-based topographic measurements or air or space-based lidar can similarly be used to quantify contour or volumetric changes. Tools such as the Digital Shoreline Analysis System (DSAS; e.g., Himmelstoss et al., 2018) have previously been used to estimate decadal and longer shoreline change rates for much of the US Arctic coastline from various shoreline data products; these analyses collectively indicate that coastal land losses in the region are generally widespread but exhibit considerable alongshore variability (e.g., Gibbs et al., 2019Gibbs et al., , 2021. However, only a select number of studies have focused on event scale evolution of beach, bluff, and tundra systems. Many of the most detailed studies on the oceanographic, geomorphic, and thermal interactions that have driven substantial coastal retreat at are at Drew Point, AK Frederick et al., 2021;Jones et al., 2018;Overeem et al., 2011;Ravens et al., 2012). At this Drew Point site, and numerous others along the northern Alaskan coastline, the adjacent beach is relatively narrow (e.g.,  and the dynamics of coastal tundra retreat in these ice-rich/fine grained sediment settings (e.g., Ping et al., 2011;Ravens & Peterson, 2018) are largely driven by the formation of erosional niches that ultimately fail and slump off. An analysis by Jones et al. (2018) showed that bluff erosion rates at Drew Point over a 9 year period were not found to be strongly correlated with annual variability in relevant environmental variables. This suggests episodic erosion via extreme storm events plays an important contribution to drastic geomorphic responses in the region, with bluff retreat up to ∼19 m/yr having been observed . Additionally, thaw slumping at annual scale is also recognized as an important contributor to the retreat of these features (e.g., Lantuit & Pollard, 2008). However, factors controlling tundra retreat may differ from site to site, with a detailed analysis at Barter Island, Alaska by Gibbs et al. (2021) indicating that annual scale variability in incident wave power is a predominant factor in bluff retreat.
Given the present and forthcoming risk to coastal Alaskan and other Arctic coastlines, a range of numerical tools have been developed to simulate relevant erosive processes. For example, XBeach, which is a process-based nearshore morphodynamic model originally developed for sandy coasts (Roelvink et al., 2009), was modified to include a heat transfer module to represent the role of reduced erodibility of frozen permafrost (Ravens et al., 2017). Other tools have focused on more detailed 3D thermo-mechanics of permafrost bluffs to enable the direct simulation of niche formations and block failure (e.g., Frederick et al., 2021). Whereas these computationally demanding modeling tools enable the incorporation of a broad range of relevant mechanical and thermal processes, the data input requirements, computational demands, and expertise required to successfully run these tools, which may limit their application for engineering-level guidance on hazard mitigation. In part due to these constraints, more empirically or geometrically based tools exist to also simulate beach and/or bluff dynamics relevant to portions of the Arctic coastline. For example, Ravens et al. (2012) developed a simpler modeling framework to predict the development of erosional niches and subsequent block failure using limited local inputs such as local wind speed and beach elevation. Another quasi-process model developed by  included a force-balance to assess bluff stability and was able to account for bluff failure and topple mechanisms. In conjunction with field datasets, their work indicated that bluff erosion was highly episodic and driven in large part by periods with elevated wave energy during the ice-free season.
A primary focus on the geomorphic evolution of coastal environments within the Arctic has been on systems with limited exposed beach and where nearly vertical and/or high elevation ice-rich tundra bluffs directly abut the ocean. While this morphologic classification encompasses portions of the Alaskan Arctic coast, there is in fact a broad range of coastal geomorphic landforms in this region that include coastal barriers and deltas (e.g., Gibbs & Richmond, 2015;Gibbs et al., 2019). For example, Point Hope consists of a cuspate foreland system that juts out into the Chukchi Sea. The shoreline along the northern section of Point Hope has been eroding at approximately 1-2 m/yr over the last few decades , similarly leading to loss of the adjacent tundra that has contributed to the relocation of the town in the 1970s and has threatened culturally significant sites such as ice cellars and sod houses that cannot be easily moved. Whereas permafrost bluffs can be upwards of 10 m vertically In this study, we focus on a 3.5 km stretch of coast along the northern coast of Point Hope to decouple the relative roles of environmental forcing and beach evolution on driving alongshore variable tundra retreat. A combination of geomorphic data from satellite imagery and lidar are used to quantify spatio-temporal patterns in beach and tundra retreat. These changes are put in context with local measurements and global model hindcasts of wave and water level properties. Given that the beach and tundra are independent geomorphic features, thus precluding the extrapolation of shoreline change rates directly to vulnerability of the tundra, morphologic modeling is additionally completed to understand historical and predict future potential landscape losses. Based on the unique geographic and geologic properties of the site relative to Drew Point and other more northern Arctic sites, we take the approach of modifying an analytical model developed for coastal foredune erosion along sandy outer coast systems to adapt relevant processes specific for simulating the erosional dynamics of the Point Hope coastline. Specifically, the analytical model of Palmsten and Holman (2012) (henceforth PH12) is (a) implemented as a cross-shore profile model and adapted to account for shoreline retreat and hindcast erosion at the site, (b) used to investigate these relative morphologic and environmental drivers on tundra scarp retreat, and (c) applied to understand future risk to the ice cellars, archeological resources, and other infrastructure at Point Hope. With the use of this modeling tool we are explicitly hypothesizing that collision of the tundra by high total water levels (TWLs) is the primary driver of backshore tundra erosion at Point Hope. This modeling effort is used in conjunction with available field and environmental datasets to quantify the alongshore drivers and time periods of tundra retreat at Point Hope and to test this hypothesis.
In the manuscript, additional details of the field site are further described in Section 2. Numerous relevant field datasets are presented in Section 3 to characterize the morphology and environmental characteristics of the region. The development of the numerical model and modeling results are described in Section 4. A discussion and conclusions are provided in Sections 5 and 6, respectively.

Study Site
Thousands of years ago, Iñupiat people established a village on the Lisburne Peninsula on the northwestern Alaskan Chukchi Sea coast, near to where the modern-day City of Point Hope resides ( Figure 1). Point Hope is a micro-tidal site with a coastline comprised of modestly sloping beaches (∼0.05-0.1 m/m) adjacent to vegetated tundra that rise approximately 3-5 m in elevation above the local mean sea level. The beach is made of heterogeneous coarse-grained materials, primarily shingle, with coarse gravel and sand (Kindle, 1909). Whereas the beach is frozen during the winter, the beaches thaw to a depth of ∼1.5-3 m deep in the ice-free months which can allow for marine erosion of beach sediments-(North Slope Borough, 2017). Sea ice cover within the Chukchi Sea limits wave energy to reach the shoreline for portions of the year, particularly late spring and early summer (Sturtevant et al., 2004). The average open water period typically spans from around June to November, though models project continued lengthening of the open-water season in the Chukchi Sea which could amplify the wave energy that contributes to coastal erosion and sediment transport in the decades to come Jones et al., 2018;Serreze et al., 2016). Future sea level rise (SLR) projections for the region are uncertain, although analysis from long-term tidal measurements about 400 km south from Point Hope at Nome, AK indicate still water levels have been increasing at 3.89 mm/yr between 1992 and 2021 (Sweet et al., 2022). Over 100 years this rate would contribute to an increase in the mean water level by 0.4 m. Currently air temperatures at Point Hope range from −1 to 10°C in the summer and −17 to −23°C in winter and about 10-12 inches of precipitation, including 36 inches of snow, occurs annually at the field site (Brubaker et al., 2010). Both air temperature and total precipitation are also projected to increase into the future at the field site (Brubaker et al., 2010).
A 2021 erosion exposure assessment found that the Point Hope airport is the only engineered infrastructure on the peninsula that is likely to be directly vulnerable to erosion by 2079, though realignment of the airport runway is already planned to reduce risks to transportation infrastructure on the remote peninsula (Buzard et al., 2021). However, numerous culturally significant sites are currently in the zone of active erosion, such as sod houses and ice cellars (siġluaqs) located along the highest parts of peninsula, immediately adjacent to the ocean on the north shore ( Figure 1). Erosion rates of 0.3-2.4 m/yr (Overbeck et al., 2020) threaten this infrastructure. In attempts to slow erosion along the north shore, the City of Point Hope has established makeshift coastal defenses in recent years, using materials such as armored rocks, sand filled drums, bulldozing of sediment, and large sandbags that have been placed in front of the airport and the siġluaqs (North Slope Borough, 2017). 10.1029/2022JF006813 5 of 26

Background
The beach landscape is continually changing in response to wind, wave, and tidal forces that transport sediment along the coast, offshore to the shoreface, or landward via aeolian processes or overwash. Whereas these morphologic changes to the shoreline can occur on a daily basis (e.g., Phillips et al., 2019), changes to any landscape features landward of the beach are usually more episodic in nature and driven by energetic conditions during storms. Impacts to dune or bluff systems are therefore often found to be controlled by the (TWL; e.g., Ruggiero et al., 2001;Stockdon et al., 2007), where TWL is defined as: where SWL is the local still water level and R u is the vertical excursion of wave runup on the beach. Here the SWL reflects a time-averaged water surface that includes tidal effects and any non-tidal residuals, such as storm surge, as would be measured at a tide gauge. Wave runup is a dynamic process that includes components from both wave setup and swash (e.g., Holman & Sallenger, 1985), both of which are typically empirically predicted. Numerous formulations exist to estimate wave runup, but among the most common approaches is that of Stockdon et al. (2006), which defines the 2% runup exceedance level ( 2% ) as: Where is the foreshore slope, is the offshore significant wave height (also denoted for non-deepwater waves), and is the offshore wavelength.

Still Water Levels
To calculate TWLs at the field site, a continuous record of SWL is necessary. In part to develop baseline environmental characteristics at the study site, a field campaign was carried out in Summer and Fall 2021. An RBR SoloD pressure sensor was deployed on a t-post in the low energy Marryat Inlet ( Figure 1) to characterize local SWL fluctuations. The pressure sensor was deployed from 24 July 2021 to 7 September 2021 and measured continuously at 16Hz. The raw pressure time series was corrected with atmospheric pressure measurements from public data available from the local airport. Data were then averaged in 6 min intervals to characterize the time-variable SWL, as shown in Figure 2a. Oscillations driven by lunar tide at the ∼12 hr periodicity within the lagoon are measurable but less than about a decimeter. However, the SWL is not static and varied considerably over the multiple month period of data availability, with a range of 0.77 m. While these data give insights into the SWL dynamics local to Point Hope, their duration limits extrapolation of hazards to longer time periods.
The closest tidal water level gauge for characterizing the SWLs over the long-term (>years) is located approximately 100 km to the south-east of Point Hope at Red Dog Dock (RDD), AK (NOAA Station 9491094). Nearly continuous SWL measurements are available from this site since 2004, with no data from the site prior to 2003 and intermittent data dropouts after sensor installation. The SWL time series from RDD and Marryat Inlet show moderate correlation (R = 0.565; Figures 2a and 2b) and qualitatively show agreement in patterns at weekly scale, although can differ particularly at hourly intervals. The qualitative coherence between these two locations indicates regional scale forcing patterns dominate the SWL signal at daily to weekly scales that are consistent across the broader region. Root mean square (RMS) and maximum differences between the two sites are 0.16 and 0.64 m, respectively. These differences are normally distributed ( Figure 2c) and appear to be related in part to a muted tidal response within Marryat Inlet relative to RDD, temporal lags in tidal and non-tidal residual processes due to spatial separation between sites, and controls driven by local basin characteristics at the locations of the SWL measurements. In general, the SWL is more dynamic at RDD, with both higher and lower SWLs relative to the measured data at Marryat Inlet. Given that the lagoon itself is sheltered, increased SWL dynamics observed at RDD may also be more representative of the SWL fluctuations along the open coast of Point Hope. Based on the general similarity in SWL signatures between the two available datasets, the RDD data set is used to characterize longer term fluctuations in the SWL at the field site for this analysis. For any of these temporal gaps in the RDD record, the outputs from the global HYCOM model (Chassignet et al., 2006) are filled in for the time series from a node close to the study site. HYCOM shows similar agreement with Marryat Inlet dynamics (R = 0.453), having a similar RMS difference of 0.15 m and maximum difference of 0.54 m for the three-hourly data that is available from the model (Figures 2a, 2e and 2f). Note that SWL datasets from both RDD and HYCOM represent the water level under sea ice and therefore include measurements during ice covered periods.
From the combination of RDD and HYCOM a single, interpolated hourly SWL time series was generated for the time periods from 2004 to 2019, as shown in Figure 3, for further analysis. Note that for consistency with the morphology data, all water level data was converted from native tidal datums into the NAVD88 orthometric datum using the Alaska Tidal Datum Calculator (https://dggs.alaska.gov/hazards/coastal/ak-tidal-datum-portal. html). There were still some vertical datum offsets following conversion, specifically related to HYCOM matching up to the measured data. To avoid systematic vertical biases in the data products, all datasets were modified to have the same time-averaged vertical mean value as the RDD data for further comparison.

Wave Properties
Limited in-situ wave records exist for the Chukchi Sea, in part resulting from sensor risk due to seasonal ice cover. Short term (days to months) wave measurements in the general region have been completed to assess specific scientific questions (e.g., Hošeková et al., 2020) or to support general validation of regional scale wave models (e.g., Pingree-Shippee et al., 2016). During the same time period as the RBR pressure sensor deployment, a Spot- ter wave buoy was deployed offshore of the northern Point Hope coastline in 2021. Bulk wave statistics were recorded on half hourly intervals from the sensor from the period from 23 July 2021 to 12 November 2021. Some anomalous peak wave periods (T p ; e.g., >10 s) were noted for a few hours within the Spotter time series that were removed from record and assumed to be related to strong wind, current effects, and/or landing birds. No additional filtering of these data was completed. These data were then linearly interpolated to fixed hourly intervals to be consistent with other available wave data. As shown in Figure 4, the largest recorded H s in this interval was 3.4 m and the average H s was less than 1 m.
Similarly to the approach of the SWL time series, the measured wave data can be used to inform the suitability of other longer term environmental products for the Point Hope field site. There are many regional or global hydrodynamic models that do include outputs of bulk wave parameters in high spatial detail within the Chukchi Sea region. These include forecasting and hindcasting tools such as WaveWatch III (Tolman, 2009), ERA5 (Hersbach et al., 2020), and the Wave Information Studies (Bryant et al., 2016) that each has different grid/output resolution and parameterizations of relevant physical processes. For this application, based on the public availability of long-term hindcast data, model outputs from the ERA5 model were obtained from the closest output grid node relative to the northern shoreline of Point Hope. A comparison of these short-term measurements to the ERA5 hindcast is provided in Figure 4. When wave directions (D) were oriented from the northwest (270° < D < 360°) according to ERA5, the mean error (RMS) between the field measurements and modeled wave heights and wave periods were 0.12 m (0.23 m) and 0.2 s (1.0 s), respectively. There is about a 15° mean discrepancy in the measured and modeled D from this quadrant (Figures 4h and 4j). Conversely, when considering all wave directions mean differences in H o , peak wave period (T p ), and D were 0.39 m, 0.1 s, and 52°, respectively. These larger errors are primarily related to ERA5 predicting wave occurrence from the southeast, whereas the wave buoy measurements from the northern coast of the spit ( Figure 1a) shows that the Point Hope peninsula shelters waves from this direction. Overall, both wave height and wave period from directions with large open-water fetch show considerable agreement ( Figure 4). This justifies the general use of ERA5 modeled waves for Point Hope for characterizing longer-term controls on beach/tundra erosion, as long as southerly waves are excluded from the analysis. Time series of , T p , and D from the period from 2004 to 2019, corresponding to the time period of available morphology data, are used for longer term analysis ( Figure 3).

Total Water Levels
By replacing R 2% , as defined by Equation 2, for R u in Equation 1, the TWL can be calculated for the field site when the foreshore beach slope is known. A time series of TWL using a representative of 0.1 m/m is shown in Figure 3e. Assuming this constant beach slope, collision of runup with the tundra would be expected 419 hr total in the period of interest between 2004 and 2019, or 0.9% of the time. Conversely, with a of 0.06 m/m the frequency of collision goes down by about a factor of 4 (94 hr, 0.2% of the time) due to slope effects on wave runup (e.g., Stockdon et al., 2006). The TWL concept is used as a driver of tundra retreat in subsequent model development, to be described in Section 4.  Stockdon et al. (2006), with TWL shown in panel e. Wave collision occurs with the tundra when the total water levels exceed the red line, which represents the tundra toe elevation, in panel e. Gray dots represent all data points for waves with northerly directions (270° < D < 90°). Blue dots indicate southerly wave directions that were excluded from the analysis based on the local shoreline orientation.

Annual to Decadal Scale
In this study we utilize all available imagery products that meet quality control criteria to quantify alongshore variable coastal retreat at Point Hope. Gibbs et al. (2019) derived a number of shorelines from an existing NOAA topographic sheet (t-sheet) and aerial photography. Overbeck et al. (2020) derived an additional more recent shoreline from aerial photography. These shoreline vectors used a combination of high water line (HWL) and instantaneous land-water interface (LWI) shoreline proxies, which are compatible to use together in assessing shoreline change at the decadal scale. Two historical archive satellite images were accessed from the Global Enhanced GEOINT Delivery (G-EGD) platform to provide additional recent shorelines, and two satellite images were tasked in fall 2021 to establish the most recent shoreline and tundra positions. All of the shorelines used are cataloged in Table 1 and shown in Figure 5. The LWI was determined to be an acceptable shoreline proxy due to the micro-tidal nature of the site and was manually digitized in ESRI's ArcMap for each G-EGD and tasked satellite images at a scale of 1:1000. Cross-shore oriented transects were cast with 50 m alongshore spacing for a 3.5 km stretch of the Point Hope peninsula (Figure 1). The manually derived shoreline polylines were intersected with these cross-shore transects to define the shoreline position at fixed alongshore increments for each time period.   For the present analysis, long-term SCRs were calculated via linear regression through all of the shoreline data points at each transect using DSAS for the time period between 1950 and 2021. In addition to a DSAS-derived SCR, applicable interval end-point shoreline change rates (EPR) were calculated for shorter time periods from the imagery-derived shorelines to better assess the spatio-temporal complexity of shoreline dynamics within the region of interest. These EPRs are calculated at each DSAS transect location as the SCR between two individual shorelines divided by the time between those shoreline collections. No tidal stage corrections were implemented for any inter-annual SCR calculations given that the timing of all such satellite passes are not known (Table 1) and environmental time series to define the SWL and TWL conditions are not available for the pre-2004 shoreline periods. Thus, these outputs are meant to primarily inform general spatial patterns in shoreline behavior.
The When considering these sub-decadal recent time periods the alongshore variability appears to operate in part on ∼500 m wavelength intervals, as apparent for the 2012 to 2016 and 2016 to 2019 time periods in Figure 6. Additional manual digitization of the vegetation line for images where multispectral data were available, as shown in Figure 5, indicates that there is considerable alongshore variability in the landward retreat of the vegetation line between 2016 and 2021.

Storm Time Scale
Additional satellite imagery was tasked through G-EGD to understand the magnitude of morphologic changes during high energy conditions. Based on tracking of the local weather patterns and the use of forecast waves and water levels from WaveWatch III, an energetic storm in late August 2021 was identified as an event that may have caused tundra collision. Imagery was tasked prior to the event on 19 August 2021 and after the event on 1 September 2021.
Using the tasked imagery and the local environmental products, this specific storm was investigated in detail. Shorelines derived from these images are shown in Figure 5. In between these overflights, the waves peaked on 26 August 2021, with H s reaching up to 3.50 m with an associated T p of 8.5 s according to the wave buoy. The ERA5 model similarly predicted a storm during this period, but only with a maximum H s of 2.89 m and peak wave period of 8.3 s. Utilizing the ERA5 and RDD time series (consistent with the model implementation), the instantaneous TWLs at the time of the 19 August 2021 and 1 September 2021 shorelines were 1.3 and 1.6 m, respectively, given an average local beach slope of 0.076 m/m derived from available 2018 airborne lidar data (to be described). Correcting for these TWL differences at the imagery collection time indicates a net average (maximum) estimated loss of −1.9 m (−4.9 m) of shoreline. Given that the satellite imagery has horizontal uncertainties of 3-6 m and additional uncertainty induced by the digitization procedure (e.g., Gibbs et al., 2019), this scale of potential shoreline changes at the storm time scale measured by satellite are less than or equivalent to the uncertainties of those measurements. Nearly no net change in the manually derived vegetation lines from 19 August 2021 and 1 September 2021 further indicates that the tundra was not severely impacted during this event.
Unfortunately, there are no other clear satellite images available at the storm timescale from other time periods for the region to further constrain tundra impacts at this shorter timescale.

Profile Change
Airborne lidar was acquired over a subset of Point Hope in 2004, 2018, and 2019 that provides insights on the coastal profile change of both the beach and tundra. The geographic coverage of each of these overflights varied, although none of them covered the entire Point Hope peninsula. These ground-filtered data were downloaded natively as.las files from NOAA Digital Coast (https://coast.noaa.gov/dataviewer) and linearly interpolated onto cross-shore transects with cross-shore grid resolution of 0.1 m every 50 m in the alongshore at the same transects as the shoreline locations (Section 3.2.1). Examples of coastal profile evolution over this timeframe are shown in Figure 5 for transects located every 250 m in the alongshore (see Figure 1 for transect locations). Note that the only sites where useable lidar data of the entire beach-tundra system was available across all three surveys is at 9 cross-shore transects, all located in the vicinity of the ice cellars between transects 50 and 58 (see Figures 1  and 5). The morphologic changes between 2004 and 2018 are analyzed in depth at these nine transect locations. A larger number of profiles, spanning from near the spit tip (Transect 13) to east of ice cellar locations (Transect 58) had data in both 2018 and 2019. This annual time scale is analyzed separately.
Between 2004 and 2019 there has been considerable retreat (up to 40 m) of the beach at shoreline. For the purposes of this work we define the shoreline as 1.5 m NAVD based on data availability at this contour elevation across the available transects and datasets, although this elevation is slightly higher than the local mean high water elevation (MHW; 1.12 m NAVD). . The highest topographic points, where ice cellars and sod houses had formerly existed in Figures 7d and 7e, for example, have been lost to the sea in this timeframe. SCR were calculated from the lidar data as the time evolution of the 1.5 m NAVD contour. For the purposes of this work the tundra was defined as being any location landward of the 3.5 m NAVD contour, with the tundra change rate (TCR) calculated using this contour location. This elevation for the toe of the tundra is based on the approximate regionwide values of maximum curvature on the upper beach, consistent with how the dune toe elevation is commonly defined (e.g., Stockdon et al., 2007). The beach slope was calculated as the slope between the 1.5 and 3.5 m NAVD contours. The data show that the steepest beach sections ( > 0.1 m/m) have observed TCRs that exceed the SCR (Figure 8). Generally, the locations with the highest shoreline change rates (SCR < −2 m/yr) have TCR that are lower than the SCR, indicating that the SCR is not the sole factor driving the retreat of tundra. There has been an average and maximum of 50.1 m 3 /m and 61.8 m 3 /m of net volumetric loss of the tundra above the 3.5 m NAVD contour, as measured between 2004 and 2018 and considering landscape changes up to 5 m landward of the 2018 3.5 m contour.
Whereas at decadal scale (2004-2018) the tundra has retreated landward close to the ice cellars, the lidar data at annual scale between 2018 and 2019 indicates a more muted response. Over this time scale, SCR varied between −5 m/yr and +2 m/yr and TCR varied between −1.6 m/yr and 0.6 m/yr. Although it is of note that morphologic changes to the beach/tundra profile at this time scale are also likely to include some effects related to ice push, placement of sand bags for coastal protection, anthropogenic movement of materials, scarp face thaw or slumping, and/or vertical and horizontal lidar position uncertainty.

Sediment Size
Surface sediment samples were collected on the beach face and the tundra scarp on the Point Hope spit. Figure 9 shows the grain size distributions from two individual samples that are assumed to be representative to the site. The beach sample was composed nearly entirely of gravels and had a median grain size (D 50 ) of 14.9 mm. The tundra sample was bimodal with about 22% of the sample consisting of silts and clays, 70.8% sands, and 7% gravel. The D 50 for the tundra sample was 1.4 mm.

Model Formulations
While the available observational data indicate that there has been widespread coastal retreat, there remain gaps in terms of understanding the timescales and drivers of tundra erosion. In the absence of additional datasets, a numerical modeling analysis is completed to supplement the available data to inform the relative oceanographic and morphologic controls on tundra retreat at the study site. The PH12 model builds off the wave impact theory model of Larson et al. (2004) for coastal foredunes and utilizes the assumption that the volume losses of the dune (∆ ) are linearly related to the momentum flux impacting the dune face, according to: whereby is a model coefficient related to the erodibility of sediment, t is the duration of exposure, toe is the vertical elevation of the tundra scarp toe (see Figure 10), and T is the wave period. Consistent with other applications and with Equations 1 and 2, the TWL is the sum of the SWL and the R 2% calculated using Stockdon  (2006). Despite less process complexity relative to other coastal retreat models (e.g., XBeach), the PH12 model has shown skill at simulating spatio-temporal dune retreat rates in both laboratory and field settings. (Cohn et al., 2021;Palmsten & Holman, 2012;Splinter & Palmsten, 2012). The primary tunable model coefficient is which has reported values of between 0.0025 and 0.07 for sandy dunes.
The numerical implementation of PH12 assumes that a vertical scarp is developed during erosion events and that scarp is not subject to avalanching. The original PH12 model assumed that the landward trajectory of toe fell along a linear slope ( ) that was proportional to . For the purposes of this implementation at Point Hope, it is assumed instead that toe is constant-in part based on the low elevation of the tundra which limits and application of the tool beyond the storm time scale. PH12 also by default assumes a runup factor of 1.26 whereby the outputs of Equation 2 are increased by 26% before the TWL computation, reflecting the fact that runup can be higher on the backshore against dune/scarp features due to local slope effects relative to that predicted close to the shoreline (e.g., the lower portion of the beach). This assumption is kept for the purposes of this model implementation. The additional assumption about the persistence of small scarps is also maintained, consistent with observations from the field (Figure 1f).
In the absence of detailed thermal data for Point Hope, for the purposes of the present work we neglect potential thermal constraints on the erodibility of sediment which are known to be important in permafrost settings for scarp collapse and erodibility (e.g., Frederick et al., 2021). However, it is of note that at this site there is no wave energy that reaches the shoreline when the sea is iced over and thus thermal effects are indirectly accounted for when the ocean water temperature is close to or below freezing. Additionally, given the low-lying nature of the tundra landscape at Point Hope, the vertical distance from toe to the top of the tundra is typically less than 2 m, resulting in a small potential scarp and differing from the morphodynamics of coastal permafrost bluffs on the Beaufort Sea coast.
The PH12 model has previously been utilized for simulating dune erosion at the storm time scale and have not included any parameterizations of morphological beach changes. However, the shoreline is evolving at multiple meters per year in some locations at Point Hope. In this present application, a constant SCR is applied at the shoreline (X shore ) to account for the aggregate effects of longshore transport gradients and subaqueous cross-shore processes at the shoreline in order to extend the timescales of useability of the analytical model. Between X shore and the cross-shore location of the scarp toe (X toe ), a linear slope is then calculated to define for input to Equation 2. In the absence of more complicated couplings, the SCR is treated as chronic and continual erosion that drives a landward migration of the shore at hourly time scales. Given that shore remains at a fixed elevation, this results in a general steeping of f during time periods where the scarp toe remains fixed. However, toe also migrates landward during erosional storm events, which subsequently reduces f. This time variability in feeds back into the potential for scarp retreat through feedbacks on R 2% and freeboard above the scarptoe. There is no pre-condition to the model on the TCR; instead X toe naturally evolves as a function of the oceanographic forcing and feedbacks of wave runup from time-variable changes to .

Model Hindcasts
Two sets of model hindcasts were run at the sites of interest. First a set of cases are run from 2004 to 2018, which are initialized with the lidar topography from 2004, for the 9 transects where useable data are available. The  model is forced with wave properties from the nearest ERA5 node, as shown in Figures 3a and 3b and shoaled to 25 m water depth using linear wave theory. During time periods where the ERA5 mean wave direction was from the south (90° < D < 270°), H s was assumed to be 0 m. Wave runup is calculated according to Equation 2 for each hourly time step of the model and is added to the SWL time series to define the TWL (Equation 1). For the initial analysis the lidar-derived SCR at 1.5 m NAVD for each interval was input into the model to define the constant retrogradation of X shore . toe is assumed to be constant at 3.5 m NAVD. is recalculated each hour to define R 2% based on the length scale between X shore and X toe . A model default C s value of 0.011, as determined by Cohn et al. (2021) for modeling of a coastal dune system in response to a hurricane and Nor'Easter event, was assumed. No model calibration or tuning of the C s parameter was completed. Outputs from these model simulations are used to constrain the time periods and mechanisms of tundra retreat in close proximity to the ice cellars.
To further confirm the ability of the model to simulate coastal behavior across time and space scales, additional hindcast simulations were completed for the time period from 2018 to 2019. The model was forced with the same environmental input sources as the previous cases, although the broader geographic coverage of lidar data for these two time periods permitted cases to be run for about 40 transects. The SCR from this time period was similarly derived from the lidar data for self-consistency.

Model Exploration
Many of the 2004 lidar profiles showed a topographic high associated with ice cellar locations that was close to the shoreline. It is unclear if details of the high antecedent morphology in 2004 contributed specifically to those sites being erosional hotspots along the coastline. To gain some insights into these feedback effects, an additional set of simulations was completed that initialized the model with the 2019 measured topography, but utilized the 2004 to 2018 time series of environmental forcings. These outputs are compared to the initial set of decadal scale simulations to comment on how future erosional behavior may differ from past erosional behavior given the present shape of the beach-tundra landscape. These cases are also run with modified SCRs, sea level rise rates, and changes in wave energy to understand potential drivers of future change at the site.

Deterministic Model Hindcasts
Simulations were completed independently for the periods from 2004 to 2018 and 2018 to 2019. Examples of the time series of TWLs, beach slope, and tundra retreat metrics, as simulated according to the model, are shown in Figure 11. For illustrative purposes we focus on outputs related to representative Transects 51 and 58 to describe the overall timing of erosion at the site. For the decadal scale simulations, the model indicates that at example transect 58 there were 282 hr between 2004 and 2018 (122,713 hr total in the record) where the tundra was in the collision regime. Due to slope and morphology variations, this duration of predicted wave impact was longshore variable. Example transect 51 was estimated to have been in the collision regime 134 hr for the same time period. These hours of potential erosion were  often consecutive during individual storms, although numerous episodic erosion events did occur in the record that led to a net tundra sediment loss of up to 30 m 3 /m.
In cross-shore profile view, the model predicted the loss of a local topographic high-the site of former ice cellars-along the majority of coastline where data were available for both 2004 and 2018 (Figures 12a-12i). There is visibly an under-estimation of the erosion at some transects, with the former topographic high feature not entirely eroded by the model (e.g., Figure 12f). Conversely, there is over-prediction of horizontal tundra retreat and volume losses in other select transects (Figures 12a and 13). However, overall the model generally compared favorably to the data. The model over-predicted erosion at the 3.5 m contour on average (root mean square error, RMSE) by 2.0 m (3.0 m) (Figure 12a). The model does appropriately predict that the TCR is lower than the SCR for the zones of fastest shoreline retreat and that the TCR is faster than the SCR for zones of the slowest shoreline retreat (Figures 8 and 13). Noteably, this consistency between the field datasets and model results highlights that, with equivalent environmental forcing, the local SCR is a primary factor controlling tundra response. Volumetric comparisons are also similar between the model and field, with a mean underprediction and RMSE of 1.5 m 3 /m and 6.1 m 3 /m. Note that the model does only account for wave-driven erosional processes, whereas the measurements do indicate some aggradation of the tundra surface that is likely to be due to ice effects or anthropogenic activities (e.g., Figures 12a, 12c and 12d) that influence these volumetric error statistics. In the case of the annual simulations, both the field measurements and numerical model indicate that there was nearly no retreat of the tundra ( Figure 13) with an average error (RMSE) of <0.1 m (0.1 m) between the two estimates. The model simulated similar behavior between all transects. Both Transects 51 and 58 had similar initial beach slopes, but the higher SCR at Transect 58 ultimately resulted in 69% and 30% more volumetric tundra erosion and horizontal tundra retreat, respectively, relative to that simulated at Transect 51. This again highlights the importance of local beach properties, especially SCR, on vulnerability of the tundra to collisional impacts.
The average conditions that led to erosion for the Transect 51 over the decadal time scale was a H s of 3.5 m, T p of 9 s, and SWLs of 1.55 m; these conditions being well above the climatological average for the region (e.g., Figure 3). Figures 14b-14d shows the probability distributions of wave and water level climatology and those conditions that drove erosion, indicating that the confluence of these factors is the typical driver of tundra retreat at the site. These conditions do not occur in the shorefast ice period given that no or minimal wave energy reaches the shoreline during this time period. The conditions that can cause erosion may happen as early as July, according to the model, and has the highest probability of occurrence in November (Figure 15).

Exploratory Model Simulations
Decadal scale simulations in which each model was initialized with 2019 measured topography instead of that from 2004 indicated that all portions of the coast would still be expected to be erosional at decadal scale, independent of the initial morphology. However, there was an average 27% reduction in volumetric tundra loss using the updated profiles-suggesting that the present day morphology (e.g., reduced tundra height and differing beach slopes) is generally less prone to erosion than the morphologic configuration of previous decades. However, there were some simulated locations where volume losses increased (by up to 10.4 m 3 /m) in these new simulations as well, indicating that the shape of the landscape plays an influential role in the potential for impacts.
SCRs also control the magnitude of volumetric tundra retreat and corresponding landward retreat of the beach-tundra boundary. The model predicts that in the case of a static shoreline (0 m/yr SCR) that there would still have been a loss of −5.6 m 3 /m of sediment from the tundra between 2004 and 2018, although these rates increase to −16.2 m 3 /m, −34.3 m 3 /m and −41.4 m 3 /m for SCR of −1 m/yr, −2 m/yr, and −3 m/yr, respectively ( Figure 16). Uniform 10%, 20%, and 30% increases in H s are expected to drive increases of the −1 m/yr erosion condition to −18.9 m 3 /m, −21.5 m 3 /m, and −24.1 m 3 /m, respectively (Figure 16a). NOAA does not report a SLR rate for the RDD station due to the relatively limited record, although a SLR rate of 3.89 mm/yr is reported from Nome, AK (Sweet et al., 2022). The SWL time series here between 2004 and 2018 suggests a similar rate is likely appropriate for RDD and Point Hope, with a linear regression through the SWL series indicating a rate from this short time period of around 6.3 mm/yr. Further increases in sea level were also added to the model to explore potential feedbacks on the tundra retreat. Additional increases in sea level increasing the frequency of TWL collision with the tundra (Figure 16b). Therefore, there is a trend of increasing tundra erosion with increasing SLR for fixed SCR rates. If an additional 3.89 mm/yr of SLR was added, approximately doubling the current rate of SWL increases, the tundra volume losses are expected to increase by 14.4%, 10.7%, 1.3%, and 0.9% for SCR of 0 m/yr, −1 m/yr, −2 m/yr, and −3 m/ yr relative to the baseline case. If SLR rates were increased to 7.78 mm/yr then these volume losses would be expected to increase to 31.2%, 22.4%, 2.7%, and 1.8% relative to the no additional SLR scenario. These results suggests that SLR will have a larger influence on increasing volumetric erosion along relatively static shorelines relative to shorelines that are already rapidly retreating, which is largely due to the increased relative number of hours of tundra collision with and without SLR for low SCR cases. This observation highlights the critical role of SCR on determining the scale of tundra retreat, independent of future changes in SLR or environmental forcing.

Geomorphic Evolution of Point Hope
Satellite-derived SCR from available imagery and chart products indicate that erosion rates of over −1 m/yr have been ongoing since at least the 1950s along the northern coast of Point Hope (Figures 5 and 6). The shoreline data does not show any clear trends in SCR acceleration since the 1950s. However, imagery available at shorter (sub-decadal) time scales indicates that these changes are not temporally or spatially constant ( Figure 5). The magnitude of shoreline change was often largest in proximity to the spit tip (Figure 6a). To the east of the spit tip, there were also regular rhythmic shoreline oscillations with a wavelength of ∼500 m during the period from 2012 to 2021 (Figure 6b). These oscillatory shoreline and shoreline change patterns are consistent with the presence of megacusp embayments that are characteristic of many coastal settings and which may contribute to enhanced erosion of backshore dune and bluff features (Castelle et al., 2015;Thornton et al., 2007). Whether zones of high shoreline change and possible cusp formation are associated with bathymetric templating (e.g., Castelle et al., 2015;Cohn et al., 2021), offshore geological controls (e.g., McNinch, 2004), heterogeneity of bed sediments (e.g., Guest & Hay, 2019), or self-organized processes (e.g., Coco et al., 2000) at Point Hope is not clear from the available data. The coherent downcoast migration of zones of progradation and erosion (Figure 6b) could be related to the propagation of erosional hotspots (e.g., Ashton et al., 2003), indicating that the locations most at risk along the spit are likely to vary with time. This in part explains why short-term (<5 years) SCR can vary substantially from the long-term SCR rates derived from aerial imagery products. Shoreline retreat is also likely at the storm timescale, although even moderate storm events do not demonstrate clear erosion patterns beyond the scale of uncertainty from remote sensing (Sections 3.2 and 4.3.1). The relative roles of cross-shore beach erosion and longshore processes in these shoreline dynamics, across timescales of storms to decades, is as of yet unknown at this site without additional analysis or modeling. Although given that there is a coupling between the shoreline dynamics and magnitudes of tundra retreat (Figure 8), the loss of high ground at the site is intrinsically linked to these complex beach dynamics in some capacity. Whereas longshore transport processes are likely to be a major factor in beach erosion at the site, high TWLs that drive cross-shore redistribution of sediments are the primary driver of retreat of tundra scarps at the storm time-scale. The available field data and observations suggest the presence of only modest scarps (Figure 1f), the height of which are largely limited by the low-lying nature of the local topography. The minimal height of the scarps limits the potential for block failure mechanisms of permafrost at this site. It is also thought that multiple meters of thaw can occur on the beach at the site (North Slope Borough, 2017), likely distinguishing the dynamics of the peripheral Arctic circle landscapes, such as Point Hope, from those at more northern latitudes that have a longer ice coverage duration and differing geomorphic and geologic framework. This poses limitations for the use of tools specifically developed to simulate these block failure mechanisms and their contributions to land loss in these permafrost laden coastal landscapes (e.g., Frederick et al., 2021;Jones et al., 2018;Overeem et al., 2011) at Point Hope. The available geomorphic observations at Point Hope also show some aggradation of the tundra surface close to the beach-tundra interface (Figures 7a-7c, 7e-7g and 7i) which is likely to be a combination of ice push, coastal management by the community, or wind-driven growth effects. While these effects do influence volume change calculations of the tundra losses and do complicate numerical prediction of the local landform based on TWL and SCR concepts alone, it is of note that these modifications are small relative to the net sediment losses of the coastal landform on the beach system at decadal scale (e.g., Figure 7). Upwards of 60 m 3 /m of sediment was removed from some cross-shore transects at Point Hope between 2004 and 2018, with this large erosion including the loss of many former ice cellars (Figure 7). Whereas the ice cellars are purposely built into regions of high silt and clay content, the available grain size data indicates that the current tundra scarp face in regions where these cellars are not located may contain substantial quantities of sands and gravels (∼80% volumetrically). Similarly to bluff and scarp systems from other geographic locations, it is possible that the erosion of the tundra may contribute substantially at short time scales to beach sediment supply with these coarser clasts (e.g., Runyan & Griggs, 2003). The lack of silts and clays within the collected beach sample (Figure 9) may be an indication that the finest grain sizes are transported offshore and not retained within the active littoral system. While the beach sample also had minimal sand size material (Figure 9), aerial imagery (e.g., Figure 1d) and field observations suggest that there is wide heterogeneity in beach composition throughout the area of interest, some of which is likely to include coarse sand size material that appears to be characteristic of the tundra composition ( Figure 9). This wide variability in beach grain size also has important implications for the shape of the beach, in part due to well documented relationships between and D 50 (e.g., McFall, 2019). These complex linkages between these grain size properties, alongshore variable SCR and , and tundra vulnerability likely play a critical role in when, where, and why there have been zones of enhanced erosion during certain time intervals in the historical record.
Locations with faster SCRs generally observed large landward migration of the beach-tundra interface between 2004 and 2018 (Figures 7 and 8). However, the ratio of the TCR and SCR is below 1 for the fastest SCR (SCR < −2 m/yr). Although the tundra did rapidly retreat in these locations, the lower TCR could in part relate to these zones having among the shallowest , which translates to lower wave runup relative to steeper coastal sections, as per Equation 2. Conversely, the TCR exceeds SCR where exceeds 0.1 m/m (in 2004) and where SCR > −1.5 m/yr. These trends may also be related to time-evolving SCR, in part associated with longshore transport gradients and potential megacusp migration.

Environmental Drivers of Coastal Change
Appropriate characterization of the environmental forcings is critical to understanding the environmental drivers of erosion. Although daily, astronomically-driven tidal ranges of about 10 cm at Point Hope, there was upwards of 60 cm of SWL fluctuation observed in Marryat Inlet over the July-September period associated with regional scale meteorological and oceanographic forcing (Figure 2). The longer available water level record from RDD showed a 3.9 m SWL range between 2004 and 2019. This range has important implications on the potential for tundra impacts and therefore SWL processes are necessary to be incorporated into local hazard assessments despite this being a micro-tidal system. While the RDD data shows close agreement with the SWL dynamics at Marryat Inlet and appears to adequately characterize the timing and magnitude of SWL fluctuations at the field site in the absence of local long-term measurements, there was similar agreement with the global HYCOM model. Future work could constrain the drivers of variability in SWL at this regional scale from these spatial model outputs.
The SWL often makes up a major component of the TWL (Figures 3d and 3e), but waves remain the primary driver of coastal change through the generation of spatial gradients in sediment transport (Figure 14). In the summer and fall 2021 field season, storms with H s exceeding 3 m were observed ( Figure 4a) and wave periods were generally between 3 and 9 s throughout this time period. When the dominant wave direction was from the north, there is considerable similarity in the locally measured waves and those from the ERA5 model hindcast.
Although the unique geometry of the field site, whereby the Point Hope peninsula sticks out over 10 km into the Chukchi Sea, is not resolved on the relatively coarse (0.25°) global model wave grid. Thus, the closest output grid node does include waves coming from the south that in reality are sheltered by the peninsula. There were additionally some biases for waves from the northeast, likely also associated with poor characterization of the land-water boundary with coarse grid cells. Although in general, the <0.25 m RMSE for 270° < D < 360° and close agreement in measured and modeled wave periods gives confidence in these hindcast products for reasonable definition of wave energy in proximity to Point Hope. When considering waves only from the north, the average significant wave height according to ERA5 was 1.3 m, with 1668 hr and 20 hr of waves exceeding H s of 3 and 5 m, respectively, between 2014 and 2019. There was strong seasonality in these wave conditions, with ice limiting any waves to reach the study site in most years in the months of December to April and wave energy generally increasing from May to October/November (Figure 3).
The presented framework assumes that erosional impacts to the tundra are dependent on TWLs. Since the model reasonably accurately simulates historical tundra retreat, this suggests that high TWLs are a dominant process in this setting, as was hypothesized (Section 1). The model suggests that it is typically the confluence of energetic wave and SWL properties that ultimately contribute to tundra retreat (Figure 14), the combination of factors which are relatively rare in the environmental record (Figures 3 and 14). This is in part what appears to drive relatively mild tundra impacts on an annual basis, with large erosional events that intermittently driving large scale tundra retreat ( Figure 11). These findings are in line with those of Jones et al. (2018) that suggests that mean environmental forcings in the Arctic are not the driver of tundra retreat and that instead the extreme storm events are likely most at cause. The confluence of conditions during the ice-free season makes the fall, specifically November at Point Hope according to this analysis, the period with the most number of hours of potential tundra collision.
Given that there were observed differences in the RDD and ERA5 products relative to local site measurements and the unique geometry of the Point Hope site, the generated long-term oceanographic record likely does not characterize the details of all extreme storms correctly along the northern shore of Point Hope (e.g., Sections 3.1.2 and 3.1.3). This includes potential biases from coarse resolution of the shoreline and potential inadequate definition of shorefast ice and their relative effects on waves in ERA5 (e.g., Hošeková et al., 2021). Thus, while the data indicates the relative rigor of these products, it is important to note that the timing and magnitude of predicted erosional impacts is dependent on these datasets for which the ultimate uncertainty on their use at Point Hope is unknown beyond the summer and fall 2021 data collection period. However, the possibility to utilize modeled environmental data-with some understanding of the uncertainty of those outputs-does enable a more complete understanding of present and future risks to Point Hope data that the available sparse local in-situ data alone cannot provide.

Predictive Skill of Hindcasting Coastal Change at Point Hope
The analytical model of PH12 has been successfully used in numerous locations around the world to simulate sandy coastal foredune erosion at time scales of hours to months (Cohn et al., 2021;Splinter & Palmsten, 2012;Splinter et al., 2018) with the assumption of a static shoreline. In this present application to Point Hope, the longer time scales of interest have necessitated added treatment of a moving shoreline. This additionally allows for a time-evolving beach slope that leads to additional feedbacks that can enhance or decrease vulnerability for backshore erosion related to wave runup effects. Slope effects contributed to numerous years of coastal retreat predicted by the model from 2010 to 2014 (Figure 11b), followed by numerous years thereafter associated with decreased risk following beach slope shallowing. Although limited morphologic data exists to confirm these trends on the storm-to inter-annual time scale, the available lidar data does indicate a net shallowing of the beach profile in the vicinity of the ice cellars from 2004 to 2018-consistent with these numerical predictions.
Despite widely different geographic and geologic characteristics of coastal foredunes and the intermittently frozen tundra of Point Hope, a default C s value used for coastal foredunes was found to be appropriate for this work. This likely suggests the general similarity in the erodibility of the sediments within the tundra relative to dune materials, which are similarly usually composed of sands and have vegetation present. Additional sensitivity analysis testing of C s indicated that this parameter does strongly influence the magnitude of simulated net tundra volume erosion and retreat (not shown). The model results are instead generally much more sensitive to the oceanographic conditions and SCR input to the model relative to changes in the C s coefficient. As the offshore oceanographic conditions are similar across the northern shoreline of Point Hope, this would suggest that spatial variability in SCR has a predominant control on the local magnitude of tundra retreat.
Although mostly empirical devices are used to simulate the timing and magnitude of erosion, the analytical model implemented here shows some skill across multiple time scales and alongshore locations, with minimal required model inputs or tuneable parameters. In its present implementation, there is a lack of incorporation of thermal properties which are generally important to consider for most Arctic sites (e.g., Frederick et al., 2021). The present application may not benefit from additional complexity in these thermal model parameterizations given that the largest magnitudes of erosion are expected to occur in the fall coincident with the period of largest thaw of seasonally frozen soils. However, these dynamics could change into the future with changes in ocean and atmospheric temperatures, storm patterns, and ice coverage. Additional complexity for thermal dynamics on landscape erodibility could be added through time and space modification of the C s parameter.
A further limitation to the current analytical approach in this modeling framework is a pre-definition of the tundra scarp toe. The definition of this toe elevation is consistent with needs for coastal dune applications (e.g., Splinter & Palmsten, 2012), although can be poorly defined in sites that do not have sharp changes in curvature between the beach and the hinterland feature. The maximum curvature based definition of the toe at this site from the available data also appears justifiable, although it is important to note that there are sensitivities to the model based on this parameter as there are with similar dune applications (e.g., Wernette et al., 2018). Other more sophisticated process-based modeling tools do not inherently require this definition and can also directly simulate changes to the beach and tundra landscape (e.g., Arctic XBeach) which may be necessary to resolve for certain detail oriented applications, such as the design of engineering structures. However, these more complex tools also often have limited reliable predictive skill beyond the storm time scale due to compounding feedback effects (e.g., Cohn et al., 2021;Vousdoukas et al., 2011). Significant model tuning and the addition of morphologic stabilization routines that limit runaway erosion (e.g., Roelvink et al., 2019) may enable the scalability of these simulations for annual or longer scales. Continually improving computational power may also increase the computational scalability/accessibility of some of these tools. While more complex process-based numerical tools do exist, the adapted PH12 model application to this site appears to be appropriate for the specific time scales and mechanisms for this site based on the model skill across time scales and transects.

Implications for Infrastructural Impacts and Future Site Changes
The model results generally indicate that both the antecedent morphology and environmental characteristics each play a role in controlling the timing and magnitude of coastal tundra retreat. The presence of topographic highs close to the beach, such as were present in 2004 where the ice cellar mounds were present, may serve to reinforce higher than average beach slopes by maintaining the tundra toe in a more fixed position and that in turn drives above average erosion. While the shape of the landscape undoubtedly plays a role in the style of erosion that occurs, SCRs dominate the scale of tundra retreat. SCR have historically been both spatially and temporally variable at Point Hope ( Figure 5). It is unclear how SCR will change in the future, although it is likely to be influenced by increases in wave energy due to reduced sea ice coverage and by further modifications to sediment supply (Mahoney et al., 2014). The model results suggest that wave energy increases alone are sufficient to increase TCR, even in the absence of future changes in SCR. The possibility of the co-occurrence of increased SCR, increased wave energy, and sea level change is likely to lead to compounding effects that drastically increase present day volumetric losses to the tundra system ( Figure 16). This puts the remaining ice cellars at Point Hope at continued short-term risk of loss of these culturally significant and immoveable infrastructural features in the absence of preventative measures (e.g., hard or soft infrastructure solutions). The only engineered infrastructure at direct risk from erosion on the northern portion of Point Hope is the airport runway. Realignment of the Point Hope Airport to reestablish the required buffer area around the runway was estimated to total $20-30 million, as of 2015, and is expected to be constructed in the next few years (Alaska DOT, 2015). This is a significant investment for a community with a population of only about 800 people, with these infrastructure improvement/re-location needs only likely to increase throughout the Arctic in the context of rapidly changing coastal conditions. In contrast to existing plans for airport realignment, culturally significant locations, such as the Old Tikiġaq site and the Ipiutak archeological site, which is a U.S. National Historic Landmark located northeast of the airport, will require alternative means of protection if they are to not be lost to the ocean (Dawson et al., 2017). Beyond Point Hope, many coastal communities across Alaska experiencing the impacts of erosion face relocation or the staggering costs of engineered solutions to adapt (Mason et al., 2012).

Conclusions
Coastal systems throughout the world are in a period of rapid change. Understanding the drivers of these changes and potential changes into the future, especially in rapidly changing Arctic systems, is critical for defining present and future risk in order to develop suitable adaptation strategies. In the case of Point Hope, AK, chronic erosion of the landscape poses major risk to infrastructure and historic sites, such as the community ice cellars. While shoreline erosion is continual at inter-annual scale, unique field datasets and numerical modeling suggest that the tundra retreat itself is episodic in nature and driven by the confluence of high SWLs and/or wave energy that lead to collisional impacts. These impacts have historically been most probable in late fall before the onset of shorefast ice, but as the ice-free season lengthens it can be expected that more storms may lead to wave impact on the tundra scarp. While this has the likelihood of continuing to enhance erosion, there remains a strong coupling between the beach morphology evolution and the response of the tundra to climate-driven changes in permafrost and wave action. Thus future changes to the peninsula must inherently account for these morphologic changes to the beach which will likewise be driven by changes in both wave magnitude and direction. Model simulations incorporating variable SCRs suggest that indeed these erosional trends will be exacerbated by forthcoming changes in sea level and wave energy.