More Than Marine Heatwaves: A New Regime of Heat, Acidity, and Low Oxygen Compound Extreme Events in the Gulf of Alaska

Recent marine heatwaves in the Gulf of Alaska have had devastating impacts on species from various trophic levels. Due to climate change, total heat exposure in the upper ocean has become longer, more intense, more frequent, and more likely to happen at the same time as other environmental extremes. The combination of multiple environmental extremes can exacerbate the response of sensitive marine organisms. Our hindcast simulation provides the first indication that more than 20% of the bottom water of the Gulf of Alaska continental shelf was exposed to quadruple heat, positive hydrogen ion concentration [H+], negative aragonite saturation state (Ωarag), and negative oxygen concentration [O2] compound extreme events during the 2018–2020 marine heat wave. Natural intrusion of deep and acidified water combined with the marine heat wave triggered the first occurrence of these events in 2019. During the 2013–2016 marine heat wave, surface waters were already exposed to widespread marine heat and positive [H+] compound extreme events due to the temperature effect on the [H+]. We introduce a new Gulf of Alaska Downwelling Index (GOADI) with short‐term predictive skill, which can serve as indicator of past and near‐future positive [H+], negative Ωarag, and negative [O2] compound extreme events near the shelf seafloor. Our results suggest that the marine heat waves may have not been the sole environmental stressor that led to the observed ecosystem impacts and warrant a closer look at existing in situ inorganic carbon and other environmental data in combination with biological observations and model output.

exposed to just one single stressor (Breitberg et al., 2015;Kroeker et al., 2021;Thomsen et al., 2013).At the same time, abrupt extremes will have different consequences for the adaptation potential of organisms to global changes than evolutionary response to the gradual long-term changes, such as warming, ocean acidification, and deoxygenation (Bell et al., 2021).This is even more complex for subsequent or compound extreme events that can lead to conflicting selection pressures and decrease the genetic diversity and adaptation potential (Gaitán-Espitia et al., 2017).
The Gulf of Alaska marine ecosystem has been exposed to both marine heat waves and ocean acidification extreme events.However, compound extreme events with multiple environmental conditions outside of their historic variability envelope have not been documented yet.Living marine resources from the Gulf of Alaska not only sustain economically important seafood and tourism industries, but they also play a crucial role in the way of life of Indigenous communities.Recent marine heatwaves have triggered devastating and lasting responses at various trophic levels, from primary producers to commercially and subsistence caught fish species in different regions across the Gulf of Alaska marine ecosystem (Barbeaux et al., 2020;Bellquist et al., 2021;Ferriss & Zador, 2022;Piatt et al., 2020;Rogers et al., 2021;Suryan et al., 2021;Von Biela et al., 2019;Weitzman et al., 2021).The longest marine heat wave to date, called "the Blob," occurred between 2013 and 2016 (Bond et al., 2015;Di Lorenzo & Mantua, 2016;Hobday et al., 2018).This 2013-2016 marine heat wave was followed by another marine heatwave in 2018-2020 (Amaya et al., 2020;Barkhordarian et al., 2022).The 2013-2016 marine heat wave was initiated by a strong atmospheric ridge over the northeast Pacific in the winter of 2013/2014 that weakened the Aleutian Low and surface winds and caused sea surface temperature changes that project on the Pacific Decadal Oscillation (PDO) pattern (Bond et al., 2015).A combination of reduced Ekman transport of cold water from the north, decreased upper ocean mixing, and weakened surface heat loss led to a warmer surface ocean.The 2018-2020 marine heat wave (Amaya et al., 2020;Barkhordarian et al., 2022) was primarily driven by a weak North Pacific High, which induced weak surface winds and less evaporative cooling.It started in the summer during a highly stratified period, inducing strong temperature anomalies at the surface (Amaya et al., 2020).
Observations and hindcast model output suggest that the Gulf of Alaska ecosystem has also been exposed to ocean acidity extreme events (Bednaršek et al., 2021;Hauri et al., 2021).Evidence of severe dissolution in pteropod shells, which are an indicator species for ocean acidification, were found in the Gulf of Alaska with concomitant unusually acidified conditions (Bednaršek et al., 2021).These ocean acidity extreme events were likely triggered by a combination of increased upwelling of carbon dioxide (CO 2 )-rich water in the Alaskan gyre and ocean acidification from rising atmospheric CO 2 (Hauri et al., 2021).The upwelling strength of the gyre is driven by decadal variability of the local wind stress curl that depresses sea surface height and has been described as the Northern Gulf of Alaska Oscillation (NGAO, Hauri et al., 2021; Figure 1).
There is overwhelming evidence that marine ecosystems and their associated services are under threat because of local (e.g., over-fishing) and global (e.g., heat waves) changes (Cooley et al., 2022).These pressures combine in a unique way at each location and identifying the major ocean stressor or combination of multiple stressors driving the biological response is critical for the implementation of solutions.For example, addressing local effects of ocean acidification would require a combination of global CO 2 mitigation and local adaptation solutions (e.g., management to increase ecosystem resilience; IOC-UNESCO, 2022).However, attributing observed biological changes over time in an ecosystem to one or multiple stressors is not an easy task (Widdicombe et al., 2023).
Here, we used output from a regional ocean biogeochemical model that simulated the environmental conditions in the Gulf of Alaska from 1993 through 2021 to study the occurrence and drivers of extreme and compound extreme events at the surface and near the shelf seafloor.

