Introduction

Climate variability has been a subject of interest for ecologists primarily because variations in climate often have a strong impact on ecological systems1. Marine resources, such as fish stocks, have been shown to be strongly influenced by climate variability2,3,4, with changes in productivity resulting in huge consequences for socio-economic systems relying on such resources5,6. In the Anthropocene, with the impending threat of climate change, understanding the impact of climate variability on marine ecosystems and resources has become even more central, since climate variability at interannual to decadal timescales can alter the magnitude of ongoing long-term climate change7,8. Hence an integration of climate information into modeling of exploitable resources is necessary not only to understand ecological processes but also to forecast future states of the system at interannual to decadal timescales9. The latter is particularly fundamental for management, since forecasting fish abundance depending on decadal climate variability is necessary to devise timely interventions to ensure sustainable use of resources10. Nevertheless, the application of climate models to predict ecosystem processes at decadal timescales remains a challenge11,12,13.

In many cases the impact of climate on fish stocks has been studied through experiments and modeling, and empirical relations have been established. Climate has been shown to influence fish directly or indirectly through recruitment, food availability, fecundity, growth, and migration14,15,16. Still, climate variables are rarely included in the management-oriented modeling and forecasting of fish populations17. This is not only due to historically large impact of fishing mortality on commercial stock biomass18,19,20 but also due to forecasting being complicated by frequent transient and non-stationary properties of climate impacts on fish stocks21,22,23. Moreover, the fish experience the cumulative impacts of different drivers; fishing pressure and climate can have combined effects inducing non-linear dynamics in fish stocks. Strong synergistic effects can lead to management failures and abrupt collapses of socio-ecological systems22,24,25. With anthropogenic climate change superposed on natural climate variability, including environmental variables in the modeling and forecasting of fish stocks is becoming increasingly important from both scientific and management points of view5,26.

One of the key limitations impeding the integration of climate information in fisheries forecasts originates from inadequate representation of shelf seas in coupled global circulation models (GCMs) providing future climate information. Pioneering-approaches have used bioclimate envelope models27,28 or detailed ecosystem and population dynamics models29,30 forced with climate projections from GCMs to examine the impact of climate change on fisheries. However, GCMs lack a proper representation of shelf-sea dynamics, mainly due to their coarse resolution, and they have limited representation of trophic interactions and associated energy transfers. Other approaches have combined GCM output with highly resolved physical–biological shelf-sea models accounting for trophic interactions6. These approaches focus on long-term (>30 years) changes and thus do not provide information on decadal (1–10 years) fisheries forecasts. Moreover, since decadal forecasting usually involves an ensemble of predictions, high computational costs associated with the aforementioned approaches also motivate exploration of novel approaches towards fisheries predictions using GCM-based decadal climate predictions.

The prospect of decadal prediction of fish stocks emerging from decadal predictability of the physical environment is enticing. Specifically in the North Atlantic, where decadal variability of the physical environment is highly predictable using GCMs31,32,33. This prospect emerges not only from the influence of Atlantic inflow on both hydrography34,35 and marine ecosystem of the North Atlantic shelf seas, such as the North and Barents Sea15,16 but also from the impact of anthropogenic warming on marine ecosystems13. In these climate-driven marine ecosystems, statistical climate–fisheries models36 provide a promising approach for transforming GCM-based ensembles of decadal climate predictions into reliable fisheries forecasts.

In this article we assess predictability of two Atlantic cod (Gadus morhua) stocks in the northeastern North Atlantic shelf seas. Atlantic cod is a commercially, historically, socially, and ecologically important species. There are many stocks of Atlantic cod and they are widely distributed on the shelf seas off the northern North Atlantic. Some of these stocks have been severely reduced in the last decades, largely due to unfavorable climate and intense fishing20,22. Thus, being able to predict the stock biomasses is important to guide sustainable management decisions. We investigate two stocks with opposite status: (1) the North Sea cod, a stock close to the upper temperature limit of distribution of this species, over-exploited for many decades and in a very low productive state over the last 20 years, and (2) the Northeast Arctic cod, residing in the Barents Sea close to the species’ lower temperature limit and recording record-high biomass levels in recent years37,38,39. The average age at first maturation is 3 for North Sea cod and 7 for Northeast Arctic cod.

