Haida Gwaii is a remote archipelago located off Canada's mainland Pacific coast ( Figure 1 ) and the recognized Aboriginal title lands of the Haida Nation (Haida Nation and Canada, 2024a, 2024b). The archipelago is characterized by cool, nutrient-rich marine environments which support extensive Nereocystis luetkeana and Macrocystis pyrifera (Lindstrom, 2023 ;syn. Macrocystis tenuifolia) kelp forests (HMTK Participants et al., 2011; Lee et al., 2021). The boundaries of the study region were defined to include the Hlk̲inul ChiiG̲as.sgii SG̲aagiidaay Protection Management Zone designated in the Haida Gwaii Marine Plan (MaPP, 2015) to protect an extensive kelp forest. Our study region of 800 km² on the northeast coast of Moresby Island ( Figure 1 ) spans a sea surface temperature gradient of ~3°C and a strong exposure (fetch, wind and tidal current) gradient, allowing us to represent a wide range of conditions while minimizing the impact of confounding variables when comparing across vast distances. Sea otters were extirpated from Haida Gwaii during the maritime fur trade (1785-1840) and have not yet returned (Lee et al., 2021). This region is characterized by dense, large Macrocystis forests throughout ( Figure 1c ), with smaller dense patches of Nereocystis forests around Cumshewa Head ( Figure 1a ). These kelp forests are detectable from medium- to high-resolution satellites allowing us to monitor kelp forest change back to the 1970s (Gendall et al., 2023). Haida people note the importance of the Cumshewa and Flagstaff areas within the study region as crucial ecological and cultural areas for Macrocystis harvest for the herring-roe-on-kelp fishery (in Haida: k’aaw), however, have observed kelp declines in the early 2000s (HMTK Participants et al., 2011; MaPP, 2021) highlighting the need for increased monitoring and management (MaPP, 2015).

The time series comprises two components: (i) a historical snapshot of century-old kelp presence recorded in British Admiralty charts, and (ii) kelp canopy area derived from satellite imagery collected from 1973 to 2021. The British Admiralty charts data were created between 1867 and 1945 (hereafter historical kelp; Costa et al., 2020). Charts with scales ranging from 1:6,080 to 1:500,000 were scanned, georeferenced, and manually digitized to capture hand-drawn kelp features. In cases where overlapping kelp features were present on multiple charts, in the same region, the highest-resolution charts were prioritized. Each kelp feature was assigned a reliability category based on the depth of the seafloor below the feature, with all features within the study region classified as 'Very High' reliability due to their occurrence in depths shallower than 40 m. These charts, originally designed to map potential navigational risks for mariners, provide only a snapshot of kelp presence and should not be used to infer absences, as not all kelp forests may have been mapped. As such, the historical kelp distribution (Costa et al., 2020) was visually compared with the kelp distribution derived from the satellite imagery time series and was not included in any statistical analyses.

To construct a satellite imagery time series (1973-2021) of floating kelp forest canopy area (hereafter kelp area), we adopted the methodological framework outlined in Gendall et al. (2023). This framework includes: 1) imagery compilation and assessment, 2) imagery preprocessing, 3) classification, and 4) validation (see Supplementary Methods for detailed description of each step). To collect and assess the imagery, we compiled a dataset of archived remote sensing imagery from 1973 to 2021, with resolutions ranging from high (0.5 m) to medium (60 m; Table 1 ). Images were selected based on criteria ( Supplementary Table S1 , see Supplementary Methods ) controlling for cloud cover, tidal conditions ( Supplementary Figure S1 ), haze, waves, glint, and acquisition month (Schroeder et al., 2019; Cavanaugh et al., 2021; Gendall et al., 2023). The final image dataset consisted of 52 suitable images spanning 48 years. Imagery preprocessing included geometric corrections, atmospheric corrections, masking of land, deepwater and soft substrate (Gendall et al., 2023). Steep sloping bathymetry result in a narrow band of depths that are habitable to kelp before becoming too deep, limiting the extent of kelp forests to thin, fringing strips close to shore making detection difficult from medium-resolution satellites (Gendall et al., 2023). To mitigate uncertainties related to mapping small kelp forest canopies at varying resolutions, we excluded these areas of steep sloping bathymetry, exceeding 11.4%, from the analysis (Gendall et al., 2023). However, these steep sloping nearshore areas that were removed only accounted for 0.01% of the overall kelp area in the study region (Gendall et al., 2023).