Model Set-Up
This study was based on an ocean hindcast simulation  with the Gulf of Alaska configuration (Figure 2; Hauri et al., 2020Hauri et al., , 2021) ) of the three-dimensional physical Regional Oceanic Modeling System (ROMS, Shchepetkin & McWilliams, 2005) coupled to the 3PS Carbon, Ocean Biogeochemistry, and Lower Trophic marine ecosystem model (3PS COBALT; Stock et al., 2014;Van Oostende et al., 2018).The main characteristics of this configuration were the 50 terrain-following depth levels and an eddy-resolving 10.1029/2023AV001039 3 of 18 horizontal resolution (4.5 km) that resolves the regional coastal upwelling and downwelling.3PS-COBALT simulates the cycles of nitrogen, carbon, phosphate, silicate, iron, calcium carbonate, oxygen, and lithogenic material with 36 state variables.More details on the model setup can be found in Hauri et al. (2020) and references therein.(Hauri et al., 2021).The amount of variance associated with the first EOF of the model is ∼24%.(Sculley, 2010) are illustrated by lines.The K-MEANS algorithm was carried out with the scikit-learn Python package (Pedregosa et al., 2011).Our current analysis focuses on the gyre (red line) and shelf (magenta line) regions.

Initial and Boundary Conditions
Physical ocean initial and boundary (daily) conditions (temperature, salinity, and sea surface height) were taken from the Global Ocean Physics Reanalysis (GLORYS; Hamon et al., 2016) until the end of 2020 and came from Global Ocean Physics Analysis and Forecast (2023) (accessed in January 2023, https://doi.org/10.48670/moi-00016)after that.Output from the Japanese 55-year re-analysis (JRA55-do) version 1.5 (Tsujino et al., 2018) was used as forcing for wind, air-surface temperature, pressure, humidity, precipitation, and radiation at a three-hourly resolution.Data from the NOAA Greenhouse Gas Marine Boundary Layer (Lan et al., 2023) was used as forcing for atmospheric partial pressure of CO 2 (pCO 2 ).Initial and boundary conditions for nitrate, phosphate, oxygen, and silicate were provided by the World Ocean Atlas 2018 (Boyer et al., 2018).Initial and boundary conditions for iron were provided by a yearly climatology of the GLORYS biogeochemistry reanalysis product (Perruche, 2019).The iron atmospheric soluble deposition was derived from Geophysical Fluid Dynamics Laboratory global products and reduced by 50% to match the order of magnitude provided by Crusius et al. (2017).Riverine freshwater discharges were the same as in Hauri et al. (2020) and derived from Beamer et al. (2016) and Hill et al. (2015).Riverine iron input was set to 100 nM to be more consistent with the value available in Lippiatt et al. (2010).All other river biogeochemical variables are described in Hauri et al. (2020).The chlorophyll to carbon ratio of the different phytoplankton classes has been updated to better simulate the ecosystem of the Gulf of Alaska (Table 1).A 10-year spin-up was run based on a 4-year (1993)(1994)(1995)(1996) simulation repeatedly run in a loop.After 10 years the biogeochemical variables reached an approximate seasonally varying steady-state.

Model Evaluation
Modeled surface and subsurface fields have been evaluated in previous studies, suggesting that the model simulates the general physical and chemical oceanographic patterns (Dobbins et al., 2009;Hauri et al., 2020Hauri et al., , 2021)).Our ROMS simulation captured well the large-scale observed seasonal and interannual SST variability patterns (Figure 3) as is typically expected for hindcast simulations forced by atmospheric reanalysis and bulk heat and moisture flux parameterizations (Doney et al., 2007).Point by point comparison of model output with in situ dissolved inorganic carbon (DIC), total alkalinity (TA), and oxygen data from various cruises (Figure 4) shows a relatively high Pearson correlation (>0.8), with statistically significant p-values (<0.01).Most normalized standard deviations fall close to 1, with low normalized root mean square error, suggesting that the model represents the observations reasonably well.

Definition of Extreme, Compound Extreme, and Triple or Quadruple Compound Extreme Events
Due to the strong natural variability of the Gulf of Alaska, we chose to use relative thresholds to define the extreme events, assuming that organisms might be adapted to the wide range of variability in their environment.This analysis was based on the daily output of the 29-year-long simulation.First, a 14-day running mean was applied to the daily model output (Figures 5a and 5b).Smoothing over 7 days or 30 days did not change the results, which was also found in Le Grix et al. (2021).As a next step, the seasonal anomalies were computed for each grid cell, where the smoothed data was deseasonalized by removing a daily resolution climatological cycle constructed from averaging the available 29 years of model output (Figure 5c).The daily relative thresholds were defined as the 95th percentile (Figure 5d) for temperature and [H + ], and the 5th percentile for aragonite saturation state (Ω arag ) and [O 2 ] of the distribution of daily anomalies for the full multi-decadal time period (Burger et al., 2022).An event was considered as extreme if the anomalies cross the threshold (Figure 5e).This method implies that per a given grid cell the frequency of extreme events is 5% across the timeseries.A compound extreme event was defined as the co-occurrence of two (or more) extreme events in time and space (i.e., the same grid cell).Over 10,000 days (and approximately 1,000 independent time points with 14-day running mean) per grid cell were used to compute percentiles between 1993 and 2021, making the analysis statistically robust.A fixed temporal baseline was used because the strong natural climate variability of the area can mitigate or accelerate Note.The Chl/C for the large phytoplankton was based on Gulf of Alaska in situ values (Coyle et al., 2019;Strom et al., 2010).Chl/C for medium was reduced by the same amount (28%) as large phytoplankton to remain smaller and is now at the higher range of the value given by Sathyendranath et al. (2009).Chl/C for small phytoplankton and diazotrophs was reduced to what was found by Sathyendranath et al. (2009).