In order to provide decadal predictions of cod stocks, we use a linear regression model to transform dynamical prediction of sea surface temperature (SST) into the prediction of total stock biomass (TSB). TSB was used because it reflects the integrated impact of climate and fishing36. The dynamical prediction of SST is provided by 10-year long initialized forecasts (and hindcasts) from a decadal prediction system based on the Max Planck Institute Earth System Model (MPI-ESM, see “Methods”). The initial conditions for decadal hindcasts are taken from an assimilation experiment which assimilates observed atmospheric and oceanic information into MPI-ESM. In order to isolate the prediction skill due to external forcing, an ensemble of non-initialized historical simulations of the same size as that of initialized hindcasts and driven by observed external boundary conditions is also analyzed (see “Methods”). Our forecasts for the period 2020–2030 suggest continued unfavorable environmental conditions for the North Sea cod with no significant recovery under any of the three different fishing mortality scenarios. For the Northeast Arctic cod, assuming fishing at the current sustainable level is continued, we forecast a decline in TSB in the coming decade compared to the last decade, attributed mainly to a decline in temperature.

Results

Variability in cod stocks and their physical environment

The time series of SST in the North Sea and Barents Sea Opening highlight key differences in the two regions: While SST has an increasing trend both in the North Sea and in the Barents Sea Opening, the absolute values are very different from each other, highlighting that the cod stocks reside at the two extreme ends of the thermal habitat available for cod40 (Fig. 1a). The TSB time series of the two stocks also show opposite development: North Sea cod has declined continuously since the 1960s, with very low and stable biomass levels since the beginning of the twenty-first century (Fig. 1b). Northeast Arctic cod exhibits multi-annual to decadal variability for the same period, with a recent record high level of TSB (Fig. 1b). However, multi-decadal variability in North Sea stock biomass has been reported for a longer period41. Fishing mortality (F) trends are similar in these two stocks, increasing in the central period of the time series and recently declining as stricter management measures started to be enforced (Fig. 1b). Interestingly, while the decline in fishing mortality of Northeast Arctic cod seems to have resulted in an increase in TSB, in the North Sea, the cod stock did not manage to recover even after the management measures were in place. This has been attributed to the effect of an interacting driver (i.e. warming) which has inhibited the productivity of North Sea cod22.

In the North Sea SST, the magnitude of warming over the period 1960–2019 (1.68 °C) is more than twice the year-to-year variability (σ = 0.65 °C), indicating that the increasing temperature trend is part of how the North Sea has changed under natural and anthropogenic forcing, and thus the trend cannot be excluded from the analysis. Temperature increase corresponds to decrease in TSB, thus there is a negative correlation between the two variables. Linearly detrended North Sea temperature maintains the same negative effect on TSB the following year (r = −0.48, p = 0.0025, see Supplementary Figs. S1 and S2 for detailed statistical analysis). Interestingly, the fishing mortality of 2–4-year-old cod does not exhibit a monotonous trend and does not show a strong correlation with TSB (r = −0.19, p > 0.05). This weak signal might partly be due to the fact that a decline in fishing mortality in the last years did not correspond to an increase in TSB (Fig. 1b). The low correlation exhibited by fishing mortality may limit its usage as a predictor for TSB using linear models and could indicate a time-varying F–TSB relationship typical of systems presenting discontinuous dynamics.

Fig. 1: Observed variability in temperature and cod stocks.
figure 1

a Time series of annual mean sea surface temperature from observations and the assimilation experiment (see “Methods”) for the North Sea (NOR, red lines), subpolar gyre (SPG, orange lines), and the Barents Sea Opening (BSO, purple lines). b Time series of total stock biomass (TSB, green and blue circles) and fishing mortality (F, green and blue lines) for the cod stocks in the North and Barents Sea. The inset in a shows regions over which the temperature is averaged (boxes) and the ICES ecoregions (color-filled) for delimiting cod stocks.

Full size image

The TSB of Northeast Arctic cod does not exhibit a long-term trend (Fig. 1b). This stock exhibits multi-annual to decadal variability manifested as multiple cycles of decline and increase. Similar low-frequency variability is visible in the surface temperature of the North Atlantic subpolar gyre (SPG), suggesting a possible linkage. Statistically, this linkage is supported by the high correlation between the surface temperature of the SPG and TSB of Northeast Arctic cod (r = 0.78, p = 0.0435) with the SPG-temperature leading TSB by 7 years (Supplementary Figs. S1 and S2), and consistent with previous work36. Dynamically, this linkage points to the influence of SPG circulation on the properties of Atlantic water crossing the Greenland–Scotland ridge heading towards downstream shelf seas34,35,42.