To improve the detectability of kelp from non-kelp features in the imagery, we generated band indices/ratios specific to kelp ( Table 1 ; Supplementary Table S2 ), which were linearly enhanced and integrated with the original bands for input into the classification process. Satellite imagery was classified utilizing eCognition Developer Software (Trimble Germany, 2011) with enhanced indices/ratios and bands as inputs. The classification results were validated by cross-referencing with available field and historical data (Gendall et al., 2023; Supplementary Table S3 ). When considering all sensors included in the validation, the overall global accuracy ranged between 88% and 94% ( Table 1 ; Supplementary Table S3 ), where most errors were associated with sparse, fringing, or partially submerged kelp canopy (Gendall et al., 2023).

The historical distribution of kelp from approximately a century ago (1867-1945; Figure 2a ) closely mirrored the distribution of the maximum kelp area observed in satellite imagery from 1973 (8.18 km²; Figure 2b ). However, since the historical kelp distribution was represented as artistic, hand-drawn features on maps, direct comparisons of kelp area were not possible.

Our satellite imagery analysis revealed a regional declining trend of 5 ± 2% (± SE) per decade between 1973 and 2021 (p<0.05, R² = 0.11, df = 1,33), with significant variability driven by regional drivers like the ENSO, PDO, SST anomalies and MHW metrics ( Figure 3 ; Table 2 ). Warmer conditions coincided with lower kelp area, with environmental factors explaining up to 29% of the variance in kelp area ( Table 2 ; Supplementary Table S6 ). Our best models identified the ENSO (one-year averages), SST anomalies (one-year to two year averages) and MHW metrics (MHW days and MHW cumulative intensity, one- to two-year sums) as equally good predictors (ΔAICc < 2) of regional kelp area ( Table 2 ; Supplementary Table S6 ).

The time series exhibited periodic losses and recoveries with kelp area initially declining by 75% between 1973 and 1977 ( Figures 2b, c , 3a ) coinciding with the PDO phase shift ( Figure 3b ). This PDO phase shift was marked by an annual average SST anomaly of 0.26°C and a 10 day MHW in may of 1977 ( Figures 3b–d ). Subsequent fluctuations were tied to major El Niño events (e.g., 1983, 1998) and the 2014-2016 ‘Blob’ MHW. In 1998, one of the strongest El Niño events on record occurred, with SST of 1.23°C above average, a total of 193 MHW days and a MHW cumulative intensity of 348.67°C ( Figures 3b–d ). This event aligned with the time series minima, where only 5% of the kelp area (0.43 km²) remained ( Figure 2e ), however recovered to 33% by 1999 ( Figure 2f ).

In the study region, the 2014-2016 ‘Blob’ MHW had minimal impacts on regional SST in 2014, but led to SST anomaly of 0.92°C with a total of 179 MHW days in 2015, and in 2016 a SST anomaly of 0.91°C with a total of 114 MHW days ( Figures 3b–d ). This event corresponded with the second and third lowest kelp areas at 13% and 17% in 2015 and 2016, respectively ( Figures 2g , 3a ). Following the ‘Blob’, SST remained above average for the duration of the study period (2016-2021), however, not as extreme, allowing kelp forests to recover to 48% by 2021 ( Figures 2h , 3a ).

We employed a spatially constrained cluster analysis to group local 1 km segments into subregions based on fetch, local SST, wind, and tidal current. Fetch, wind, and tidal current exhibited high correlations (fetch and wind r = 0.78, fetch and tidal current r = 0.86, wind and tidal current r = 0.80, df = 237, p<0.001). Consequently, only fetch and SST were included in the cluster analysis. The pseudo-F-statistic analysis (Milligan and Cooper, 1985) defined five as the number of optimal subregions (pseudo-f-statistic: 4 groups = 272.92, 5 groups = 303.61, 6 groups = 279.9233; Supplementary Table S7 ). Gray Bay (subregion 1) and Flagstaff (subregion 2) were characterized by warm SSTs and moderate fetch, wind and current speeds ( Figure 4a ; Supplementary Table S7 ). Cumshewa East and West (subregions 3 and 4, respectively) were characterized by cooler SST, however, Cumshewa West had lower fetch, wind and tidal current (i.e. was more sheltered) than Cumshewa East. Mathers Creek (subregion 5) was characterized by warmer SSTs and low fetch, wind and current speeds.