Table 1
Table Showing the Chlorophyll to Carbon Ratio (Chl/C, in g/g) for Each COBALT Phytoplankton Group as Used in Hauri et al. (2020) and in the Current Model Version the apparent rate of ocean acidification (Hauri et al., 2021), making a moving baseline hard to apply over these relatively short periods of time (less than 30 years).Fixed baselines have been used in several recent studies of extreme events (Barkhordarian et al., 2022;Burger et al., 2022;Desmet et al., 2023;Laufkötter et al., 2020;Le Grix et al., 2021) and are generally well suited to investigate the effect of extreme events on organisms (Oliver et al., 2019).No minimum duration threshold for extreme events was set as the minimum duration is unknown (Collins et al., 2019;Le Grix et al., 2021).

[H + ] at Fixed Temperature
[H + ] was computed at a fixed seasonal temperature to highlight the importance of temperature in triggering [H + ] extreme events.At each grid cell, [H + ] at a fixed temperature (H + Tfix ) was computed based on the model annual-mean temperature climatology calculated at each grid cell and modeled DIC, TA, and salinity (PyCO-2SYS python package; Humphreys et al., 2022).PyCO2SYS uses the same equations as COBALT to resolve the carbonate cycle.

Extreme Event Association With Various Climate Indices
To investigate potential relationships of extreme and compound extreme events with different empirical modes of variability (short: "climate modes") we followed the approach of Holbrook et al. (2019) and Le Grix et al. (2021).We computed the frequency of marine heat, positive [H + ], and negative Ω arag and [O 2 ] extreme and compound extreme events during the positive, negative, and neutral phases of the NGAO and a newly defined Gulf of Alaska downwelling index.The positive phases of each index were associated with days when the indices were above 50% of their maximum, whereas the negative phases of each index were associated with days the indices were below 50% of their minimum, and the neutral phases correspond to the days when the indices were between 50% 10.1029/2023AV001039 6 of 18 of their maximum and 50% of their minimum.To determine if a given climate mode influenced extreme event frequencies, we compared the frequency during the positive and negative phases to the frequency during the neutral phase at each grid cell.Finally, to verify the statistical significance we randomly shuffle the temporal order of each index and recomputed the frequency 1,000 times.The result was considered statistically significant if the observed frequency was higher than 95% or lower than 5% of the shuffled cases.

Definition of the Areas Used for the Extreme and Compound Extreme Events
The areas used in this study were defined based on the spatial pattern of the modeled de-seasonalized sea surface height in the Gulf of Alaska.The seasonal signal was removed from monthly model output and a classical K-MEANS algorithm was applied (Sculley, 2010).
The K-MEANS algorithm was configured with four clusters, including the land mask with an assigned dummy value.A silhouette analysis was performed to choose this number of clusters.Prince William Sound, Cook Inlet, and the Juneau areas were discarded since the model has not been specifically designed and evaluated for these near shore areas.The area close to the model lateral boundary (within 50 km) was also discarded to avoid the direct effect of the boundary conditions on the results.The K-MEANS algorithm defined the following four main areas (Figure 2): (a) the center of the gyre that corresponds to the lowest mean SSH and therefore the maximum upwelling intensity, (b) the buffer zone between the center of the gyre and the shelf, (c) the continental shelf area that corresponds to the area with generally positive SSH (coastal downwelling), and (d) the land zone.
The surface area affected by an event (in %) was computed every day for the shelf and gyre regions.The intensity corresponds to the spatial average (over all the grid cells affected by an event) of variable anomalies during an event, every day.The bottom waters were defined as the bottom area on the shelf with depths >50 and <250 m.Finally, the number of days affected by an event during the 2013-2016 and 2018-2020 periods were calculated based on the spatial and temporal average over the gyre and shelf areas.The average takes only grid cells into account that have been affected by at least one day of an event.

Empirical Orthogonal Function
Empirical orthogonal functions (EOF) were used to determine the modes of variability of the Gulf of Alaska by decomposing the oceanic field into a set of uncorrelated spatial modes and their corresponding temporal variations or principal components (PCs).The EOF was performed across the model domain, which covers the Gulf of Alaska with its southernmost point located offshore of Vancouver Island (48.2°N, 148°W) and reaches from the southern tip of Prince of Wales Island in the southeast (55.1°N, 133°W) to the west of Sand Point in the middle of the Aleutian Islands (55.2°N, 163°W).Following Hauri et al. (2021), the EOF was applied to the sea surface height daily model output after removing the long-term temporal trend and deseasonalizing the data.The results showed that the leading EOFs of SSH were well separated based on their eigenvalues with PC1 = 49%, PC2 = 21%, and PC3 = 12%.