After removing respective trends from time series of the SPG temperature and Northeast Arctic cod TSB, the correlation remains high (r = 0.77, p = 0.0425), suggesting a dominating signature of decadal variability. The effect of temperature is opposite on this stock compared to the North Sea, since in the Barents Sea, temperature has a positive impact on cod biomass. These opposite impacts of temperature on biomass reflect the different temperature regimes in which the stocks reside22,43. In the case of Northeast Arctic cod, fishing mortality of 5–10-year-old cod is strongly correlated with TSB (r = −0.88, p = 0.00351). This correlation is higher than the one between temperature and TSB, and peaks at lag-2 years (Supplementary Fig. 2). For our purpose of decadal prediction of cod biomass this finding has two implications. First, the predictability horizon for TSB from a statistical point of view would be shorter with fishing mortality as a predictor compared to temperature. Second, the higher explanatory power in fishing mortality might constrain the uncertainty in the first few years of forecasts.

Statistical models for cod prediction

Once the predictors for the two cod stocks are identified, we assess various cross-validated statistical models (see “Methods”) to analyze the retrospective skill arising from the impact of temperature and fishing on the TSB and to select a model to issue forecasts. We test three different models, two simple linear regression models based on temperature and fishing mortality separately and one multiple linear regression model based on temperature and fishing mortality as explanatory variables.

As expected from the correlational analysis, the results for the North Sea cod and the North-East Arctic cod are quite different. For North Sea cod, the linear model using just fishing mortality has no predictive power (Fig. 2a and Supplementary Fig. S3 for analysis of skill from detrended variables). When the impacts of fishing and temperature are modeled together, the skill is comparable to the linear model based on temperature alone, suggesting that no additional information is gained by adding fishing mortality. For the Northeast Arctic cod, although the fishing-only model provides a better fit to the TSB data (adjusted R2 = 0.77) than the temperature-only model (adjusted R2 = 0.62), the difference in skill between these two models is not statistically significant (p = 0.15, Fig. 2b and Supplementary Fig. S3 for skill from detrended variables).

Fig. 2: Statistical models for cod prediction.
figure 2

a Cross-validated correlation skill from three linear regression models (two simple and one multiple) for TSB of North Sea cod based on North Sea surface temperature (TNOR) and fishing mortality (F). b Same as (a) but for Northeast Arctic cod and using SPG temperature (TSPG) as one of the predictors. The number in square bracket is the prediction horizon in years for each model. The dots show median skill and the whiskers show the 95% confidence limits (see “Methods”).

Full size image

Out of the three models, the model which uses both temperature and fishing has the best fit to the TSB (adjusted R2 = 0.84) and is the most suited model considering the information gained by combining SST and F (Table S1). However, both the fishing-only and the combined fishing and temperature-based models do not allow for a longer prediction horizon than the temperature-only model. This is because fishing mortality leads TSB by 2 years while temperature leads TSB by 7 years. Since our focus is on long prediction horizons, we choose the temperature-only model for the hindcast period, and for the forecast period (2020–2030), we complement the temperature-based forecasts of cod biomass with forecasts from the combined fishing and temperature-based model.

Decadal prediction of the physical environment

We now assess the prediction skill of North Sea and SPG temperature in the MPI-ESM (see “Methods” for a detailed description of this model and the decadal prediction system). In general, the skill degrades as the prediction horizon moves farther from the year of initialization (i.e. at longer lead times). However, in the North Sea, prediction skill remains high until lead year-10, and is matched by the skill from the historical simulations (Fig. 3a). This can be explained by the long-term linear trend in the underlying time series (Fig. 3b), which is present in all lead year time series. This points to the long-term trend (driven by anthropogenic external forcing) in the North Sea temperature as the source of prediction skill. Noticeable exception is the SPG where the skill is largely intact irrespective of the trend, and is higher for initialized hindcasts than historical simulations (Fig. 3c, d and Supplementary Fig. S4).

Fig. 3: Dynamical predictions of temperature.
figure 3

Anomaly correlation coefficient (ACC) as a function of lead year for initialized hindcasts (red), lagged-persistence (blue), and historical simulation (magenta dot) for annual mean (a) North Sea surface temperature and (c) subpolar gyre surface temperature with respect to the assimilation experiment for the period 1971–2019. Time series of (b) North Sea surface temperature and (d) subpolar gyre surface temperature anomalies (with respect to 1970–2019 mean) from the assimilation experiment (circles), initialized hindcast (dark colored line), forecast and historical+RCP8.5 simulation (light colored line). The solid lines in (b and d) are the respective ensemble mean predictions (or simulations) and the shading is the entire range of the respective 16-member ensemble. The regions for computing area averaged surface temperatures of the North Sea and subpolar gyre are shown in Fig. 1a. The lagged-persistence (LP)-based skill is provided for 1-, 4-, 7- and 10-year lags. The shading and whiskers in (a and c) depict 95% confidence intervals.

Full size image