Five subregions exhibited varying responses to climatic events, largely driven by differences in SST. Warmer more exposed subregions like Gray Bay and Flagstaff experienced drastic declines, with complete kelp loss in Gray Bay after 1977 ( Figures 4a, b ). Flagstaff showed a long-term decline of 9 ± 2% per decade (p<0.001, R² = 0.34, df = 1,33), with minimal recovery after the loss in the late 1970s ( Figures 2c, d , 4c ). Three-year SST averages best explained the variation (42 %) in the Flagstaff subregion (Table 2; Supplementary Table S8 ). By the early 21st century, approximately 90% of the original kelp area in the Flagstaff subregion had disappeared ( Figure 4c ).

In contrast, the cooler Cumshewa East and Cumshewa West subregions displayed resilience, with no significant long-term declining trends ( Figures 4a, d, e ). Kelp area in both subregions still fluctuated with environmental conditions, with ENSO (one-year and two-year averages) and MHW metrics (MHW days and MHW cumulative intensity, one- to three-year sums) explaining up to 25% (Cumshewa East) and 34% (Cumshewa West) of variation in kelp area (Table 2; Supplementary Tables S9, S10 ). However, unlike the northern subregions, these cooler subregions demonstrated a greater capacity for recovery following events, such as the 1977 PDO shift, the 1983 El Niño, and the 1998 El Niño ( Figures 4d, e ). Kelp area declined drastically during the 2014-2016 MHW in both subregions, however, Cumshewa West rebounded to approximately 77% (2017-2021), while Cumshewa East only reached 28% ( Figures 4d, e ).

The Mathers Creek subregion, despite being geographically warmer like Flagstaff and Gray Bay subregions, exhibited resilience similar to Cumshewa East and Cumshewa West ( Figures 4a, f ). In the Mathers Creek subregion, the PDO (one-year averages) and ENSO (one- and two-year averages) explained up to 24% of the variation in kelp area (Table 2; Supplementary Table S11 ). This subregion supported the lowest proportion of kelp area (maximum of 0.21 km² in 1990; Figure 4f ), which was largely located around the mouth of Mathers Creek in the cooler portion of the subregion.

At the local scale, variability of kelp persistence was significantly explained by local conditions of SST, wind, tidal current and depth where less wind, shallower depths and lower tidal currents coincided with higher persistence ( Supplementary Table S12 ). Specifically, our best model included SST, wind and depth, which explained 41% of variation in kelp persistence ( Supplementary Table S12 ). In Gray Bay, shallow water segments (2-5 m depth) experienced the warmest SST and exhibited no kelp persistence (0%; Figure 5a ). Segments further offshore in Gray Bay showed low persistence (3-5%) and coincided with areas where kelp forests had entirely disappeared during the late 1970s PDO shift. In the Flagstaff subregion, kelp canopy persistence closely mirrored the SST gradient, where the warmer northern segments displayed low persistence (3-31%) with no recovery, and cooler southern segments demonstrated high persistence (67-95%), regularly recovering after climate events when cross referenced with the time series. In the more resilient subregions, the more sheltered nearshore segments in Cumshewa East (3), Cumshewa West (4), and Mathers Creek (5) maintained persistent kelp forests. Further offshore within the shallower mid-channel areas of Cumshewa East (3) and Cumshewa West (4), a core kelp area known as Fairbairn Shoals remained persistent throughout the time series.

Where kelp was lost, field surveys suggested the habitat transitioned to turf-dominated reefs rather than urchin barrens. In particular, we observed only five individual urchins within the study region (n = 48), and in subregions Gray Bay and Flagstaff, where kelp was lost, the understory community was predominantly composed of short turf seaweeds and branched brown algae such as Desmarestia spp. on subtrate primarily composed of a mix of granule and cobble ( Figures 5b, c ). In the Flagstaff subregion, only one site where kelp was lost supported an understory kelp community; this site was closest to the remaining kelp forests. In contrast, where kelp canopy persisted, an abundant understory of kelp species were found on a variety of substrate types from pebble to boulders and bedrock ( Figures 5b, c ).