Surface Heat and Acidity Compound Extreme Events in the Gulf of Alaska
The Gulf of Alaska ecosystem experienced the first heat and positive [H + ] compound extreme events during the 2013-2016 marine heat wave (Figures 6a and 7a).While several marine heat waves occurred throughout the simulation, the 2013-2016 and 2018-2020 marine heat waves lasted longer and were more extreme than the previous ones, which aligns with previous observational and modeling studies (Amaya et al., 2020;Di Lorenzo & Mantua, 2016).The 2013-2016 marine heat wave was preceded by a 7-year-long cold phase (Danielson Figure 5. Definition of extreme events.Raw daily sea surface temperature (SST) model output from the gyre area (a) was smoothed with a 14-day running mean (b).Anomalies (c) per grid cell were then used to compute the relative thresholds based on the 95th percentile for temperature (d).If an event crosses the relative threshold, it is considered as extreme (e).et al., 2022), with frequent low Ω arag (Figures 6c and 7c) and high [H + ] (Figures 6b and 7b) extreme events that affected a large proportion of the gyre area (Figures 7a-7c; Hauri et al., 2021).Negative Ω arag extreme events also affected up to 95% of the surface water on the shelf at times (Figure 6c), whereas positive [H + ] extreme events were not as pronounced there (Figure 6b).The rapid succession of positive [H + ], negative Ω arag , and heat extreme events that started in the early 2000s was followed by several marine heat and positive [H + ] compound extreme events, frequently exposing the ecosystem to environmental conditions outside their norm.During the 2013-2016 marine heat wave, the gyre area experienced on average 19 (±15, STD) days of these compound extreme events and up to 88 days for some grid cells.The marine heat and positive [H + ] compound extreme events were more prevalent on the shelf, with a spatial mean of 55 (±32) days and a maximum of 146 days.During the peak of the 2013-2016 marine heat wave, when the intensity was highest and the entire shelf area was affected by extreme heat, 60% of the shelf area experienced a marine heat and positive [H + ] compound extreme event.The 2018-2020 marine heatwave already started in late 2018 (Figures 6a and 7a) and therefore earlier in the Gulf of Alaska than in other parts of the North Pacific, which corresponds with observations by Amaya et al. (2020).The number of days affected by marine heat and positive [H + ] compound extreme events over the gyre area increased twofold during the 2018-2020 marine heatwave compared to the 2013-2016 marine heat wave.During the 2018-2020 marine heatwave, up to 50% of the gyre region was affected by marine heat and positive [H + ] compound extreme events, which lasted on average 41 (±32) days, and up to 4% of the shelf ecosystem was exposed over an average of 59 (±39) days.In contrast, short extreme heat and negative Ω arag compound extreme events only occurred during the 2018-2020 marine heat wave and affected small portions of the shelf (<10%, Figure 6e).

Temperature Induced Early Onset of Heat and Positive [H + ] Compound Extreme Events
While [H + ] and Ω arag both follow a secular trend as a result of ocean acidification (Hauri et al., 2021), the intensity and extent of positive [H + ] and negative Ω arag extreme events were spatially and temporally decou- pled as a result of their sensitivity to different climate change and natural variability-affected physical drivers (Figures 6-8).The natural climate variability of the offshore upwelling, described by the local climate index NGAO, directly affects extreme and compound extreme events at the surface (Figure 8).For example, during negative NGAO events (2006-2012, Figure 1) strong offshore upwelling of cold, CO 2 -rich subsurface water in the Alaska gyre leads to both positive [H + ] and negative Ω arag extreme events, although negative Ω arag extreme events were more intense and covered a wider area in the gyre (Figures 7b and 7c; Hauri et al., 2021).Negative Ω arag extreme events also occurred on the shelf (Figure 6c), whereas positive [H + ] extreme events were less common there (Figure 6b).It is important to understand the underlying cause of positive [H + ] extreme events.Because of the complexity of seawater acid-base chemistry, elevated acidity can occur either because of more CO 2 -rich water or because of warmer temperature with fixed inorganic carbon and alkalinity.At fixed seasonal temperature, positive [H + Tfix ] extreme events increased in the gyre (Figure 7d) and expanded onto the shelf (Figure 6d), following a similar pattern as the negative Ω arag extreme events.This suggests that during a negative NGAO phase, offshore upwelling of cold subsurface water (Figure 9) dampened the signal of [H + ]-rich upwelled water and thus the strength and extent of positive [H + ] extreme events in the gyre and on the shelf.
Elevated temperature during the positive NGAO phase amplified elevated [H + ] concentrations to trigger the first surface marine heat and positive [H + ] compound extreme events that were particularly strong in the shelf region.During a positive NGAO phase (e.g., marine heat waves in 2013-2016 and 2018-2020, Figure 1) upwelling was less intense, which led to warmer (Figure 9) and less acidic conditions in the gyre.This natural mechanism preconditioned the water and caused the early onset of the 2018-2020 marine heat wave in the Gulf of Alaska.Simultaneously, it weakened positive [H + ] and negative Ω arag extreme events in the gyre region during the positive NGAO phase (Figures 8b and 8c).However, elevated temperature in the positive NGAO phases (during the 2013-2016 and 2018-2020 marine heat waves) increased [H + ] concentrations enough to cause [H + ] extreme events, triggering the first marine heat and positive [H + ] compound extreme events.Temperature was warmer on the shelf than in the gyre region (Figure 3a), leading to a stronger compound event there (Figures 6e and 7e).At fixed temperature, positive [H + Tfix ] extreme events did not occur during positive NGAO phases (Figures 6e, 7e, and 8d), highlighting the major role of the temperature effect.The generally colder conditions in the gyre due to persistent upwelling activity (Figures 3a and 3b) tampered [H + ] events, even during NGAO positive phases, and resulted in generally weaker compound extreme events.The NGAO positive phases favor the occurrence of compound extreme events (Figure 8e) everywhere in the Gulf of Alaska.However, the strong percentage changes in the gyre region (Figure 8e) are due to the low frequency of compound events during neutral NGAO phases (most likely due to the low number of compound events) and should be seen as an indicator rather than an absolute value.While seasonal surface Ω arag on the shelf can be affected by the rather acidic riverine freshwater input (Hauri et al., 2020), freshwater runoff exhibits only small interannual variability and the correlation with negative Ω arag extreme events was weak (−0.05, p-value < 0.01).

