Climate attribution time series track the evolution of human influence on North Pacific sea surface temperature

We use the Fraction of Attributable Risk (FAR; [18–20]) to attribute the human contribution to observed SST anomalies:

$$\begin{align*}{\text{FAR}} = 1 - {\text{ }}\frac{{{\text{preindustrial probability}}}}{{{\text{current probability}}}}.\end{align*}$$

The preindustrial probability is the likelihood that a given event might occur in a counterfactual world without human influence on the climate, while the current probability is the likelihood for the same event in the factual world with human influences [21]. Using the ratio of these probabilities, FAR gives an estimate of the proportion of risk for a class of event (SST anomalies ⩾ observed) that can be assigned to human influence. We also calculate the related Risk Ratio (RR) statistic [21], the ratio of risk for the same class of events with and without human activities:

$$\begin{align*}{\text{RR}} = {\text{ }}1/\left( {1 - {\text{FAR}}} \right).\end{align*}$$

A RR value of 10, for example, is evidence that a given class of event is ten times more likely to occur due to human influence on the climate.

We calculated FAR and RR values for ERSSTv5 SST observations [22] from 1950–2022 averaged over six areas: the North Pacific basin poleward of 20° N, and five northeast Pacific regional ecosystems (Eastern Bering Sea, Gulf of Alaska, British Columbia Coast, and Northern and Southern California Current; figure 1). Our regional focus is motivated by the scale over which fisheries are typically managed.

Zoom In Zoom Out Reset image size

We calculated FAR and RR from individual ensemble members for 23 Coupled Model Intercomparison Project Phase 6 (CMIP6) models ([23]; table 1), using normalized SST anomalies relative to the pre-1950 time series (1854–1949 for observations, 1850–1949 for CMIP6). This reference period was selected to be long enough to summarize a wide range of conditions related to internal variability, before the 1960s emergence of the observed anthropogenic trend in global temperature [24]. For each SST observation, preindustrial probabilities were estimated for each CMIP6 model as the proportion of SST anomalies from preindustrial runs producing anomalies as large as or greater than the observed SST anomaly. Present-day probabilities were calculated for each model using combined historical and SSP5-8.5 runs. Population dynamics for exploited fishes are affected by climate conditions at multiple life stages [25, 26], so average conditions over several years may be more important for fisheries yield than annual conditions [27]. We therefore calculate FAR for both annual and three-year running mean values of SST.

Model bias and non-independence may make the equal weighting of climate models suboptimal [28], so we weighted CMIP6 models to reward model skill and independence [29, 30]. For each model i, the weight w_i is calculated as:

$$\begin{align*}{w_i} = \frac{{{e^{ - \frac{{D_i^2}}{{\sigma _D^2}}}}}}{{1 + \mathop \sum \nolimits_{j \ne i}^M {e^{ - \frac{{S_{ij}^2}}{{\sigma _S^2}}}}}},\end{align*}$$

where D_i is the distance of model i to the observations (i.e. the difference between model i and observations in scaled units of the weighting variable), S_ij is the distance of model i to model j, M is the total number of models, and σ_D and σ_S are the shape variables; the numerator is the skill weight and the denominator is the independence weight. Values of σ_D and σ_S were selected using a perfect model setup (see appendix for details). We used four weighting criteria: the bias and trend in SST climatology, interannual SST variability, and interannual SST autocorrelation. Selection of these variables was motivated by their status as fundamental physical properties that are important both for FAR calculations and for biological responses. Values for each weighting variable were scaled by the median across all 23 CMIP6 models for each region, so that each variable would carry roughly the same importance. The four variables were then averaged for each model-region combination to calculate D_i , the difference score for each model i. The difference for each model pair i,j (S_i,j) was calculated as the mean of the inter-model differences for the same variables, after each variable was scaled by its median. Multi-model estimates for preindustrial and 'present' conditions were then generated with a binomial Bayesian regression model that weighted probabilities from each CMIP6 model using the normalized weights (see appendix for details). To evaluate the relationship between observations and CMIP6 output using this weighting approach, we plotted ERSST anomalies (in ˚C) with respect to the pre-1950 mean for each region against weighted mean CMIP6 anomalies from the model ensemble for historical, SSP2-4.5, and SSP5-8.5 runs, with the range of most likely outcomes for each CMIP6 run plotted as ± twice the weighted SD.