The observed and predicted time series of SPG temperature suggests that during the hindcast period, most of the skill in the initialized hindcasts is derived from the ability of the model to capture the decadal cooling and warming trends (Fig. 3d). The 16-member historical simulation does not capture the full extent of the decadal variability in SPG temperature. Thus, it appears that initialization of oceanic conditions is the dominant source of predictability of SPG temperature33, while the long-term trend, mainly arising from external forcing, dominates predictability in the North Sea. The robustness of the decadal prediction skill of subpolar North Atlantic SST in the MPI-ESM-LR based decadal prediction systems has been thoroughly analyzed and is consistent with other decadal prediction systems33,44.

Dynamical–statistical cod prediction

Now, we combine the dynamical prediction of temperature with the statistical temperature–cod relationship. We choose the simplest model with temperature as the explanatory variable for both the North Sea and Northeast Arctic cod to model and forecast TSB. The utilization of temperature, derived from the dynamical model, allows us to extend the predictability horizon of cod stocks. We also include forecasts using a multiple linear regression model with fishing and temperature, and we use various scenarios of fishing mortality based on current management advice from the International Council for the Exploration of the Sea (ICES).

The dynamical–statistical prediction model shows robust skill (correlation as well as mean square error skill) in simulating the North Sea cod biomass (Fig. 4a). Note that the regression coefficients for the statistical models are not calculated from the hindcast time series of temperature, but from the observed TSB and assimilated temperatures (see “Methods”). The similarity in hindcast skill obtained from initialized hindcasts and historical simulations provides another piece of evidence that the skill is mainly due to the trend in the North Sea temperature (Fig. 4b). Our forecast of North Sea temperature for the period 2020–2030 suggest a continuation of the warm anomalies (Fig. 3c) which translates into a further decline of North Sea cod (Fig. 4a).

Fig. 4: Decadal prediction of cod stocks.
figure 4

a Time series of retrospective predictions of total stock biomass (TSB) of North Sea cod using the dynamical–statistical prediction model (retrospective predictions of North Sea surface temperature combined with the linear statistical temperature–cod relationship) for the period 1971–2019 using the initialized hindcast and historical simulation of North Sea surface temperature. The observed TSB is shown by blue circles. Also provided is the forecast for the period 2020–2030 comprising three fishing mortality scenarios: status quo (FSQ = 0.50), maximum sustainable yield (FMSY = 0.30), and precautionary approach (FLIM = 0.54). The bars and whiskers show the 95% confidence limits (2.5th, 25th, 50th, 75th, and 97.5th percentiles are shown) for the respective forecasts for the whole period (2020–2030). The historical North Sea surface temperature is extended using the RCP8.5 scenario for issuing forecasts of cod biomass. b Anomaly correlation coefficient (ACC) and mean square error skill score (MSESS) for the retrospective North Sea cod TSB prediction (1971–2019) with respect to observations. The whiskers show 95% confidence limits. c, d Same as (a, b) but for Northeast Arctic cod and using initialized hindcasts and historical simulation of SPG temperature. For the forecast, assumed fishing mortality scenarios are FSQ = 0.42, FMSY = 0.4 and FLIM = 0.74. The shadings in (a and c) show 95% confidence limits.

Full size image

In order to make our predictions of cod biomass usable in fisheries management, we provide both an SST-based forecast (for 2020–2030) and forecasts under different fishing scenarios (using the SST+F model). In particular we chose three scenarios: an FMSY scenario, in which the biomass is fished at the maximum sustainable yield (FMSY = 0.3), an FSQO scenario in which F is the mean over the last three years (FSQO = 0.5), and an FLIM precautionary scenario which is the maximum F applicable before collapse (FLIM = 0.54). The predicted total biomass of North Sea cod shows similar trends under all these scenarios, modulated in magnitude by fishing. Lower fishing initially favors a stock increase, but the constant increase of temperature leads to a further decline of the stock over time, keeping the stock in a low productivity regime. This indicates that deteriorating environmental conditions will hinder a substantial stock recovery, even with strong limitation on the fishery.

