Modulation of ocean acidification by decadal climate variability in the Gulf of Alaska

Uptake of anthropogenic carbon dioxide from the atmosphere by the surface ocean is leading to global ocean acidification, but regional variations in ocean circulation and mixing can dampen or accelerate apparent acidification rates. Here we use a regional ocean model simulation for the years 1980 to 2013 and observational data to investigate how ocean fluctuations impact acidification rates in surface waters of the Gulf of Alaska. We find that large-scale atmospheric forcing influenced local winds and upwelling strength, which in turn affected ocean acidification rate. Specifically, variability in local wind stress curl depressed sea surface height in the subpolar gyre over decade-long intervals, which increased upwelling of nitrate- and dissolved inorganic carbon-rich waters and enhanced apparent ocean acidification rates. We define this sea surface height variability as the Northern Gulf of Alaska Oscillation and suggest that it can cause extreme acidification events that are detrimental to ecosystem health and fisheries. Decadal wind patterns enhance upwelling of carbon-rich waters and increase apparent surface ocean acidification rates in the Gulf of Alaska, according to a 34-year regional ocean biogeochemical hindcast model simulation.

O cean acidification is a slow, multi-decadal, and global environmental press that is primarily driven by the oceanic uptake of anthropogenic CO 2 from the atmosphere. On shorter temporal and spatial scales, ocean acidification interacts with pulse disturbances from the dynamic, coupled ocean-atmosphere system 1 . For example, while it is commonly assumed that the surface ocean CO 2 partial pressure (pCO 2 ) trend approximately follows the atmospheric pCO 2 trend, this assumption may break down temporarily due to internal Earth System variability 2,3 . In addition, localized natural or humandriven processes, such as remineralization of organic matter, upwelling of CO 2 -rich waters, input of low alkalinity freshwater from rivers, glaciers, or sea ice, and input of nutrients and carbon from land sources can contribute to acidification across much shorter timescales [4][5][6][7] .
Ocean acidification involves the rising concentrations of CO 2 in seawater and a decrease in pH and the carbonate ion (CO 3 2− ) concentration, which in turn leads to a decrease in the saturation state (Ω) of biogenic CaCO 3 minerals, such as aragonite (Ω arag ) and calcite. These CaCO 3 minerals chemically dissolve when Ω <1. However, many marine organisms that make shells from CaCO 3 (calcifiers) are sensitive to a decrease in Ω well above the thermodynamic threshold of Ω =1 8 . For this reason, ocean acidification negatively affects important socio-economic subsistence and commercial fisheries, particularly in subpolar regions such as the Gulf of Alaska 9,10 . The Gulf of Alaska commercial fisheries are of economic importance to Canadian and US harvesters and have a high nutritional value for subsistence fishing communities, including the Native Alaskan and Canadian populations [10][11][12] . Salmon, crab, mollusks, and other shellfish make up a substantial part of the highly productive subsistence and commercial fisheries in the Gulf of Alaska and have been shown to be vulnerable to the effects of ocean acidification [13][14][15] .
In the Gulf of Alaska, surface waters are naturally low in CO 3 2− and Ω, compared to tropical and temperate surface waters, due to the increased solubility of CO 2 at low temperatures, ocean mixing patterns, substantial riverine and glacial acidic freshwater inputs, and iron-limited biological carbon drawdown in summer 6,[16][17][18][19] . Indeed, Earth System Models project that under the high emissions scenario (Representative Concentration Pathway (RCP) 8.5), surface waters across large areas of the subpolar North Pacific will become undersaturated with respect to aragonite year-round by the end of this century 20 . The physics, chemistry, and biology of the Gulf of Alaska ecosystem are also strongly influenced by decadal-scale fluctuations of ocean conditions such as characterized by the Pacific Decadal Oscillation (PDO) and the North Pacific Gyre Oscillation (NPGO) indices [21][22][23] . For example, the NPGO closely correlates with interannual fluctuations in ocean salinity, nutrients, Chlorophylla, and O 2 in the Northeastern Pacific 23 and could therefore also affect the inorganic carbon chemistry.
Here, we illustrate how climate-driven regional ocean fluctuations can modulate ocean acidification rates on sub-decadal to decadal timescales with model output from a 34 year (1980−2013) long regional ocean biogeochemical hindcast simulation for the Gulf of Alaska 24 . These fluctuations can exacerbate the seasonal periods of stress from ocean acidification and lead to ocean acidification extreme events earlier than if ocean acidification were the only driver.