Expected return time (the inverse of annual probability) is an established approach for evaluating the likelihood of extreme climate events [31]. To evaluate the changing incidence of extreme SST events over time, we calculated the expected return time for the warmest SST observation in each region at a range of North Pacific warming levels (preindustrial, 1950 to 0.5° warming, 0.5° to 1.0° warming, 1.0° to 1.5° warming, 1.5 °C to 2.0 °C warming). The CMIP6 multi-model probability of extreme events was estimated with a weighted Bayesian regression using annual probability (proportion of annual SST anomalies ⩾ the largest observed) as the response variable, the North Pacific warming level as a categorical explanatory variable, and the model weights described above.

Timing of 0.5 °C warming increments for the North Pacific was estimated for observations using a Bayesian smoothed regression fitting year (as the response variable) to warming relative to 1854–1949 (explanatory variable). For CMIP6 models, a weighted regression of the same formulation was employed to generate a multi-model estimate, using a single weighting criterion based on skill and independence for modeling the observed warming trend (see appendix for detailed Methods). Warming rates for CMIP6 were estimated separately for SSP2-4.5 and SSP5-8.5.

We use the Gulf of Alaska commercial sockeye salmon (Oncorhynchus nerka) fishery as a case study to show how attribution time series can contextualize the human role in observed climate change impacts on ecosystem services. Gulf of Alaska sockeye salmon support a commercial fishery with an average ex-vessel value of $64.4M USD during the last decade [32], and are an important cultural, subsistence, and recreational resource. Commercial harvest data were obtained from the Alaska Department of Fish and Game [32]. Productivity in this species (adult offspring per spawner) shows strong responses to ocean temperature at multiple life history stages, such that three-year running mean SST is a better predictor of productivity than annual SST [25, 27].

This analysis followed four steps. First, we confirmed the negative relationship between ocean warming and sockeye productivity [27] with a Bayesian regression model relating catch to three-year running mean SST, weighted by the average year of ocean entry for each harvest year (see appendix for details). We then put this negative consequence of warming into the context of human-induced climate change by fitting a second Bayesian regression model to FAR values corresponding to three-year SST. Then, since FAR values fix at 1 and become uninformative past a certain level of warming, we used a categorical model to compare fishery harvest for salmon entering the ocean under the strongest evidence for anthropogenic impacts (FAR ⩾ 0.98) with ocean entry years with more moderate evidence for anthropogenic impacts (FAR ⩽ 0.91), with the cutoffs for the two groups defined by a bimodal distribution in the data. Finally, we evaluated changing climate risk for this fishery by calculating probability density functions for three-year running mean Gulf of Alaska SST relative to an identified sockeye salmon critical threshold (see Results). These probabilities were derived from CMIP6 historical/SSP5-8.5 runs for different North Pacific warming increments as defined for each model with Bayesian regression. We then defined the probability density by resampling the multi-model population of anomalies from each warming increment, weighting each model based on the weighting criteria described above.

The natural log of catch was used as the response variable and the models also included an autocorrelated error term to account for serial dependence in catch. We limited analysis to catch years 1989–2022 in order to avoid the complexity of nonstationary SST-salmon relationships earlier in the time series [27, 33].

Bayesian models were fit using Stan 2.21.0 [34] and R 4.1.2 [35]. For all parameters, the potential scale reduction factor ($\hat R$) was less than 1.01, effective sample sizes were greater than 500, and there were no divergent transitions in the No-U-Turn-Sampler Markov chain Monte Carlo algorithm. We evaluated chain convergence and model fits graphically (e.g. with trace plots), and via posterior predictive checks [36].

Raw SST for CMIP6 models showed strong and variable bias compared with observations (appendix). However, weighted CMIP6 anomalies showed good agreement with observations; the envelope of plausible historical hindcast values captured the range of observed anomalies as expected, and recent warming is consistent with the weighted projections (figure 2).

Zoom In Zoom Out Reset image size