For assessing the retrospective prediction skill of Northeast Arctic cod biomass, we combine the statistical model with lead-year-4 initialized hindcasts of SPG temperature from MPI-ESM. Beyond lead-year-4, the dynamical hindcast skill degrades and is comparable to the skill from the historical simulation (Supplementary Fig. S4). Our dynamical–statistical prediction model performs well in reproducing past variability in the TSB of Northeast Arctic cod (Fig. 4c, d). Both the 1970s decline as well as the recent decadal shift in the TSB is captured by the initialized hindcast, quantified by the mean square error skill score (Fig. 4d). The correlation skill associated with historical simulation is lower but not statistically different from the hindcast skill (p = 0.164). However, the variability in the reconstructed TSB time series of Northeast Arctic cod using historical simulation is suppressed (Fig. 4c). This reconstructed time series fails to capture the recent decadal shift in the Northeast Arctic cod stock, which, as discussed above, likely follows variability in SPG temperature and is not captured by the historical simulation. This lack of variability in the reconstructed TSB time series using the historical SPG temperature is reflected in the mean square error skill score (MSESS, Fig. 4d), which suggests that this type of prediction is not significantly better than predicting a long-term mean value for the TSB.

For the Northeast Arctic cod, future predictions based on initialized hindcasts suggest a climate-driven decline of biomass in the coming decade compared to the present stock size (Fig. 4c). The RCP8.5 scenario based forecast, however, suggests a biomass level close to the long-term mean. Since the historical simulation-based hindcasts of cod biomass do not capture the full extent of past variability (MSESS for historical simulation is not significantly different than climatology), the extent of future decline in cod biomass based on RCP8.5 scenario is likely underestimated. The purely climate-driven decline in initialized forecasts is larger than in the FSQO- and FMSY-based fishing scenarios but is comparable to the decline under the FLIM scenario. Given the recent management history of this stock, the FLIM scenario is very unlikely. This could be explained by the fact that even if we are just using climate to predict cod stocks, the forecast is based on TSB levels wherein the impact of fishing is implicitly included. Cold periods in the past also coincide with periods of high F (around FLIM). This influences the forecast made using just the temperature because the statistical part of the dynamical-statistical model is trained on past TSB values. This explains why models with both fishing and temperature, where fishing is relatively low (FSQO = 0.42 and FMSY = 0.4) can maintain the stock at a higher biomass level. The forecast declining tendency (compared to the present level) in TSB of Northeast Arctic cod in all scenarios is due to the delayed (advective) impact of 2010–2016 cooling of the SPG (Fig. 1a). The future prediction of the Northeast Arctic cod is thus similar to that of North Sea cod concerning fishing mortality, indicating that a sustainable fishing pressure is necessary to maintain the stocks, but very different concerning productivity, highlighting again how climate has opposite impact on the two stocks in the next 10 years. These results provide evidence that GCM-based initialized decadal climate predictions can be deployed for prediction of marine resources through climate–ecosystem linkages.

Discussion

Sustainable management of fish stocks in the eastern North Atlantic shelf seas requires a reliable assessment of their future abundance. Incorporating environmental information in such assessment models has not always shown an improvement in prediction skills due to large uncertainties associated with recruitment–climate relationship, and also because these uncertainties might increase in a warming climate11,21. Here, we show that cod stock abundance, represented by TSB, can be successfully predicted on a decadal scale. We assess the feasibility of decadal predictions of cod stocks in the North- and Barents Sea using climate predictions from the MPI-ESM. Such an extended prediction relies on two conditions: (a) that there is a robust relationship between cod and the physical environment and (b) that the physical environment is predictable at multiyear lead times. For the North Sea, we find strong negative correlations between temperature and cod biomass, which can be explained by non-linear dynamics of the stock20,22. Ocean warming has been indicated as an important factor affecting cod in the North Sea through direct and indirect mechanisms, such as high temperatures causing low recruitment and changes in prey availability15,23,45. Fishing, on the other hand, has brought the stock close to collapse and now fishing restrictions may not be able to make the stock recover due to the detrimental effect of warming22.

We find that the long-term trend in surface temperature explains a large part of variance in the North Sea cod biomass, and consequently the high hindcast skill is largely due to the trend (externally forced). Since the detrended interannual variability in the North Sea surface temperature is not skilfully predicted by the MPI-ESM-LR (Supplementary Fig. S4), the 2020–2030 forecast for the North Sea cod biomass is mainly indicative of the long-term trajectory of the cod biomass and not of year-to-year variations around the trend. Also, future work on predictions for North Sea cod could take into account observations showing that the decline in cod abundance in the North Sea is much more pronounced in the southern North Sea than in the northern part, and there may be separate populations of cod within the North Sea management area.