The inclusion of historical data and medium resolution satellites revealed substantial kelp forest losses that would have gone undetected using only high-resolution satellite imagery. Specifically, we discovered a major episodic decline in the 1970s of kelp forests that had likely persisted for over a century – a finding that emerged only through our analysis of pre-1984 data sources, including nautical charts (1867-1945) and early Landsat imagery. The creation of long-term time series are crucial, as shorter time series (< 20 years) can mask important patterns in kelp forest resilience due to high interannual variability (Bell et al., 2020; Wernberg et al., 2019). While most remote sensing of kelp canopy on the Pacific Coast of North America focuses on 30 m Landsat imagery from 1984 onwards (e.g., Cavanaugh et al., 2011; Bell et al., 2020; Hamilton et al., 2020), researchers have successfully integrated diverse historical records to establish long-term baselines. Examples include the use of early 20th century kelp census data in Washington (Pfister et al., 2018), historical surveys from unpublished theses and aerial photography dating to the 1930s in South Australia (Carnell and Keough, 2019), Darwin's observations on the Voyage of the Beagle in Southern Chile (Mora-Soto et al., 2021), and nineteenth-century nautical charts in the Salish Sea (Mora-Soto et al., 2024b) – the latter revealing similar patterns of kelp loss in warmer regions. This underscores how incorporating historical data sources, despite their varying quality, resolution and methodologies, is essential for establishing accurate pre-industrial ecological baselines and understanding the true magnitude of changes that may otherwise be obscured by recent data alone.

Our findings emphasize that regional drivers, such as the PDO, ENSO, SST anomalies, and MHW metrics, are critical in explaining kelp forest dynamics, with impacts extending beyond immediate SST stress, with multi-year metrics explaining more variation in kelp area than single-year averages. Many studies corroborate that these regional drivers at multiple scales, from MHW metrics to large scale climate cycles like the PDO and ENSO, impact kelp forest conditions, where warm conditions are detrimental to kelp forests (e.g., Cavanaugh et al., 2011; Bell et al., 2020; Arafeh-Dalmau et al., 2019; Wernberg et al., 2019; McPherson et al., 2021; Wernberg, 2021). The extended multi-year impacts likely result from complete mortality of whole kelp plants rather than surface die back, as mortality reduces the reproductive capacity of the remaining kelp forest even after favorable conditions return (Cavanaugh et al., 2011; Pfister et al., 2018; Schroeder et al., 2019; Starko et al., 2022).

Regional drivers not only influence SST but also interact with other oceanic conditions exacerbating not only changes in kelp forests but the ecosystem as a whole. For instance, the significant kelp declines in the 1970s coincided with the 1977 PDO shift, which brought warmer SSTs, a deepening of the thermocline, low nutrient availability and increased ocean stratification across the North Pacific (McGowan et al., 2003; Parnell et al., 2010). This shift triggered a cascade of declines in other coastal systems, particularly in California such as reduced plankton biomass, larval fish abundance, fisheries landings, and similarly, Macrocystis biomass (McGowan et al., 2003; Parnell et al., 2010). Notably, the documented loss of kelp forests in this study predates the widespread coral reef declines attributed to the 1983 El Niño (Oliver et al., 2009) challenging the notion that kelp forests are not sensitive indicators of broader climate change impacts (Reed et al., 2016).

Local heterogeneity in factors like SST can either ameliorate or exacerbate the responses of foundation species like kelp forests to regional drivers, leading to complex patterns of decline and resilience (Russell and Connell, 2012; Starko et al., 2022). In our study region, cooler areas, like Cumshewa East and West, contribute to this resilience, as kelp forests there exhibited a stronger ability to endure and recover after adverse regional conditions. This spatial pattern of temperature-dependant responses has been documented elsewhere on the Pacific Coast, with losses concentrated in warmer areas during the 2014-2016 ‘Blob’ MHW in Barkley Sound (Starko et al., 2022) and the southern Salish Sea (Mora-Soto et al., 2024a). The complex coastline of British Columbia and Alaska creates high environmental heterogeneity (Starko et al., 2019; Cavanaugh et al., 2021), making it crucial to consider local factors such as temperature, wave exposure, and trophic dynamics (e.g., Watson and Estes, 2011; Rogers-Bennett and Catton, 2019; McPherson et al., 2021; Starko et al., 2022, 2024) to understand kelp forest variability and accurately forecast species responses to climate change moving forward (Russell and Connell, 2012).

In addition to temperature, local patterns of exposure, wind and current can either ameliorate or exacerbate the response of kelp forest to regional drivers. In more exposed locations, storms can tear up large swaths of kelp forests like the intense storms during the 1983 El Niño in Point Loma, California, that led to a loss of 560 ha of kelp forests (Dayton and Tegner, 1984). Recent work shows that storm tracks have intensified on the coast of BC from the 1960s onwards (Abeysirigunawardena, 2010) and with the significant negative relationship between kelp persistence and wind in the study region, storms could have compounded or caused some of the loss or variation not explained by regional or local drivers in this study.