Quadruple Compound Extreme Events on Shelf Seafloor Influenced by Anthropogenically Amplified CO 2 Concentrations in Deep Water Intrusion
The environment in the bottom waters of the shelf area also experienced extremely high [H + ], and extremely low Ω arag and [O 2 ] conditions during the 2013-2016, and more intensely during the 2018-2020 marine heat waves, causing quadruple compound extreme events.Sub-surface marine heat extreme events occurred throughout the timeseries but became more prevalent during the 2013-2016 marine heat wave and affected up to 80% of the shelf seafloor (Figure 10a).This subsurface marine heat wave started in 2015 and therefore later than at the surface, which is consistent with in situ observations (Danielson et al., 2022).The timing of our modeled heat waves on the shelf seafloor is similar to the timing of events found in an ocean reanalysis study (Amaya et al., 2023), although they used the 90th percentile and considered a larger shelf area (bottom grid cells shallower than 400 m), which could have led to the differences in the spatial extent and duration of the heat waves.
Triple positive [H + ], negative Ω arag , and negative [O 2 ] compound extreme events covered >20% of the shelf seafloor area in 2005 and increased in frequency, intensity, and spatial extent over time (Figure 10f), and especially during the 2013-2016 and 2018-2020 marine heat waves.These triple compound extreme events affected up to 38% and 47% of the shelf seafloor area over an average of 108 (±65) and 135 (±68) days, respectively, between 2013-2016 and 2018-2020.These triple compound extreme events were triggered by intrusion of deep, salty, CO 2 -rich, and O 2 -poor water, which can be approximated through the intensity of coastal downwelling and therefore is indicated by sea surface height variations.The second mode of variability of sea surface height in our model domain, here defined as the Gulf of Alaska downwelling index (GOADI), can be used as a proxy for the coastal downwelling strength (Figure 11).A positive phase indicates increased sea surface height and therefore strong downwelling and weak or no intrusion of deep and colder water onto the shelf.A negative phase shows relaxation of downwelling and intrusion of deeper and colder water onto the shelf.Higher salinities on the shelf have been linked to relaxation of downwelling and deep-water intrusion onto the shelf (Danielson et al., 2022;Ladd et al., 2005;Stabeno et al., 2016).Modeled positive salinity events are highly correlated with positive [H + ] with a Pearson correlation coefficient (Pcc) of 0.65 (p-value < 0.01), and with negative Ω arag (Pcc = 0.63, p-value < 0.01) and [O 2 ] (Pcc = 0.8 p-value < 0.01) and exhibit the same temporal variability even though their intensity and affected area can slightly differ.During a negative phase (weak downwelling activity), the frequency of positive [H + ] and negative Ω arag /[O 2 ] triple compound extreme events strongly increased (Figures 12b-12f).In contrast, the heat extreme events were dampened during this phase (Figure 12a).For example, the intense 2016 marine heatwave along the shelf seafloor was weakened as the downwelling relaxed toward the end of the year (downwelling index flips from strongly positive to negative), which suggests that intrusion of cold water alleviated this sub-surface marine heat wave.During a positive phase (strong downwelling activity), positive [H + ], and negative Ω arag /[O 2 ] triple compound extreme events are suppressed while the intensity and occurrence of marine heatwaves on the shelf seafloor tends to increase.Assuming GOADI can be represented by a first-order autoregressive [AR(1)] process, we observe a relatively short memory of the index with an inverse damping rate λ −1 of The more intense and frequent triple events were mainly triggered by the secular trend of ocean acidification, as there were less positive [H + ] (area increased by 1% per year) and negative Ω arag (0.85%/year) extreme events at the beginning of the timeseries.Earlier model results have shown that surface Ω arag has statistically significantly decreased by −0.064 ± 0.002 decade −1 within 100 km of the coastline, and faster in the center of the gyre (−0.080 ± 0.002, Hauri et al., 2021).While the NGAO is responsible for the enhanced ocean acidification rate in the gyre, shelf waters followed the trend expected from increasing atmospheric pCO 2 .We did not find any evidence that the intensity and frequency of deep-water intrusion has changed over the years.negative Ω arag , purple: negative oxygen, positive [H + ], heat), and quadruple compound extreme events (negative oxygen, positive [H + ], heat, negative Ω arag ).The color bars indicate the intensities of the marine heat, acidification, and low oxygen extreme events in red, blue, purple, and green, respectively.The benthic area on the shelf is defined as the area with a bottom depth between 50 and 250 m depth.The intensity of the events was normalized to the maximum occurring intensity.
The interaction between marine heat waves (remotely triggered) and the local natural variability of deep-water intrusion with anthropogenically amplified CO 2 triggered heat, positive [H + ], negative Ω arag , and negative [O 2 ] quadruple compound extreme events on the shelf seafloor.While both marine heat waves (2015-2016 and 2019-2020) affected >80% of the seafloor on the continental shelf, widespread (26% of shelf area over an average of 40 (±32) days) events with compound marine heat, positive [H + ], and negative Ω arag /[O 2 ] quadruple extremes only occurred during the second 2018-2020 marine heat wave (Figure 10f).During this time, the downwelling index was generally negative, and triggered intrusion of deeper, colder, CO 2 -rich, and O 2 -poor water onto the shelf.While the colder water dampened the intensity of the bottom water heat wave, it was not enough to avoid the occurrence of quadruple compound extreme events across >20% of the shelf area.