The strong positive correlation between temperature and Northeast Arctic cod biomass is justified through the effect of temperature on life history traits of this stock37. While the details of how the temperature influences Northeast Arctic cod are well described16,37, the importance of the pronounced decadal variability in the SPG46, which lends predictability to the Northeast Arctic cod, is worth highlighting here. We hypothesize that the volume of Atlantic water, modulated by the SPG strength, entering the Barents Sea plays an important role. The hydrography of Norwegian–Barents Seas is related to the Atlantic inflow across the Greenland–Scotland Ridge47. When the SPG circulation is weak, the proportion of subtropical waters in the Atlantic inflow through the Faroe Shetland Channel increases35,48. The resulting increase in the volume of Atlantic water in the Barents Sea can influence the extent of sea-ice in this region, which can lead to increased productivity through extended periods of increased primary production and also due to expansion of feeding grounds. This hypothesis is consistent with the present understanding of the relationship between Atlantic heat transport and extent of sea-ice in the Barents Sea49,50 and its predictability using global coupled models51; however, the stationarity of this relationship needs to be further explored.

Interestingly, when the respective time lags between SST and cod are taken into account, the annual mean SSTs in the SPG region explain around 65% variability in the Northeast Arctic cod biomass while the local SSTs at the Barents Sea opening explain only around 12% variability. The SPG temperature is characterized by pronounced decadal variability46 while local SSTs at the Barents Sea opening prominently reflect the high-frequency atmospheric variability52 and the strong surface warming trend characteristics of these latitudes. However, the SPG signal is present in subsurface waters at the Barents Sea opening (Figs. S5 and S6). Thus local SSTs fail to capture the variability in ecosystem variables, such as the TSB, which integrate high-frequency atmospheric variability and resemble decadal temperature variability of the SPG.

A 7-year prediction horizon in Northeast Arctic cod stock has been shown to emerge from observations of SSTs in the North Atlantic but excluding the fishing mortality, and such a prediction horizon is also consistent with the length of the life cycle of Northeast arctic cod36. In the present study, we extend the predictability horizon further to a decade by using dynamically predicted SPG temperature as a predictor. Further value in our results is derived from the fact that our forecasts are based on a 16-member ensemble dynamical–statistical prediction system (see “Methods”) and various fishing mortality scenarios, which take into account the uncertainty associated with future evolution of the climate system and fishing pressure. We have also been able to identify the source of decadal prediction skill in cod stocks in the two cod habitats. In contrast to the North Sea where the externally forced trend dominates, our results emphasize decadal variability in SPG temperature as the dominant source of prediction skill in Northeast Arctic cod biomass. The predictions based on historical simulations do not capture the full extent of the decline in the cod stock in 1970s and its increase from 2005 to 2014, and hence, in terms of MSESS, these predictions do not match or outperform the predictions based on initialized hindcasts.

The approach used in this study, although novel, has certain caveats. First, the underlying climate variability that influences Northeast Arctic cod biomass has a low-frequency character. Thus, prediction skill and its uncertainty estimation is based on the assumption that the training period is representative of the climate variability associated with the subpolar North Atlantic. In case this is not so, then the skill might drop. Second, the utilization of ICES stock assessment outputs (total biomass and fishing mortality) as observations is a concern. These quantities are model outcomes, and are not entirely independent53. Third, the linear models examined here are applicable to cod stocks in our regions of interest, where the underlying oceanic variability and its impact on marine ecosystems is well understood and the stocks situated near the extremes of the species’ overall distribution range. Our models do not cover the complex issues such as those related to the impact of temperature on carrying capacity and lifetime reproductive output. This could be the subject of future work. Finally, we have assumed that the statistical models and the variables analyzed here implicitly account for possible ecosystem processes. While ecosystem processes such as species interactions are definitely important in shaping fish stocks, they are often not taken into account in management processes54, although they are to some extent taken into account in management of Barents Sea capelin (Mallotus villosus)55.

Our study attempts to bridge the gap between environmental and fisheries prediction. Through the present work, we demonstrate how decadal prediction of climate can be used to provide extended prediction horizons for fisheries combined with various fishing scenarios. Various incentives as well as the lessons learnt from past failures have motivated this effort. Foremost is the added value that such predictions can bring to the sustainable management of fish stocks. For example, at present, many fish stocks, including those considered in this article, are managed by setting annual quotas based on annual assessments of present stock size and short-term predictions (1–2 years) combined with harvest control rules based on target exploitation rates. Reliable predictions of fish biomass on a decadal scale could enable the adjustment of future catch targets (exploitation rates) to account for climate-driven fluctuations in productivity56,57. Also, predicting catch levels on a decadal scale will be important to the fishing industry, as investments in vessels, processing plants, etc. are made with a time horizon of several decades.

Climate-informed fishery management is also poised to benefit from rapid advances in multiyear prediction of other fishery-related variables such as net primary production by Earth system models58. In the North Atlantic, proper representation of open ocean-shelf connections in such models would attract further research in decadal predictions of fish stocks towards realizing a climate resilient sustainable fisheries management.

Methods

Dynamical model