Multiple lines of evidence, the lack of urchins, and presence of turf reefs, in our 2021 field survey suggest that temperature, rather than herbivory, likely drove kelp forest declines in the study region. Despite the presence of urchin barrens (HMTK Participants et al., 2011; Lee et al., 2021), and the absence of sea otters since the late 1800s (Lee et al., 2021) in Haida Gwaii, our field survey found no urchin barrens and very few individual urchins in areas of kelp loss and persistence. The absence of urchin barrens is particularly notable as they are known to persist for years to decades due to urchins’ ability to slow metabolic activity when food resources are scarce (Spindel et al., 2021). Additionally, the patchy mix of unconsolidated substrate – consisting of sand, pebbles, cobbles and boulders, is known to dissuade urchins from foraging (Laur et al., 1986) and may limit the distribution of urchins across the study region. Notably, in Gray Bay, where kelp was never present, both turf and eelgrass dominated shallow areas on sand and pebbles. This distribution pattern aligns with known habitat requirements, as canopy forming kelp require larger hard substratum such as cobbles, boulders or bedrock to provide secure attachment points (Druehl, 1978; Gregr et al., 2019; Starko et al., 2022), suggesting that substrate type likely plays a key role in determining kelp forest distribution across the region. Additionally, in areas where kelp disappeared, we found predominantly turf and branched algae (Desmarestia spp.) – a transition pattern commonly observed from temperature driven kelp losses in other regions of the globe, such as Western Australia (Wernberg, 2021), Atlantic Canada, France and Denmark (Filbee-Dexter and Wernberg, 2018). While our conclusion about herbivore impacts in this system is supported by two lines of evidence, our single time-point highlights the need for further research to fully understand trophic dynamics in this ecosystem.

The regional and local drivers included in this study did not account for all variation observed in kelp area and persistence. Variability in other potential drivers, such as nutrients (e.g., Zimmerman and Kremer, 1986), light availability (e.g., Deysher and Dean, 1986), salinity (e.g., Druehl, 1978) and local stressors from human activities like pollution or overharvesting, are known to impact kelp forests (e.g., Foster and Schiel, 2010; Pfister et al., 2018), but data were unavailable for the region. Warmer SST often coincides with lower nutrients and salinity (Druehl, 1978), which generally slows kelp growth, weakens tissue and diminishes reproduction (Zimmerman and Kremer, 1986). As such, nutrients likely play an important role alongside SST in regulating kelp dynamics but were unavailable in this region. Although light data were unavailable, the significant relationship we found between depth and kelp persistence likely reflects the influence of light availability, as light attenuation increases with water column depth (Timmer et al., 2022).

While coastal development, point-source pollution and overharvesting can drive kelp declines (Foster and Schiel, 2010; Pfister et al., 2018), anthropogenic stressors are likely minimal in Haida Gwaii given its small population (< 5000 people) and minimal development (Statistics Canada, 2017) when compared to densely populated coastal areas with observed kelp declines like in Washington or South Australia (Krumhansl et al., 2016; Pfister et al., 2018). Although Macrocystis kelp was historically harvested for the herring-roe-on-kelp commercial fisheries in the Flagstaff subregion, BC’s strict harvest limits leave the majority of the kelp plants intact minimizing harvest impacts. While turbidity related to logging runoff on Haida Gwaii (Klein et al., 2012) could impact kelp, the persistence of healthy kelp forests at the mouth of Mathers Creek—the region's only river—suggests minimal impacts from both turbidity and reduced salinity. Instead, the river's outflow may actually buffer SST impacts as evidenced by the retention of kelp forests in the Mathers Creek subregion despite losses in the similarly warm Flagstaff and Gray Bay subregions.

Previous work suggests that interannual changes in kelp area up to 7% could potentially arise from errors associated with measuring kelp area at different spatial resolutions (Gendall et al., 2023). While measurement uncertainty is an inherent challenge in remote sensing of kelp forests, the magnitude of kelp decline we observed far exceeds the potential error threshold identified by Gendall et al. (2023). Additionally, our focus on long-term trends and persistence over multiple decades helps mitigate the influence of such errors on our overall findings (Magurran et al., 2010). This approach to analyzing temporal patterns is sound as it emphasizes directional changes over long time scales rather than year-to-year fluctuations that might be influenced by artifacts or uncertainties in remote sensing measurements. Additionally, our use of multiple temporal metrics—including both area changes at the regional and subregional scale and persistence patterns at the local scale—provide complementary lines of evidence that strengthen our confidence in the patterns of kelp forest change herein.