What Does This Mean for the Gulf of Alaska Marine Ecosystem?
The Gulf of Alaska marine community experienced dramatic shifts during 2014-2016 and 2019 period, with continued and/or lagged impacts in the ecosystem that persisted for multiple years (Suryan et al., 2021).The impact on commercially valuable groundfish varied by species.Gulf of Alaska Pacific cod (Gadus microcephalus Tilesius) experienced a 71% decline in abundance in 2017 and a commercial fishery closure in 2020 (Barbeaux et al., 2020).Conversely, Alaska sablefish (Anoplopoma fimbria) had strong year classes from 2014 to 2019, and 2016 was potentially the largest recruitment recorded (Goethel et al., 2022).Other changes throughout the ecosystem included reduced primary production, changes in zooplankton community composition, reduced forage fish abundance during the 2014-2016 period followed by a strong and persistent herring year class in 2016 (Arimitsu et al., 2021;Strom & the Northern Gulf of Alaska Long-Term Ecosystem Research Team, 2023), and increased mortality and reduced reproductive success of seabirds and marine mammals (Hastings et al., 2023;Piatt et al., 2020).
The ecological changes during and after this period have primarily been attributed to prolonged, elevated ocean temperatures, yet the consideration of O 2 and pH may also explain ecosystem changes not yet identified or understood.The warm ocean temperatures are suggested to have affected Pacific cod by exceeding optimal thermal thresholds of cod egg and larval survival (Laurel et al., 2023), and increasing metabolic rates, resulting in increased O 2 demand and consumption rates in a reduced prey environment (Barbeaux et al., 2020).However, the cod population has not yet rebounded, despite the return of pre-marine heatwave ocean temperatures and prey abundance, suggesting missing pieces to the narrative.During this time period, some groundfish species shifted to deeper (and therefore cooler) habitats, while others moved to shallower depths (presumably warmer), suggesting the presence of other environmental drivers and restrictions (Yang et al., 2019).Rockfish living deeper along the Gulf of Alaska slope, such as adult thornyhead rockfish (Sebastolobus spp.; 100-1,200 m), rougheye (S. aleutianus) and blackspotted (S. melanostictus) rockfish (300-500 m), and shortraker rockfish (S. borealis; 300-400 m), live in lower oxygen environments and are predicted to potentially move shallower (perhaps into warmer temperatures) if oxygen concentrations decrease further (Thompson et al., 2023).Lower pH could result in a decline of biomass or condition of deep-water corals, important habitat for juvenile Pacific Ocean perch (S. alutus) and numerous other commercially important rockfish.Increased acidity can also be detrimental to pteropods, a calcifying organism and common prey for pink salmon.Pink salmon are predicted to decline in growth and population size if a reduction in pteropod biomass occurs (Aydin et al., 2005).