The MPI-ESM is used in its low-resolution setup in the present study (MPI-ESM-LR59). The ocean general circulation component of MPI-ESM-LR, the Max Planck Institute Ocean Model (MPIOM60), is a free surface model with primitive equation solved on an Arakawa C-grid with hydrostatic and Boussinesq approximations. The MPIOM has a total of 40 z-levels in the vertical and the surface layer thickness is 12 m. The MPIOM setup used in the study has a rotated grid configuration (GR15) with one of the poles over Greenland. This enhances the horizontal resolution north of 50°N (15 km near Greenland). The resolution increases gradually to 1.5° towards the equator. Embedded in MPIOM is also the ocean biogeochemistry component, the Hamburg Ocean Carbon Cycle model (HAMOCC61). The HAMOCC incorporates oxygen and phosphate cycles, and defines the marine food web based on nutrients, phytoplankton, zooplankton, and detritus (NPZD)-based approach. The atmospheric general circulation component of MPI-ESM1.2-LR is the European Center-Hamburg (ECHAM62). The ECHAM is run at a horizontal resolution of T63 and with at total of 47 vertical levels and the model top is at 0.01 hPa. In MPI-ESM1.2-LR, the land surface–atmosphere interactions are simulated by the land vegetation module JSBACH63 which is embedded in ECHAM.

Decadal prediction system

We use one set of retrospective initialized decadal predictions (hindcasts) from the MiKlip project64, carried out with the MPI-ESM-LR. Ten-year long ensemble hindcasts with 16 members are started on 1st November every year from 1960 to 2019 (ref. 33). The initial conditions for each member come from an assimilation experiment (1960–2019) with an oceanic ensemble Kalman filter (EnKF) and atmospheric nudging. The oceanic EnKF in MPI-ESM-LR33,65 assimilates monthly profiles of temperature and salinity from EN4 (ref. 66). Simultaneously, atmospheric vorticity, divergence, temperature, and surface pressure are nudged to ERA40/ERAInterim re-analyses67. It should be noted that neither SST from satellite observations nor atmospheric temperature below 900 hPa are assimilated in order to allow for a model-consistent assimilation across the atmosphere–ocean boundary. The assimilation experiment as well as the initialized hindcasts use observed solar irradiation, volcanic eruptions, and atmospheric greenhouse gas concentrations (RCP4.5 concentrations from 2006 onward) as boundary conditions, taken from CMIP6 (ref. 68).

An additional 16-member historical simulations (1850–2005) of surface temperature taken from the MPI-ESM-LR Grand Ensemble69 are analyzed to compare the skill with the initialized hindcasts. The historical simulations are performed under natural and anthropogenic forcings derived from observations covering a total of 156 years (1850–2005). For comparison with initialized hindcasts, these historical simulations are extended with a future RCP8.5 concentrations from 2006 onward. Note that the difference between RCP8.5 and RCP4.5 scenario only emerges towards the mid of this century and hence we expect no significant impact on our short-term analysis if RCP4.5 scenario is used. The natural forcing includes solar insolation, variations of the Earth orbit, tropospheric aerosol, stratospheric aerosols from volcanic eruptions, and seasonally varying ozone. The anthropogenic forcing includes the well mixed gases CO2, CH4, N2O, CFC-11, and CFC-12 as well as O3, and anthropogenic sulfate aerosols. Atmospheric CO2 concentrations are prescribed and the carbon cycle is not interactive. It must be noted that this historical simulation is started from a pre-industrial control run and is not initialized from observations. Therefore, the internal variability in this model simulation may not be in phase with observations, and hence may not reproduce the observed timing of certain climatic events which are related to internal (natural) variability.

Linear regression models

In order to predict the time series of the TSB of cod stocks (CTSB), we construct a simple and multiple linear regression model with sea temperature (T) and fishing mortality (F) as predictors (independent variables) and the TSB as the predictand (dependent variable). For predicting North Sea cod, local oceanic surface temperature is used while for the Northeast Arctic cod, the SPG temperature is used. Both temperature time series are taken from the assimilation run as the area average of temperature of the first model layer (mid-point at 6 m depth). The time series of temperature from the assimilation run with MPI-ESM-LR compares very well with the widely used observations/re-analyses datasets, the AHOI dataset70 for the North Sea and HadISST71 for the SPG and the Barents Sea Opening (Fig. 1a). The TSB and F are taken from latest stock assessment reports from the ICES. The simple and multiple linear regression model fed with T and F anomalies (mean over 1970–2019 is removed from all variables) as predictors, for example, takes the form