All five regions, plus the overall North Pacific, have reached temperatures that can be attributed to human activities with high confidence (i.e. FAR > 0.98). In particular, North Pacific SST has effectively become fixed in an anthropogenic state, with FAR ⩾ 0.99 for three-year running mean SST since 2003, and for annual SST since 2013 (figure 3). The last years when regional SST was at levels equally likely with or without human influence (i.e. FAR 95% credible intervals that include zero) are clustered in the years following the well-described 1976/77 shift to a positive state in the Pacific Decadal Oscillation (table 2). All regional ecosystems showed FAR values > 0.98 for at least part of the 2014–2019 marine heatwaves that presented unprecedented climate perturbations to northeast Pacific ecosystems [37–39]. However, the trend and pattern in reaching these high levels of attributable risk varied among regions. Most notably, FAR estimates show a north-south gradient that is consistent with Arctic amplification of global warming. The two northernmost ecosystems have shown FAR > 0.5 (anomalies twice as likely due to human influence) for three-year running mean SST in every year since 1984 (Eastern Bering Sea) and in every year but one since 1993 (Gulf of Alaska). The strongest evidence for anthropogenic influence is seen in the northernmost region (Eastern Bering Sea), which exceeded FAR values of 0.99 for three-year running mean SST as early as 2003/2004. Annual SST FAR values for this system also exceeded 0.99 for each year from 2014 through 2020. Annual SST FAR in the Gulf of Alaska and British Columbia Coast exceeded 0.98 for at least two years since 2015, and these systems have experienced FAR > 0.99 for three-year running mean SST for five or six years since 2014. Recent FAR values were lowest for the lower latitude systems (Northern and Southern California Current), with values > 0.99 for annual SST in only one year (2015), and values > 0.99 for three-year running mean SST during 2014–2017.

Zoom In Zoom Out Reset image size

FAR values cannot exceed one, reducing the sensitivity of this metric at high levels of human influence. RR time series provide additional information about the evolution of human influence on climate to date: annual SST in each region has reached extreme levels that are tens to hundreds of times more likely due to human activities (figure 3). The effects of spatial scale (basin vs. regional) and Arctic amplification are also illustrated by RR time series. North Pacific RR values have exceeded 10⁵ for both annual and three-year mean SST; Eastern Bering Sea RR values have exceeded 10⁴ (for three-year mean SST), while RR values for other regions have peak values of 10²–10³ (figure 4).

Zoom In Zoom Out Reset image size

We found a rapid decrease in expected return times for extreme SST values, as anomalies that were exceedingly rare have become common events (figure 5(A)). In addition, there are important regional differences that provide context for observed SST extremes in the different ecosystems. The Eastern Bering Sea has experienced the largest annual SST anomaly to date of any of the regions we consider (5.08 SD), and projected return times indicate that anomalies this large should be expected every ∼40 years in the current climate (1.0 °C–1.5 °C warming, see below). The neighboring Gulf of Alaska has experienced the smallest observed anomaly of any of the regions (2.98 SD), and anomalies this large are expected every ∼5 years in the current climate. Strong ecosystem and fisheries responses have been ascribed to warm anomalies in the two systems [39, 40], but these results suggest that stakeholders should expect repeat events to be much more common in the Gulf of Alaska than in the Bering Sea.

Zoom In Zoom Out Reset image size

We found good agreement between observations and CMIP6 scenarios for the timing of 0.5 °C and 1.0 °C of warming (range of 2001–2003 and 2018–2021, respectively). The projected timing of 1.5 °C and 2.0 °C begins to diverge for the two SSP scenarios in a way that demonstrates the projected effects of mitigation on North Pacific warming (projected dates for SSP2-4.5 of 2039 and 2060 for 1.5 °C and 2.0 °C, respectively; and 2032 and 2044, respectively, for SSP5-8.5; figure 5(B)).

Regression of sockeye salmon catches on SST confirmed the documented negative response of the fishery to recent warming (figure 6(A)). A companion regression of catch on FAR values indicates that declines in sockeye harvests can also be statistically related to the increasing evidence for human influence on Gulf of Alaska climate (figure 6(B)). FAR values for catch years 1989–2016 varied between 0.33 and 0.91, indicating unmistakable but relatively moderate evidence for anthropogenic influence (SST values that are 1.5–11 times more likely because of human activity). For catch years 2017–2022, evidence for anthropogenic influence on temperature was much stronger (FAR = 0.98–0.99, indicating temperatures that were 60–190 times more likely due to human influence). We found a decline of 1.4SD in ln(catch) for the high-FAR years (FAR ⩾ 0.98), which corresponds to a 39% decline in raw catch relative to the 1989–2022 mean (figure 6(C)). The FAR values in the high-anthropogenic state (FAR ⩾ 0.98) correspond with a critical threshold of three-year running mean SST 2.08 SD above the pre-1950 climatology. Probability densities for three-year mean SST values at different warming increments illustrate the rapid change in the expected incidence of temperatures at or above this threshold value, and thus an estimate of changing climate risk for this fishery. Temperatures at or above the threshold were rare to uncommon events in historical climates (0% of three-year running mean SST values in the preindustrial climate, 1.4% between 1950 and 0.5 °C warming, 20% for 0.5 °C–1.0 °C of warming). But the critical threshold is close to the median value in the current climate (49% of values above the threshold for 1.0 °C–1.5 °C of warming). And the critical threshold is projected to represent cooler-than-normal conditions at 1.5 °C–2.0 °C warming (80% of values at or above the threshold; figure 6(D)).