Results and discussion
Following previous work 23 , a statistical Empirical Orthogonal Function (EOF) decomposition was performed on monthly anomalies (trends and monthly climatology removed) for SSH, salinity, temperature, pCO 2 , Ω arag , and pH in the Gulf of Alaska ( Fig. 1) to determine the leading mode of variability in this system. One of our key physical variables of investigation is sea surface height (SSH), as it dynamically connects atmospheric drivers with their oceanic responses. Climatic patterns of atmospheric winds drive divergent Ekman transport in surface waters. Variations in wind stress curl and corresponding Ekman suction create areas of low and high SSH, which are characterized by regional upwelling and downwelling, respectively 25 . We define the first EOF/PC of SSH variability in our model domain as the Northern Gulf of Alaska Oscillation (NGAO). The NGAO differs from the NPGO because this new index solely focuses on the northern Gulf of Alaska, whereas the NPGO describes oceanographic patterns of variability across the entire North Pacific. We computed the NGAO for both the satellite and simulated SSH anomalies and found they were significantly correlated with each other (r = 0.60, p < 0.05, Fig. 2b), suggesting that the model captures about 36% of the monthly anomalies in the observed SSH variability. This result suggests that our model is capable of translating large-scale atmospheric forcing into corresponding variations in SSH and, therefore, variability in the upwelling strength of the subpolar gyre 26,27 . The NGAO index (Supplementary Data 1) captured 24 and 28% of the satellite and simulated SSH anomaly variance, respectively, and the strongest SSH anomaly spatial signals were located in the subpolar gyre ( Fig. 1a, b). The maximum modeled SSH variations associated with the NGAO were 0.06 m.
The NGAO index (Supplementary Data 1) of modeled SSH anomalies was significantly correlated (p < 0.05) with the principal components of the 1st EOFs of modeled salinity, temperature, NO 3 , pCO 2 , Ω arag , and pH anomalies at the surface and at 100 m depth, which is below the modeled euphotic zone ( Fig. 2 and Supplementary Fig. 1). The subpolar gyre region displayed again the strongest spatial signal (Fig. 1c-h and Supplementary Fig. 2). The strong correlation between those variables and the NGAO index (Supplementary Data 1) of modeled SSH anomalies ( Fig. 2 and Supplementary Fig. 1) reflects the direct link between simulated variations in subpolar upwelling and ocean surface and subsurface physical and biogeochemical variables. As such, anomalously low sea level pressure over the Gulf of Alaska leads to anomalously strong wind stress curl, thereby depressing SSH via Ekman suction ( Supplementary Fig. 3) and enhancing the upwelling of saltier, colder, and more carbon-rich water to the surface. The maximum pCO 2 , Ω arag , and pH variations associated with the 1st normalized EOF are stronger in the upper thermocline at 100 m depth (pCO 2 = 129, Ω arag = 0.12, pH = 0.06), damped toward the surface (pCO 2 = 15, Ω arag = 0.08, pH = 0.01), and account for between 15 and 22% of the total variance of the detrended and deseasonalized model fields at 100 m depth ( Fig. 1f-h and Supplementary Fig. 2d-f).
The sub-decadal pulse disturbances to Gulf of Alaska inorganic carbon system fields associated with the NGAO index (Supplementary Data 1) overlay onto the longer-term ocean acidification trends caused by ocean uptake of anthropogenic CO 2 . Over the 34year simulation (1980−2013), modeled ocean surface Ω arag and pH decreased and pCO 2 increased over the Gulf of Alaska with mean regional trends of Ω arag = −0.071 +/− 0.003 decade −1 (1 STD), pH = −0.019 +/− 0.001 decade −1 , and pCO 2 = 17.4 +/− 1.4 μatm decade −1 (p < 0.05, Supplementary Fig. 4a-c). Similar patterns of interannual climate variations and secular trends in the ocean carbon system are also found in field observations in the Northeast Pacific, south of our model domain 28 . However, the modeled trends over this time period exhibited large regional variations with stronger acidification signals in the subpolar gyre than in the Gulf of Alaska coastal area. The mean central subpolar gyre (55°N   The EOFs were computed based on the illustrated spatial domain and were applied to monthly model output after first removing a long-term temporal trend using a quadratic function and second deseasonalizing the data. A 100 km wide band along the coast has been masked out for surface salinity to avoid the interference of freshwater from the coastal river and glacial discharge with the EOF analysis.  It is important to note that the subpolar oceanic pCO 2 trend is well above the atmospheric pCO 2 trend of 18.7 μatm decade −1 . In part, the regional variations in the multi-decadal trends may reflect the simulation time window, which started in a period of a neutral or positive NGAO phase and ended following a period of a negative NGAO phase (Fig. 2a, b). Ending with a negative NGAO could bias the trend estimates towards larger absolute values. As such, if the time window ended with a more neutral NGAO, a weaker slope would be expected. Wind stress curl from the Japanese 55-year Reanalysis (JRA55-do) 1.3 project 29 , which was used to force our hindcast simulation, strengthened between 1980 and 2013 in the western part of the Gulf of Alaska (Supplementary Fig. 5). The enhanced wind stress curl intensified the model subpolar gyre circulation and resulted in stronger upwelling via Ekman suction dynamics 30 upwelling colder, saltier and alkalinity-, nitrate-and carbon-rich waters to the ocean surface and accelerating the rate of ocean acidification in the subpolar gyre ( Supplementary Figs. 6 and 7) compared to what would be expected from anthropogenic air-sea carbon fluxes alone.
The combination of the long-term acidification trends and subdecadal climate variability-driven pulses can push the ecosystem across biologically relevant thresholds of Ω arag and pH, and pCO 2 and cause extreme events earlier than if ocean acidification was the only driver. During negative NGAO phases in 2007−2008 and 2011−2013 (Fig. 3a), simulated surface Ω arag decreased to levels below the biologically relevant threshold for pteropods (Ω arag = 1.3) 31 , whereas, during positive NGAO phases, surface water conditions remained mostly favorable for pteropods. Similarly, by choosing an illustrative threshold for pH of 8.0 and pCO 2 of 450 μatm, it becomes apparent that extreme events occurred about 8 years earlier (pink dots in Fig. 3a, b) due to the climate-driven sub-decadal pulses to the system than if the Gulf of Alaska had experienced only the more gradual press of ocean acidification alone (solid black line in Fig. 3a, b). Figure 3 also illustrates that the correlation is not perfect. As the NGAO index (Supplementary Data 1) declines in 1987 and 1988, pCO 2 becomes, as expected, anomalously high. However, this biogeochemical signal does not continue as the NGAO minimum deepens in 1989 and 1990, for reasons not yet understood.
In particular, the 2007−2008 and 2011−2013 ocean acidification extreme events in the subpolar gyre, which were facilitated by a negative NGAO pulse superimposed on the underlying press of ocean acidification (Figs. 2 and 3), may have already had a direct impact on the ecosystem. Pteropods collected from the subpolar gyre in 2015 show severe shell dissolution (personal communication Nina Bednaršek 32 ). This is concerning since these calcifying mollusks are an indicator species for habitat suitability related to ocean acidification, and comprise important components of the diet of various salmon species at different life stages [33][34][35] . While negative phases of the NGAO may have always led to extreme events, ocean acidification has pushed the system closer to chemical thresholds that are harmful to certain organisms. As a consequence, ocean acidification may cause these extreme events in ocean chemistry to become more frequent, Fig. 3 Occurrences of ocean acidification extreme events. Plots showing the relative occurrence of a saturation state for aragonite (Ω arag ), b pH, and c pCO 2 (μatm) for an area in the middle of the subpolar gyre (55°N to 57°N, 147°W to 150°W) from 1980 through 2013. The blue line shows the Northern Gulf of Alaska Oscillation (NGAO) index at a monthly frequency, smoothed with a 12 month running mean to highlight the effect of sub-decadal pulse disturbances on top of the decadal press from rising atmospheric CO 2 . The dashed blue line separates the positive from the negative NGAO phases. The solid black line shows the 10th percentile of the distribution of Ω arag and pH, and the 90th percentile of the distribution of pCO 2 after applying a 10-year low pass filter to illustrate when an arbitrary threshold would be crossed as a result of the long-term press of ocean acidification. The pink dots show the 10th percentile of monthly Ω arag and pH, and the 90th percentile of monthly pCO 2 model output to indicate when an arbitrary threshold would be crossed as a consequence of the NGAO (or pulses) in addition to seasonal periods of ocean acidification stress. The dashed red line is an arbitrary threshold to mark extreme events, when the sub-decadal pulse disturbances push the system across these thresholds (Ω arag = 1.3; pH = 8.0; pCO 2 = 450 μatm).
intense, and longer in future years, thereby putting an additional strain on the marine ecosystem.
Our results show that climate-driven sub-decadal scale natural variability greatly influences the rate of change of the inorganic carbon system over a given time period and could either amplify, dampen, or even temporarily reverse the apparent rate of ocean acidification. For example, if an observational time series effort began in the middle of a negative NGAO phase and ended at the peak of the next positive NGAO phase (e.g., 1988−1993), then the sub-decadal variability would counteract the underlying acidification during the observational record, potentially even leading to a temporary positive trend in pH and Ω arag (Supplementary Fig. 4). The satellite SSH data were available through 2019, showing that the strong negative NGAO phase in the early 2000s was followed by a positive phase until the present day. It is very likely that if our model simulation covered the period between 1980 and 2020, the ocean acidification trend over the full record would be slowed in the subpolar gyre and accelerated on the shelf relative to the trend we found for the 1980−2013 time period. This is illustrated in Supplementary Fig. 4, which shows an apparent trend across a shorter time period starting and ending with a positive phase of the NGAO (1980NGAO ( −2005. This finding demonstrates that in some regions of the Gulf of Alaska, in particular the subpolar gyre, even 34 years of high-resolution climate-quality observational data would be insufficient to accurately quantify the rate of anthropogenic ocean acidification. Ocean acidification will also interact with other anthropogenic stressors affecting the Gulf of Alaska marine ecosystem, such as marine heatwaves that are predicted to become more frequent and will eventually keep constant pressure on the ecosystem or even overlap in time and space. For example, the 2014−2016 marine heatwave in the Northeast Pacific, also known as 'blob', is the globally largest and longest marine heat wave to date 36 , and had severe consequences on the marine ecosystem and dependent economies [37][38][39][40] . This marine heatwave resulted from a combination of a press (climate change) and specific pulses linked to varying climate modes (such as NPGO, PDO, and El Niño) 36,41 and started right at the end of our hindcast simulation and therefore directly followed the ocean acidification extreme event. The resilience of the marine ecosystem may depend on the ability of the ecosystem to absorb interruptions to trophic linkages. As anthropogenic stressors move the climate system into previously unseen states, the ecological consequences of multiple simultaneous stressors are likely manifested in nonlinear and surprising ways. Ecosystem changes will unfold over biological time scales that can decouple from the time scales of the pulse forcing 40 . Additional emphasis needs to be put on how marine heatwaves and ocean acidification extreme events affect single species and the ecosystem dynamics as a whole.

Conclusions
Going forward, hindcast simulations such as the one presented here could prove useful to add historical context to observed shorter biogeochemical and biological time series and specific extreme events. This study shows that in certain regions, such as the greater Gulf of Alaska and subpolar gyre, natural (sub-)decadal variability can strongly dampen or accelerate the apparent rate of ocean acidification. Improving our understanding of the influence of climate teleconnections on the biogeochemistry across the global ocean will help us with the interpretation of disturbance events and changes to the ecosystem. In addition, such simulations could be used to determine regions with large interannual and decadal-scale variability, where chemical and biological time series could be established to more efficiently detect the early crossing of biological and ecological thresholds.

Methods
The hindcast simulation (1980−2013) 24 was conducted with a 4.5 km resolution Gulf of Alaska setup of the three-dimensional Regional Oceanic Modeling System (ROMS, 42 ) coupled to the marine ecosystem model 3PS-Ocean Biogeochemistry and Lower Trophic (3PS-COBALT, 43,44 ). The Simple Ocean Data Assimilation ocean/sea-ice reanalysis 3.3.1 (SODA 45 , available as 5-day averages) provided physical initial and boundary conditions for ocean temperature, salinity, currents, and sea surface height (SSH). Surface forcing, including 3 hourly winds, surface air temperature, pressure, humidity, precipitation, and radiation was taken from the Japanese 55-year Reanalysis (JRA55-do) 1.3 project 29 . Freshwater discharge from numerous rivers and tidewater glaciers was implemented as forcing from a hydrological hindcast simulation 46 and was brought into the coastal domain as point sources via the coastal wall at all depths 47 . A spatially uniform climatology of riverine nutrients, DIC, and TA was based on available observations 24,48,49 . For initial and boundary conditions, DIC and TA from the Global Ocean Data Analysis Project (GLODAPv2.2016b, 50 ) data set were normalized to 1980 and adjusted to each following year using anthropogenic CO 2 estimates for the Gulf of Alaska 51 . The hindcast simulation is based on a 10-year spin-up of the 1980 conditions. A 30 year-long control simulation (perpetual 1980 forcing) was also conducted, which did not show any statistically significant trend in the variables discussed here. The 1980−2013 GOA-COBALT hindcast simulation was quantitatively evaluated with available satellite Chl-a and coastal in situ hydrographic observations, including inorganic carbon system parameters 24 . Overall, the model slightly underestimates nearshore DIC drawdown by primary production in spring and overestimates Chla concentrations across the domain. Offshore, the model simulates observed DIC and TA reasonably well ( Supplementary Figs. 8 and 9). The reader is referred to ref. 24 for more information on the simulation and model evaluation.
EOF analysis. We used EOF analysis to determine the dominant modes of variability by decomposing the oceanic fields into a set of independent and uncorrelated spatial modes (EOFs) and their corresponding temporal fluctuations or principal components (PCs). This is done by computing the eigenvectors of the covariance matrix 52 . EOF analysis was applied to monthly model output after first removing a long-term temporal trend using a second-order polynomial of the form z*x 2 + b*x + c and second deseasonalizing the data by subtracting the monthly climatology from the detrended data set. A second-order polynomial was used because not all variables showed a linear trend. The EOF analysis was conducted with the Climate Data Toolbox for MATLAB 53 . The EOF principal component time series were normalized by the PC's respective standard deviation and the spatial maps are displayed in physical units. The leading EOFs are well separated from the next modes based on their respective eigenvalues (SI Table 1). However, for some variables, such as pCO 2 and pH, the 1st EOF explains only a small part of the total variance.
Pearson correlations. To calculate the Pearson correlations between the first Principal Components of the different variables, we accounted for autocorrelation, and determined the effective number of degrees of freedom following the methodology described in ref. 54 .