$$\begin{array}{l}{C}_{{\mathrm{{TSB}}}}(y)={\beta }_{o}+{\beta }_{1}T(y-{L}_{T}),\\ {C}_{{\mathrm{{TSB}}}}(y)={\beta }_{o}+{\beta }_{1}T(y-{L}_{T})+{\beta }_{2}F(y-{L}_{F}),\end{array}$$

where CTSB is the statistical TSB prediction at year y, LT and LF are the lags in years at which the respective correlations between TSB and T or F are maximum, βo is the intercept, and β1, β2 are the slopes obtained from fitted observations.

Cross-validation of statistical models

In order to identify the best performing model, we applied the 80–20 cross-validation method. The regression coefficients are computed between time series of temperature from the assimilation run and the observed cod biomass. In the first step, the respective temperature and cod biomass time series are divided into training and testing sets by randomly selecting with replacement blocks of 80% of the parent time series as the training set and the remaining 20% as the testing set. The regression coefficients are calculated from the training set and applied to the testing set. Correlation coefficients are then calculated between the predictions made with the training set and observations as well as between the testing set and observations. This process is repeated 1000 times, and each time the 80% training set is selected randomly. The 95% confidence interval for the training and test set is the 2.5th and 97.5th percentile range of the respective 1000 correlation coefficients. Note that the lag (L) in the above equation is calculated separately for each predictor before testing various simple and multiple linear regression models based on these predictors. This procedure gives the uncertainty bounds presented in Fig. 2a, b.

Dynamical–statistical predictions

For hindcasts and forecasts, the regression model is trained on output from the assimilation run (and fishing mortality for the multiple regression model) and the resulting regression coefficients are applied to temperatures from the initialized hindcasts and historical simulation (and the fishing mortality scenarios for multiple regression models). The statistical model is fed with anomalies of each variable and the mean is added to the predicted TSB anomalies at the end. Mathematically this takes the form

$$\begin{array}{l}{C}_{\mathrm{{TSB}}}^{\prime}(y)={\beta }_{o}+{\beta }_{1}T^{\prime} (y-{L}_{T}),\\ {C}_{\mathrm{{TSB}}}^{\prime}(y)={\beta }_{o}+{\beta }_{1}T^{\prime} (y-{L}_{T})+{\beta }_{2}F(y-{L}_{F}),\end{array}$$

where \({C}_{\mathrm{{TSB}}}^{\prime}\) is the dynamical-statistical TSB prediction at year y, \({{T}}^{\prime}\) is the dynamically predicted temperature (lead-year-10 predictions for the North Sea and lead-year-4 for the SPG), LT and LF are the lags in years at which the respective correlations between observed TSB and T or F are maximum, βo is the intercept, and β1 and β2 are the slopes obtained from fitted observations.

The uncertainties in regression coefficients (slopes and intercepts) are also estimated using a bootstrapping methodology. First, 1000 new predictor and predictand time series of same length as the original time series are constructed by random sampling with replacement from the parent time series, while preserving their relationship. These new time series are then used to get 1000 estimates of regression coefficients. These 1000 regression coefficients are then applied to each of the 16 ensemble members (for temperature as the predictor). The 95% confidence interval is the 2.5th and 97.5th percentile range of these 16,000 predictions. This procedure gives the uncertainty bounds presented in Fig. 4a, c

Hindcast skill and hindcast uncertainty

We use anomaly correlation coefficient (ACC) and the MSESS as measures of skill of initialized hindcasts and historical simulations against observations (stock assessment for TSB and assimilation output for temperature) for the period 1960–2019. The MSESS is defined as

$${\mathrm{{MSESS}}}=1-{\mathrm{{MSE}}}/{\mathrm{{MSE}}}_{\mathrm{{REF}}}$$

where MSE is the mean square error of prediction and MSEREF is the mean square error of reference forecast (here, climatology is used as reference)

Prior to calculating ACC and MSESS (and also prior to feeding the statistical model for TSB), the initialized hindcasts are corrected for the lead-time-dependent drift72, and lead-year-dependent climatology (mean over 1970–2019) is removed. The uncertainty in hindcast skill is determined using a block bootstrapping approach. The bootstrapping is done both in time and across ensemble members. We use a 6-year overlapping block bootstrap to account for the autocorrelation in the time series. The estimated uncertainties are not sensitive to a reasonable choices of block length that allow sufficient number of blocks for sampling. Through random resampling with replacement, 1000 new block-bootstrapped time series of predictions and observation are used to obtain 1000 new estimates of ACCs. The 95% confidence interval is the 2.5th and 97.5th percentile range of these 1000 ACCs or MSESSs. This procedure gives the uncertainty bounds presented in Figs. 3a, c and 4b, d.