Usefulness of Our Statistical Approach
The metrics applied here and in previous extreme event studies (Barkhordarian et al., 2022;Gruber et al., 2021;Laufkötter et al., 2020;Le Grix et al., 2021) should be seen as a "probability of trouble" rather than a true indicator of stress.Stress can be defined as a condition evoked in an organism by one or more environmental factors that bring the organism near or over the limit of its ecological niche (Van Straalen, 2003).The biological consequence of a particular stress response depends on its intensity and duration (Boyd et al., 2018).A threshold is the limit beyond which stress and detrimental effects are expected for any environmental parameters such as temperature, carbonate chemistry, or oxygen concentration (IOC-UNESCO, 2022).However, based on our understanding of the biological response to environmental changes, no simple or absolute thresholds are expected.First, these will vary depending on locations (e.g., local adaptation; Vargas et al., 2017Vargas et al., , 2022)).A threshold for a given parameter will also be modulated by other environmental conditions.For example, the pH threshold for sea urchin larvae is strongly influenced by temperature and a greater tolerance to low pH is observed at low temperature (Gianguzza et al., 2014).Complex threshold landscapes emerge because no common thresholds exist, and due to local adaptation and combined effects between environmental parameters (Boyd et al., 2018).Relative thresholds allow accounting for adaptation to the local range of natural variability of a given variable and assume that the organism is jeopardized by variation outside of this natural range.Relative thresholds have been frequently used to define marine heat waves (Barkhordarian et al., 2022;Gruber et al., 2021;Laufkötter et al., 2020;Le Grix et al., 2021) as well as ocean acidity extreme events (Burger et al., 2020;Desmet et al., 2023;Gruber et al., 2021).While relatively easy to calculate and based on a good theoretical framework, this approach does not fully integrate the complexities driving local species and ecosystem sensitivities, and often ignores the physiological, ecological and evolutionary consequences of subsequent or compound extreme events (Gruber et al., 2021).Consequently, we are still lacking the level of understanding to be sure that the extreme and compound extreme events computed here and in previous studies are leading to a true stress.In other words, we cannot assume that the calculated "index" can correlate to a measure of stress.Nevertheless, mapping the spatial and temporal occurrence of environmental extremes and understanding their drivers is an important exercise.First, based on what we know, it is likely that extreme events are leading to stress and when it comes to stressors, two stressors are always worse than one (by definition).So, while we cannot infer the biological response, the probability of negative effects on marine life is there.Second, knowing where, when, and why compound extreme events occur should inform mitigation strategies and improve the predictability of such events.
Translating extremes in environmental conditions to a biological response would require more than a correlation between physical-chemical and biological observations.It would require a strong dialog between different disciplines such as modeling, physical, biogeochemical and biological observations, experimental biology, physiology, ecology, and multiple stressors.For example, laboratory experiments should focus on mechanistic understanding of the impact of key drivers to resolve their mode of action and performance curves under extreme conditions (e.g., heat, positive [H + ], negative [O 2 ], and negative Ω arag stress), the role of intensity/duration, and the relative contribution of each parameter on the response.For example, temperature, pH and [O 2 ] are impacting the physiology and metabolism of marine organisms and can lead to mortality under extreme conditions.However, it is unclear what the characteristic of a marine heat extreme event is that would drive a similar response as from a given acidity extreme event.Such a mechanistic understanding at different levels of organization would allow one to understand and model their combined effects.This could lead to a "formula" that would translate the physical-chemical characteristics of the extreme and compound extreme events into a "stress index" and the biological consequences.Combining this mechanistic understanding with modeling could ultimately explain the observed ecological changes (IOC-UNESCO, 2022) and project future changes.

Summary and Conclusions
The Gulf of Alaska has been exposed to several marine heat waves, which were associated with devastating consequences for the marine ecosystem.Our modeling results suggest that ocean acidity and low oxygen extreme events have co-occurred and thereby may have played a role in the observed ecosystem effects.
A mechanistic understanding of compound extreme events and their consequences is necessary to implement local mitigation and adaptation strategies.Our results show that local climate variability in combination with remote teleconnections and secular trends of ocean acidification and warming led to the manifestation of these compound extreme events.We are now able to use the readily available indices NGAO and GOADI to better describe past and future environmental conditions, which provides the opportunity to directly study what effects these combined extreme environmental conditions may have on the ecosystem.
Prediction and projection of environmental conditions and ecosystem change are a necessary next step to support adaptation and fisheries management decisions in the Gulf of Alaska marine ecosystem.Seasonal forecasts would support quota management on an annual basis, whereas regional long-term projections would provide insights into how local and remote drivers of extreme events will manifest themselves under a changing climate.Our study covers the Gulf of Alaska marine environment, which is in the traditional and contemporary unceded homelands of the Haida, Tsimshian, Tlingit, Eyak, Dena'ina, Sugpiaq/ Alutiiq, and Unangax/Aleut Peoples.Moreover, the offices of the University of Alaska Fairbanks are located on the unceded Native lands of the Lower Tanana Dena.We are grateful to the Indigenous communities, who have been in deep connection with their land and water for now and time immemorial, for the stewardship of their environment.We also recognize the historical and ongoing legacies of colonialism and are committed to improve equity in our lives and scientific institutions we work with.CH and RP acknowledge support from the North Pacific Research Board (NPRB 2109) and National Science Foundation (OIA-1757348, OCE-1656070).SCD acknowledges support from Palmer LTER NSF OPP-2224611.MFS was supported by NSF Grant AGS-2141728.This study has been conducted using E.U.Copernicus Marine Service Information; https://doi.org/10.48670/moi-00021,https://doi.org/10.48670/moi-00016,https://doi.org/10.48670/moi-00148.We would also like to thank two anonymous reviewers and the editor for their valuable input.This is IPRC publication 1613 and SOEST contribution 11754.

Figure 1 .
Figure 1.Northern Gulf of Alaska Oscillation Index.Monthly Northern Gulf of Alaska Oscillation (NGAO) index obtained with (a) modeled and (b) satellite-based (Global Ocean Gridded L 4 Sea Surface Heights & Derived Variables Reprocessed 1993 Ongoing, 2023) sea surface height (SSH).The NGAO index is defined as the first empirical orthogonal function (EOF) and principal component of SSH variability in our model domain(Hauri et al., 2021).The amount of variance associated with the first EOF of the model is ∼24%.

Figure 2 .
Figure 2. Gulf of Alaska regional biogeochemical model domain.De-seasonalized climatology of the modeled sea surface height (SSH, m) is shown in color.Regions obtained from modeled SSH with a 4 clusters configuration using a classical K-MEANS algorithm(Sculley, 2010) are illustrated by lines.The K-MEANS algorithm was carried out with the scikit-learn Python package(Pedregosa et al., 2011).Our current analysis focuses on the gyre (red line) and shelf (magenta line) regions.

