Archaeological Site Vulnerability Modelling: The Influence of High Impact Storm Events on Models of Shoreline Erosion in the Western Canadian Arctic

Abstract Much of the Inuvialuit archaeological record is situated along shorelines of the western Canadian Arctic. These coastal sites are at substantial risk of damage due to a number of geomorphological processes at work in the region. The identification of threatened heritage remains is critical in the Mackenzie Delta, where landscape changes are taking place at an increasingly rapid pace. This paper outlines some preliminary observations from a research program directed toward identifying vulnerable archaeological remains within the Inuvialuit Settlement Region. Coastal erosion rates have been calculated for over 280 km of the Kugmallit Bay shoreline, extending along the eastern extent of Richards Island and neighbouring areas of the Tuktoyaktuk Peninsula. Helicopter surveys conducted during the 2014 field season confirmed that areas exposed to heavy erosive forces in the past continue to erode at alarming rates. Some of the calculated rates, however, have proven far too conservative. An extreme period of erosion at Toker Point in the autumn of 2013 has yielded a prime example of how increasingly volatile weather patterns can influence shoreline erosion models. It has also provided a case with which to demonstrate the value of using more recent, shorter time-interval imagery in assessing impacts to cultural landscapes.


Introduction
The archaeological record of the Canadian Arctic extends back over 5000 years (Friesen 2016) consisting in some places of extremely well-preserved material remains. The stable layer of permafrost and dry conditions of the polar desert which foster such high levels of preservation are deteriorating at an increasingly rapid rate however, largely on account of warming climate trends (Anisimov et al. 2007, Blankholm 2009, Lantuit et al. 2012. Shoreline erosion in the western regions of the Canadian Arctic has proven to be a more immediate factor in the loss of cultural materials, due in part to the predominantly marine-focussed lifeways and settlement patterns of people who have lived in the region since time immemorial (Alunik, et al. 2003). This tendency to settle in coastal locations has prompted many northern communities, including

Original Study Open Access
Tuktoyaktuk in the Northwest Territories (Johnson et al. 2003) and Shishmaref in Alaska (Mason et al. 2012), to seek new ways of adapting to the unprecedented rates of change taking place there. Understanding the historical trajectory of shoreline change can yield useful insights on where best to focus efforts toward the management of coastally situated cultural landscapes. The Beaufort Sea coast is characterized as one of the largest areas of high sensitivity shoreline in the circumpolar Arctic (Shaw et al. 1998). As a result, archaeological sites located there are particularly vulnerable to climate related impacts such as permafrost depletion (Smith et al. 2009), backshore thaw slump activation (Lantz, Kokelj 2008), increased severity and frequency of wind events and related storm surges (Pisaric et al. 2011, Small et al. 2011, and shoreline erosion (Solomon 2005). The region is notorious for its powerful northwesterly wind events, which typically increase in strength and frequency in the late autumn prior to freeze-up (Manson, Solomon 2007, Solomon 2005. Studies have shown that storm surges influenced by the severity of fall storms and the duration of open water conditions can have substantial impacts on the rate of shoreline retreat among permafrost and massive-ice laden coastlines (Solomon 2014).
Renewed interest in resource extraction, as well as transoceanic shipping and cruise-ship tourism in the Beaufort Sea area have also been forecast to have a negative impact on coastally situated archaeological sites (AMAP/CAFF/SDWG 2013, Cinq-Mars, Pilon 1991). Given the impacts that changes in prevailing climatic conditions are having on the stability of both modern communities and archaeological sites in the Canadian Arctic, it is imperative that steps be taken quickly to improve our understanding of these impacts and capacities to address them. To this end, I have developed a Geographic Information System (GIS)-based method of identifying vulnerable archaeological remains in the Mackenzie Delta region. This model has been created as part of the Arctic Cultural Heritage at Risk (Arctic CHAR) project, directed by Prof. Max Friesen of the University of Toronto in partnership with the Inuvialuit Cultural Resource Centre.
Over 150 archival air photos and satellite images of the core Arctic CHAR study area have been processed (NAPL 1950, NAPL 1972, MVAP 2007, DigitialGlobe Inc 2013, yielding coastal change rate figures for over 280 km of the Beaufort Sea shoreline. This historical shoreline change model was developed using the ArcMAP software suite versions 10.1 -10.3 (ESRI Inc) and the Digital Shoreline Analysis System (DSAS), which is an ArcMAP extension developed by the US Geological Survey (USGS 2009). The accuracy of the results have been assessed through ground truthing efforts over the course of three field seasons. This process has shown that while the calculated shoreline trends are largely accurate with regard to the modelled versus observed rates of shoreline erosion throughout the core project area, one site in particular has undergone far higher erosion rates than the GIS model predicted. This paper will outline the likely cause of this discrepancy, thereby demonstrating the value of more recent, shorter time-interval imagery in assessing impacts to cultural landscapes.

Methods
Erosion modelling efforts were focussed on the Kugmallit Bay region of the Arctic CHAR project area (Fig.1). This region was selected in part due to the abundance of known archaeological sites in the area, but also due to the availability of digital imagery with which to produce the model. While archaeologists have long used aerial photographs in other parts of the world (Bewley, Raczkowski 2002, Crawford 1929, the earliest available imagery of the Beaufort Sea coastline dates to 1950. The Royal Canadian Air Force, the National Research Council and the Department of Mines and Technical Surveys began collaborating in 1924 toward the goal of improving the accuracy of existing Canadian maps. This concerted effort did not extend to the Arctic region until 1945 however, after which three squadrons of aircraft were dedicated to the task of mapping the whole of the Canadian Arctic (Thomas 1950). Despite this tremendous effort, one of the greatest impediments to undertaking coastal change analysis in the Canadian Arctic is the limited nature of available imagery. Historical air photos were obtained from the Natural Resources Canada National Air Photo Library (NAPL). The University of Toronto Map and Data Library secured 108 air photos in support of this research, which included imagery of the shores of Kugmallit Bay extending from the eastern extent of Richards Island down to the mouth of the East Channel, east to Tuktoyaktuk and north to Toker Point. NAPL imagery was selected for the summers of 1950 (1:40,000 scale) and 1972 (1:60,000 scale), both of which were scanned at 300 dpi by NAPL staff (NAPL 1950, NAPL 1972. These two periods were selected because of the minimal cloud cover and absence of shorefast ice in the imagery, thereby facilitating more accurate assessments of past shoreline positions. More recent imagery was sourced from the Mackenzie Valley Air Photo (MVAP) project (MVAP 2005). Indian and Northern Affairs Canada initiated the MVAP project to support decision making related to the proposed development of the Mackenzie Valley Gas Pipeline Project. The results of this project have included a more accurate topographic model of the region, consisting of contour interval data as well as a digital elevation model. The project has also yielded over 4500 natural colour digital orthophotos at 1:30,000 scale, all of which have been georectified using differential GPS ground survey controls (Schwarz et al 2007). This data has been made freely available to the public through the Northwest Territories Centre for Geomatics data portal (http://www.geomatics.gov.nt.ca/dldsoptions.aspx).
Though not incorporated into the regional erosion model overall, a series of satellite images was recently obtained through a DigitalGlobe Foundation imagery grant. The grant provided over 3000 square kilometres of panchromatic (0.5 m resolution) and multispectral (2.0 m resolution) imagery from the GeoEye-1 satellite. This data extended over the core erosion modelling area and included a combination of imagery from 2009 and 2013 (GeoEye-1 2009, GeoEye-1 2013). As discussed further below, this imagery has proven quite useful for refining erosion rate estimates in the vicinity of archaeological remains.
Each of the air photos and satellite images had to be overlain and aligned with each other within the GIS software, using a process referred to as georectification. This process involved importing each of the NAPL and GeoEye-1 images into ArcMap, where geographic control points were assigned to minimize the degree of distortion between each image and a base-layer of reliably aligned control imagery. The spatial accuracy of the MVAP aerial photography led to its use as control imagery for this project. The type of transformation method used in the georectification process determines how the air photo is twisted, rotated and generally deformed to ensure an accurate fit to the control points. In this case, a spline transformation was used which optimized local accuracy near the control points at the expense of global, whole sheet accuracy (ESRI Inc. 2014).
A minimum of ten control points are required to apply a spline transformation to the rectified imagery, though the number of points applied to each image for present purposes varied from 15 to 35 depending on its level of discrepancy with the MVAP control imagery. In order to best apply the locally optimized spline transformation process, control points were first assigned to the four corners of each air photo, essentially tacking it down onto the MVAP imagery. Further points were then assigned along and near the shoreline zone as needed. The degree of error for each batch of control points was checked using a second-order polynomial transformation, yielding an average root mean square value of less than 3.0 m in every case. The rectified images were then inspected individually to ensure adequate closeness of fit along the MVAP shoreline areas, with further control points added where necessary.
With all of the imagery georeferenced, the shoreline digitization process could begin. The DSAS User Guide noted that: "Shoreline positions can reference several different features such as the vegetation line, the high water line, the low water line, or the wet/dry line.... It is strongly recommended that initial datapreparation steps be taken to reference all shoreline vectors to the same feature (for example, mean high water) before using DSAS to compute change statistics." (Thieler 2009: 10). In order to accommodate the relatively low resolution of the NAPL air photos, the more readily visible 'wet-dry' line was selected as a standardized proxy for overall coastal movement. Shorelines were digitized manually for each time period (Fig. 2) and saved as features in a geodatabase format to facilitate processing in the DSAS ArcMap extension.
It should be noted that the degree of land-water interface visibility varied greatly between time periods covered by the air photos, with the more recent imagery yielding much more confident shoreline positions than the earlier 1950 and 1972 air photos. The highly variable coastal morphology also caused problems in the digitization process. High, shadow casting bluffs and gently sloping, vegetation rich intertidal zones made accurately identifying the wet-dry line all but impossible in some instances. These and other difficulties inherent to erosion modelling in an Arctic setting will be discussed in a forthcoming publication. Suffice it to say however that while there are limitations to this method, the approach does hold merit as a means of producing a generalized model of historical coastal change rates. For examples of similar methods applied in heritage management efforts elsewhere, see also: Radosavljevic et al (2015), Reeder et al. (2012) and Robinson et al. (2010).
The DSAS software was used to generate a series of transects oriented perpendicular to the coast in 10 m intervals. These transects were then plotted over the digitized shoreline positions to calculate change rates from each transect interval of known length and temporal span, as seen in the 'Rates' section of Figure  2. The software can accommodate a range of different statistical methods to describe trends in the shoreline change data. The end point rate (EPR), which accounts for change between the oldest and youngest shoreline positions, is attractive as a straightforward measure of change. Its usefulness in this particular study was limited however due to its inability to account for changes between erosional and depositional shoreline trends through time. The abundant sources of sediment in the region and strong coastal processes with which to transport them have led to the frequent formation and removal of ephemeral coastal features like sand spits and barrier islands throughout the study area (Peletier, Medioli 2014). The linear regression rate (LRR) was selected for its ability to incorporate all three shoreline positional data sets, and thereby its capacity for accommodating changes in shoreline trend (Thieler et al 2009). The LRR was calculated for the entire study area (Fig. 3), spanning the 54 year period covered by the NAPL and MVAP imagery.

Results
The results of the shoreline change analysis can be seen in Figure 4. Linear transects were clipped to the length of the shoreline change envelope and then colour coded according to LRR values in metres per year. Positive LRR values resulted from coastlines with depositional/progradational trends, as can be seen at the tip of Topkak Spit in the inset map. Shorelines dominated by erosional trends yielded negative LRR values, as was the case throughout much of the study area. Shoreline trends from all 25620 transects within the study area yielded an average annual LRR of -0.62 m/a. This erosional trend and change rate for the region as a whole align with observations made by Solomon, who noted that the primary forcing mechanisms behind erosion of the Beaufort Sea Coastline are long term sea level rise and periodic storms (2005).
The differential nature of shoreline change taking place along the Beaufort Sea coast has been noted by a number of researchers (Harper et al. 1985, Solomon 2005. As seen in Figure 4, areas like the mouth of the East Channel which are sheltered from strong northwesterly winds typically exhibit lower erosive rates and a higher frequency of depositional trends. Areas with higher erosive rates tended to have a northwesterly exposure, such as Toker Point, Beluga Point, an area referred to as Imnaqpaaluk or 'Melting Bluffs' by Tuktoyaktuk Elders (Hart 2011), and a 4 km long stretch of shoreline west of Canyanek Inlet. The susceptibility of these areas to storm impacts is compounded by their setting among ice-rich, low elevation regions of fine grained lacustrine sediments (Pelletier, Medioli 2014).
In order to test the accuracy of coastal change rate calculations, an area with known archaeological sites was required which exhibited high rates of erosion. Toker Point, in the far eastern extent of the study area (Fig 5), was selected for this purpose. A line of surveyor stakes was established during the 2013 Arctic CHAR field season. The stakes were set perpendicular to the average strike of the shoreline at 2 m, 10 m, 20 m and 30 m distances to facilitate monitoring over the course of subsequent survey visits. Because of the high energy nature of this exposed coast (Fig. 6), stakes could not be placed at the wet-dry line. This restriction led to measurements being taken from the active edge of the eroding backshore area in the years to follow, under the assumption that rates of change would be roughly equivalent. Following the 2013 field season, shoreline change rate calculations for the Toker Point area were bolstered by the inclusion of a fourth shoreline position. Using the methods previously outlined, the 2009 GeoEye-1 imagery was digitized and a new set of LRR-based shoreline trends were calculated. This revised model yielded average LRR values of -1.46 m/a for the Toker Point area overall, -3.70 m/a for the eastern extent of the region, and -4.27 m/a for the 10 transects closest to the monitoring stakes (Fig. 7). The persistence of this historical erosion rate into recent years can be seen in the Figure 7  Upon returning to Toker Point in 2014, there appeared to have been very little shoreline movement over the course of the 2013 storm season, as the stake nearest the shoreline was still relatively close to the edge. Upon closer inspection however, it became clear that what was originally mistaken for the 2 m stake was in fact the 10 m stake. The shoreline at Toker Point (Fig. 8), had receded roughly 8.5 m over the course of a single year, far exceeding even the highest calculated retreat rates for the area. While some level of discrepancy between modelled and observed rates was anticipated for the first set of monitoring surveys, this magnitude of difference clearly signalled that the model had failed in some way. After returning from the field, the model was reassessed to determine whether or not such a large discrepancy could be attributable to an error in the methods employed.  In light of the susceptibility of the LRR statistic to outlier effects and thereby its tendency to underestimate rates of change compared to other statistical methods (Thieler et al. 2009), the EPR values for all Toker Point transects were calculated. EPR rates were all very similar to their respective LRR values for the Toker Point region overall, as well as the eastern extent of Toker Point and from transects near the monitoring stakes (Table 1). Thus, it would appear that the type of statistic used was not a factor in the underestimated erosion calculations in this particular instance, perhaps on account of the relative abundance of heavily erosive shoreline trends which influenced the extent of regional averages. Another possible factor behind the low modelled rate estimates was the conflation of all three time periods into a single LRR statistic. Calculating EPR values for each time interval individually would allow for the assessment of change through time, possibly yielding insights as to the cause of the low erosion estimate.
To determine whether Toker Point erosion rates had changed over each of the three time periods, EPR values were calculated for the 1950 -1972, 1972 -2004, and 2004 -2009 intervals individually (Fig. 9). These calculations demonstrated that the average rates of erosion acting on the east Toker Point shoreline appear to have decreased over time. It should be noted however that much like the overall study area, erosion rates do vary considerably throughout the Toker Point shoreline zone. It should also be noted that these spatially prescribed trends do show a marked level of variation from time period to time period. While this localised erosion event does not explain the 8.5 m of shoreline loss seen at the monitoring stakes in 2014, it does demonstrate that such losses are possible under the right conditions. Given that such high rates of erosion can occur in the eastern Toker Point region, the possibility remained that discrepancies in modelled versus observed rates may be related to shortcomings in the method when accounting for occasional extreme erosive events. Changes in the spatial extent of erosive impacts are linked to geological, cyrological and geomorphological conditions at any given location, as evident in preferential erosion of the ice rich landform in Figure 10. Changes through time however are most commonly attributed to long-term fluctuations in storminess and the duration of open water conditions (Solomon 2005). Winds exceeding 50 km/h have been correlated with pronounced periods of shoreline erosion within the study area (Manson, Solomon 2007), while shorefast ice which limits the erosive impact of waves typically persists in the Beaufort Sea from October through to June (Solomon 2005). In order to determine if Toker Point experienced an unusually active storm season during the autumn of 2013, wind speed and directional data from the nearest weather monitoring station at Tuktoyaktuk was obtained from Environment Canada's online climate data archives (http://climate.weather.gc.ca/).
The number of days with wind speeds in excess of 50km/h can be seen in Figure 11. Despite the two years with missing data, the region clearly experienced a series of closely spaced high wind episodes in 2013, a trend which continued into the first half of 2014 as well. Of the 123 wind events that took place in 2013, 105 (85%) were out of the northwest. Assuming that shorefast ice is present from October through until June, 28 of the northwesterly high wind events took place during open-water conditions. 20 of these open water-related wind events occurred after the 2013 survey took place, eight of which involved speeds in excess of 60 km/h, and two others exceeded 90 km/h. A further 12 days with northwesterly winds exceeding 50 km/h took place in the open-water period of 2014, prior to survey efforts in August of that year. Thus, 32 out of the possible 120 ice free days between survey trips to the Toker Point region involved potentially impactful storm events. In order to determine if high water levels accompanied the open-water wind events, tidal gauge records from offshore buoys near Tuktoyaktuk were obtained from the Fisheries and Oceans Canada -Drifting Buoys website (http://isdm-gdsi.gc.ca/isdm-gdsi/drib-bder/index-eng.htm). Much of the buoy data was missing for the 2013 and 2014 period, but high water levels, defined by an increase of 1.25 m above chart datum over a period of no less than one hour (Manson, Solomon 2007), accompanied all of the 2013 wind events from June through October. Several other high water levels from 2013 appear to correlate with wind events that took place as late as mid-November, possibly hinting at a longer open water duration than usual for 2013. Fig. 10: Recent erosion rates at east Toker Point, showing higher rates associated with exposed ground ice.

Discussion
Shoreline erosion has received a great deal of attention from geographers, geologists, engineers and public planners. For example, Shaw et al. (1998) developed a coastal sensitivity index based on observations derived from 1:50K topographic maps. The ranked variables they developed were applied to shoreline regions throughout the circumpolar Arctic in order to determine relative levels of erosion risk. The ability to accurately project such retreat rates into the future can be hindered by a number of processes however, including the fact that such events often take place in a nonlinear fashion. For example, periods of erosion may be followed by periods of stability due to the deposition of protective aprons of eroded material at the base of coastal bluffs (Shaw et al. 1998). While it is possible to predict future trends by accounting for numerous compounding factors, the irregular patterning of natural erosive processes overall can complicate efforts. This is especially true in an Arctic setting, where unique factors such as the activation of backshore thaw slumps can completely invalidate predicted rates over the course of a single warm season.
The Extreme period of erosion at Toker Point between survey visits in 2013 and 2014 has provided a prime example of how increasingly volatile weather patterns can influence the reliability of shoreline erosion models. Understanding the historical trajectory of shoreline change can yield useful insights on where best to direct heritage management efforts. Understanding the weaknesses of historical shoreline change methods is just as critical, however, in achieving reasonable insights which can be reliably applied to said management efforts. Factors such as choice of statistical methods or the manner in which shoreline positional data is processed can have implications for the reliability of results generated. In this instance, checking the accuracy of the model resulted in the recognition of a particularly impactful storm season, and the realization that such influences may be all but invisible without access to annual imagery supplements or concerted ground truthing efforts.
The Toker Point monitoring stakes had been removed between the 2014 and 2015 survey seasons, making it difficult to further assess the quality of shoreline change calculations in the area. Patterned ground in the vicinity of the monitoring stakes allowed for the extrapolation of the 2015 shoreline position from photographs (Fig. 12). The 2014/2015 storm season inflicted 4.5 metres of coastal retreat in front of the monitoring stakes, thereby corroborating the originally calculated erosion rate of -4.27 m/a, and further highlighting the insular nature of the 2013/2014 storm impacts.
From this review of shoreline erosion calculations at Toker Point, it has been demonstrated that using aggregated time periods can yield useful results for generalized, region-wide overviews of historical shoreline change rates. However, the inclusion of further recent imagery and thereby shorter intervals of analysis can result in more nuanced models of shoreline change in the direct vicinity of known archaeological remains. Current forecasts have predicted changes in the seasonality of erosion forcing mechanisms, such that protracted open water periods coincide more often with strong autumn/winter storm events (Manson, Solomon 2007). Thus, the ability to detect sporadic, high-impact erosion events can be of tremendous use to cultural landscape management initiatives both now and in the future.