Zoom In Zoom Out Reset image size

1.1. Calculating attribution statistics

1.2. Model weighting

CMIP6 model outputs for calculating FAR and RR values were weighted with an approach that rewards model skill and independence [29, 30]. For every model i, the weight w_i was calculated as

$$\begin{align*}{w_i} = \frac{{{e^{ - \frac{{D_i^2}}{{\sigma _D^2}}}}}}{{1 + \mathop \sum \nolimits_{j \ne i}^M {e^{ - \frac{{S_{ij}^2}}{{\sigma _S^2}}}}}}\end{align*}$$

where D_i is the distance of model i to the observations, S_ij is the distance of model i to model j, M is the total number of models, and σ_D and σ_S are the shape variables. In this approach the numerator is the skill weight and the denominator is the independence weight.

Model weighting was based on skill and independence for bias in climatology, and variability, autocorrelation, and the trend in SST. In addition to being basic descriptors of temperature variability and change that are important for FAR calculations, each of these weighting criteria also reflects an ecological consideration. Bias is important because nonlinear biological responses to temperature make the actual temperature (i.e. in °C), rather than the anomaly (normalized difference from preindustrial) important. The trend in SST describes the rate of change in the climate. And interannual variability and autocorrelation are important because the red noise in SST is an important driver of fisheries variability. Bias weights were calculated as the absolute value of the difference in mean annual temperature between the model and observations, trend weights as the absolute value of the difference in coefficients from linear regression models of SST on year in CMIP6 models and observations, and standard deviation and autocorrelation weights as the absolute difference between modeled and observed estimates of these quantities.

Bias, variability, and autocorrelation were calculated for 1950–2014 (the period of historical CMIP6 runs that corresponds to the highest-quality observational data), and the trend in SST was calculated for 1973–2022 (a 50 year period capturing most of the observed anthropogenic trend). Weights were calculated separately for the entire North Pacific basin and each of the five regions.

Model weights in this approach are highly dependent on the value of σ_D . Low values of σ_D can result in an overly precise estimate based on only a few models, and large values of σ_D tend towards model democracy. Model weighting results are generally less sensitive to σ_S. We used a perfect model setup to tune the values of σ_D and σ_S to achieve the desired precision in weighted model estimates. In this setup, each of the 23 CMIP6 models is in turn treated as the 'truth'. Values of D_i and S_i,j are then calculated for the remaining models. Weighted predictions are calculated for each of the four weighting variables (bias, trend, standard deviation, autocorrelation) at a range of possible σ_D and σ_S values. The lowest value of σ_D that results in at least 80% of 'true' model values falling between the 10th and 90th quantiles for the weighted prediction is then adopted for use in the study.

Based on the results of the perfect model setup, we adopted a value of σ_S = 0.4 for all regions. Values of σ_D varied between 1.16 and 1.84 for the different regions. All regions reached the desired threshold of 80% of 'true' values falling between the 10th and 90th quantiles of weighted prediction with these combinations of except for the Southern California Current, which reached 79.3% (figure A1). Combined weights (w_i ) calculated using these shape parameters showed a stronger effect of skill weight than independence weight (figure A2). Normalized values of model weights (where a weight of one corresponds to an unweighted model) showed the most tendency towards model democracy for the North Pacific and Southern California Current, and the strongest weighting for individual models in the British Columbia Coast (figure A3). The relationship between SST observations and modeled values by weight is plotted in figure A4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We used a separate round of model weighting for estimating the timing of different North Pacific warming thresholds under SSP2-4.5 and SSP5-8.5 (figure 5(B) in the Main Text). In this instance we only weighted models by the trend in anomalies during 1973–2022 (root mean square error for linear regression of observed anomalies on modeled anomalies). This single weighting variable reflected the goal of this part of the analysis, which was to project North Pacific warming. We used the same perfect model approach for this round of weighting, with the 'true' climatology for two periods (2001–2015 and 2031–2045) as the targets for prediction. Following the results from the previous perfect model approach, we only considered a value of σ_S = 0.4, and the resulting plot indicated selection of a value of σ_D = 0.82 (figure A5). Weights calculated with these two values of the shape parameters resulted in a greater influence of independence skill than was observed in the model weights that were used throughout the rest of the study (figure A6).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

1.3. Sockeye salmon catch analysis