Figure 3 .
Figure 3. Evaluation of sea surface temperature across the Gulf of Alaska model domain.Maps show the average (a) modeled (left) and satellite observed (right) sea surface temperature for the time period between 1993 and 2021.Density of probability for modeled (blue) and satellite observed (red) sea surface temperature is shown in panel (b).Modeled (blue) and satellite observed (red) time-series for the regional average are illustrated in panel (c).

Figure 4 .
Figure 4. Evaluation of model-simulated dissolved inorganic carbon, total alkalinity and oxygen.Map of stations occupied during several cruises between 2015 and 2020 overlain on water depth (a).Data from the National Oceanographic and Atmospheric Administration ocean acidification cruise in July of 2015 (Cross et al., 2019) are shown by black triangles.The P16N cruise (May-June 2015, Wanninkhof & Denis, 2016) data are illustrated with blue circles.The Long Term Ecological Research (LTER 2018-2020, Hauri & Irving, 2021) cruises data are shown with pentagons of different colors showing different seasons and years.Taylor diagrams (Taylor, 2001) showing the normalized standard deviation of the model output, the spatial correlation, and the normalized root mean square error between vertical profiles from in situ data and model output for the upper 250 m for (b) dissolved inorganic carbon, (c) total alkalinity, and (d) oxygen.

Figure 6 .
Figure 6.Timing of sea surface marine heat and acidity extreme and compound extreme events over the shelf area.Area affected by (%) of (a) marine heat (red), (b) positive [H + ] (blue), (c) negative Ω arag (purple), (d) positive [H + ] at fixed seasonal temperature (H + Tfix ), (green) extreme events, and (e) positive [H + ] and marine heat (orange), positive [H + ] at fixed temperature and marine heat (black), negative Ω arag and marine heat (purple) compound extreme events.The color bars indicate the intensities of the events.The intensity of the events was normalized to the maximum occurring intensity.

Figure 7 .
Figure 7. Timing of sea surface marine heat and acidity extreme and compound extreme events over the gyre area.Area affected by (%) of (a) marine heat (red), (b) positive [H + ] (blue), (c) negative Ω arag (purple), (d) positive [H + ] at fixed temperature ((H + Tfix ), green) extreme events, and (e) positive [H + ] and marine heat (orange), positive [H + ] at fixed temperature and marine heat (black), negative Ω arag and marine heat (purple) compound extreme events.The color bars indicate the intensities of the events.The intensity of the events was normalized to the maximum occurring intensity.

Figure 8 .
Figure 8.The role of internal climate variability in the occurrence of extreme and compound extreme marine heat and acidity events.Frequency change in event days during negative (left) and positive (right) phases of the Northern Gulf of Alaska Oscillation (NGAO, Hauri et al., 2021) compared to their frequency during a neutral phase for marine heat (a), positive [H + ] (b), negative Ω arag (b), positive [H + ] with fixed temperature (H + Tfix ) (d) extreme events, and for heat and positive [H + ] compound extreme events (e).White cells are not statistically significant (significance level set as 95%).

Figure 9 .
Figure 9. Impact of the Northern Gulf of Alaska Oscillation on summertime temperature.Sea surface temperature (°C) anomaly in summer (June-August) during negative (left) and positive (right) Northern Gulf of Alaska (NGAO) phases between 1993 and 2021.The negative and positive phases of the NGAO index are associated with days when the index was above 50% of its maximum or below 50% of its minimum.

Figure 10 .
Figure10.Timing of marine heat, acidification, and low oxygen extreme and compound extreme events at the shelf seafloor.Area (%) affected by (a) marine heat (red), (b) positive [H + ] (blue), (c) negative Ω arag (blue), (d) negative oxygen (green), (e) positive salinity (orange), and (f) triple (orange: negative oxygen, positive [H + ], negative Ω arag , purple: negative oxygen, positive [H + ], heat), and quadruple compound extreme events (negative oxygen, positive [H + ], heat, negative Ω arag ).The color bars indicate the intensities of the marine heat, acidification, and low oxygen extreme events in red, blue, purple, and green, respectively.The benthic area on the shelf is defined as the area with a bottom depth between 50 and 250 m depth.The intensity of the events was normalized to the maximum occurring intensity.

Figure 11 .
Figure 11.Gulf of Alaska Downwelling Index (GOADI) as indicator for environmental conditions on shelf seafloor.Time series (a) of the normalized principal component associated with the second Empirical Orthogonal Function (EOF) mode of sea surface height.Maps of the EOF (b) spatial patterns of the second mode of sea surface height (SSH, m).The amount of variance associated with the second EOF is ∼9.5%.The EOFs were computed based on the illustrated spatial domain and were applied to daily model output after first removing a long-term temporal trend using a quadratic function and second de-seasonalizing the data.Spatially averaging the coefficients of determination obtained from linear regressions between the normalized GOADI and the grid-point SSH anomalies shows that on average 65% of shelf SSH variance are explained by GOADI.

Figure 12 .
Figure 12.The role of internal climate variability in the occurrence of extreme and compound extreme events on shelf seafloor.Frequency change in event days during negative (left) and positive (right) phases of the Gulf of Alaska Downwelling Index (GOADI) compared to their frequency during a neutral phase for marine (a) heat extreme events, (b) positive [H + ] extreme events, (c) negative Ω arag extreme events, (d) negative [O 2 ] extreme events, (e) positive salinity extreme events, and (f) positive [H + ], and negative Ω arag /[O 2 ] triple compound extreme events.White cells are not statistically significant (significance level set as 95%).