Ocean forcings
The Nucleus for European Modelling of the Ocean (NEMO) ocean framework46, which includes an online coupling with the biogeochemical component PISCES in a 2° latitude × 2° longitude configuration47,48, was used to simulate the historical oceanic environment (hindcast simulation). This historical simulation was forced by the Drakkar Forcing Sets 5.2 (DFS5.2)49 on the basis of a corrected set of the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis – Interim (ERA-Interim) over the period 1979−2011. Salinity, temperature and biogeochemical tracer concentrations (nitrate, phosphate, iron, silicate, alkalinity, dissolved oxygen and dissolved organic and inorganic carbon) were initialized from the World Ocean Atlas climatology (WOA09)50 and previous model climatology for iron and dissolved organic carbon51. To minimize any substantial numerical drift in the simulations related to a non-equilibrated initial state, we applied a spin-up of the ocean model and biogeochemical model for 66 years, cycling twice over the DFS5.2 forcing sets48.
Overall, the model simulates basin-scale, historical SST and salinity distribution, together with seasonal and interannual (ENSO) variability with good fidelity52. Classical biases are associated with the coarse (2°) resolution, for example, the latitudinal position of the Kuroshio Current. In the tropical Pacific, there is a cold bias of −1 °C in the central equatorial zone (between 170° W and 100° W) and a warm bias of +1 °C in the eastern part of the basin (east of 90° W). Despite some local discrepancy between simulation outputs and satellite-derived chlorophyll concentration around islands and near the American coasts, simulated mean chlorophyll in the equatorial Pacific Ocean is close to observed values51,52.
For future ocean projections, we first selected several ESMs from the CMIP5 intercomparison project53 on the basis of the ability of the models to produce accurate ENSO variability in the Pacific54. The four ESMs selected were IPSL-CM5A55, MIROC56, GFDL-ESM2G57 and MPI-MR58. We then extracted atmospheric fields from these models for the period 2011−2100 under RCP 8.5 to simulate ‘business-as-usual’ climate anomalies to build forcing sets for the NEMO–PISCES ocean model.
All ESMs display large biases in their representation of Pacific climate, including the important South Pacific Convergence Zone59,60. These atmospheric biases propagated uncertainties associated with future atmospheres into the coupled, dynamical-biogeochemical oceanic framework. For example, they result in prominent distortions in the extension and position of the warm pool61 and can be expected to affect modelling of the open ocean ecosystem up to the higher trophic levels12.
To mitigate the mean state model biases in the selected ESMs, we used a ‘pseudo-warming’ anomaly approach to force the ocean model. To do this, we extracted monthly anomalies (relative to 2010) of surface atmospheric temperature, zonal and meridional wind speeds, radiative heat fluxes, relative humidity and precipitation from the ESM models over the 2010–2100 period and applied a 31-year-wide Hanning filter to remove variability on timescales less than 15 years.
Each ESM-filtered timeseries was superimposed onto the repeating 30-year historical forcing (that is, repeated three times to span the twenty-first century) to provide the forcing for the NEMO–PISCES projections. This procedure enabled us to retain a realistic climatology and high-frequency variability from observations subject to long-term trends due to climate change based on the ESMs (Supplementary Fig. 7).
For consistency, the control simulation of NEMO–PISCES was forced using the same three, repeated, 30-year historical periods to correct any long-term drift generated internally without climate change forcing.
It is important to note that use of all ESM acronyms (for example, IPSL) in the following text refers to NEMO–PISCES or SEAPODYM simulations derived from the ESM anomaly forcing, and not to the ESM models themselves.
The four NEMO–PISCES simulations of future ocean conditions produced contrasting results in terms of dynamics and biogeochemistry (Supplementary Fig. 8). In particular, there was strong warming in the IPSL and MIROC simulations and weaker warming for GFDL and especially MPI. Spatial patterns in ocean warming produced by the NEMO–PISCES simulations differed mostly in intensity rather than spatial structure.
Using NEMO–PISCES outputs to produce SEAPODYM forcings
The outputs of NEMO–PISCES were used to provide environmental forcing variables for SEAPODYM, the model used to project the responses of the key life stages of skipjack, yellowfin and bigeye tuna to climate change (Supplementary Note 7). The following physical and biochemical forcing variables were used in SEAPODYM applications: three-dimensional (3D) temperature, dissolved oxygen (O2) concentration, zonal/meridional currents and primary production, and 2D euphotic depth. Before running SEAPODYM, these forcing variables were interpolated to a regular 2° Arakawa A grid and placed in the centre of the grid cells. Primary production was then vertically integrated throughout the water column, whereas the other 3D variables were integrated within three pelagic layers, defined according to the euphotic depth to provide the mean 2D fields for each variable per layer. Selected environmental variables from the historical ocean reanalysis and from four climate-driven ocean outputs are shown in Supplementary Fig. 3.
These integrated variables were then used to force the SEAPODYM-LMTL (lower and mid-trophic level) sub-model. SEAPODYM-LMTL relies on primary production, temperature and ocean currents to simulate the biomass of six functional groups of micronekton—mid-trophic-level prey organisms of tunas (Supplementary Fig. 4)—residing or migrating through three pelagic layers within the upper 1,000 m of the water column (the epipelagic layer and the upper and lower mesopelagic layers), with depths linked to the depth of euphotic layer Z as 1.5Z, 4.5Z and 10Z (with 10Z limited to 1,000 m). The definition of these pelagic layers is derived from the diurnal vertical distributions of micronekton species62.
Optimal parameterization of SEAPODYM during historical period
The parameterization of SEAPODYM for each tuna species is highly sensitive to ocean forcing; that is, in its average state it is free from systematic biases, and it represents interannual variability and ENSO correctly. This sensitivity enables the model to reproduce observed variability within large, geo-referenced datasets of tuna catches and length distributions reflecting changes in fish abundance12. The environmental forcings in this study were obtained from the historical NEMO–PISCES reference simulations using a realistic atmospheric reanalysis based on a consistent set of atmospheric observations. Historical fishing datasets used to achieve model optimal parameterizations were compiled from the combination of data provided by the Pacific Community for the WCPO and by IATTC for the EPO. The model spatial resolution was 2° × 2°, and the resolution for time and age dimensions was one month. The skipjack tuna reference model was obtained by integrating all available geo-referenced data—catch, length-frequency of catch and tagging release–recapture data—into a likelihood function and obtaining the solution using the maximum likelihood estimation (MLE) approach (Supplementary Note 7). The initial habitat and movement parameters for bigeye and yellowfin tuna were also estimated by integrating tagging data into the model; however, the final parameterizations of the reference models for these two species were based mainly on fisheries data. The methodology and optimal reference solutions obtained for skipjack, yellowfin and bigeye tuna, and model validations with statistical metrics, are described in other publications documenting the use of SEAPODYM13,63,64,65.
The structures of the populations of the three tuna species in December 2010 (the last time-step of the reanalysis) were used to set the initial conditions for the projections starting in 2011. A second historical simulation was run to remove the effects of fishing mortality (Supplementary Figs. 9 and 10) to establish the initial conditions for the unfished tuna populations (Supplementary Fig. 10). In these latter simulations, the stocks increase and reach an equilibrium state in a time that is defined by the lifespan of the species and the estimated stock–recruitment relationship. We assume that at the end of the 30-year reanalysis (December 2010), stocks of all three tropical tuna species are at their virgin (unfished) state and influenced by environmental variability and demographic processes only.
Projections of climate change impacts on tuna
Previous studies on the impact of climate change on tropical tuna species in the Pacific Ocean produced projections based on the full-field NEMO–PISCES output from a single ESM (IPSL) under the IPCC business-as-usual scenario6,10,12,66,67. These projections were subject to biases, resulting in poor coherence between historical and projected environmental forcings and abrupt changes and biases when switching from a historical reanalysis to a projected time series12. To reduce this problem, we used an approach based on the four, bias-corrected, projected climates from NEMO–PISCES outputs (Supplementary Methods).
Simulations of the SEAPODYM tuna model were run with parameters from the reference MLE models for the three tuna species, with forcings from the four NEMO–PISCES and mid-trophic simulations, under the RCP 8.5 scenario to project tuna population dynamics until mid-century. We estimated the virgin biomass of each species in the decade 2011−2020 and computed the relative change in biomass by 2050 (2044−2053) as follows:
$$updelta _Bleft( {2050} right) = frac{1}{N}mathop {sum}limits_{t = 2011}^{2020} {left( {frac{{Bleft( {t + {Delta}t} right)}}{{B(t)}} – 1} right)}$$
(1)
where Δt is the time interval corresponding to 33 years and N is the number of monthly time steps in the selected time period (120 months between 2011 and 2020). We chose to average over 10 years at 33-year intervals to compare two distant periods with the same atmospheric variability, thus removing the possible effects of interannual variation and allowing better detection of the climate change signal.
The relative biomass change δB (2050) was computed for the EEZs of Pacific SIDS and all high-seas areas in the WCPO and EPO (Supplementary Fig. 1).
Sensitivity analyses to explore uncertainty
We analysed the impacts of climate change on skipjack, yellowfin and bigeye tuna with an ensemble of simulations focusing on the greatest sources of uncertainty in the NEMO–PISCES variables and in SEAPODYM (Supplementary Fig. 11 and Supplementary Table 21). The methods used to explore these uncertainties, and the rationale for these analyses, are explained in the Supplementary Methods.
Modelling tuna distribution under lower-emissions scenarios
The simulations based on RCP 8.5 project a redistribution of tuna biomass by 2050 as globally averaged surface temperature rises to 2 °C above pre-industrial levels by mid-century. To evaluate possible effects of a lower GHG emission scenario on tuna redistribution, we also estimated the responses of tropical tuna species to conditions similar to RCP 4.5 and RCP 2.6 by 2050.
In the absence of ocean forcings and SEAPODYM outputs for RCP 4.5 and RCP 2.6, we used estimates based on the RCP 8.5 simulations using a ‘time-shift’ approach68. This method consists of identifying the time segment in RCP 8.5 in which a key variable (for example, CO2-equivalent (CO2e)) matches the value expected for the selected RCP in 2050. Accordingly, we selected the periods in the RCP 8.5 curve when total CO2e concentrations in the atmosphere reached those projected for RCP 4.5 and RCP 2.6 in 2050 (Supplementary Fig. 12). On the basis of this method, the equivalent of RCP 4.5 in 2050 is reached in 2037 under RCP 8.5, and the equivalent for RCP 2.6 in 2050 is reached in 2026.
An important assumption of this method is that the dynamical pattern corresponding to a given change of global temperature is independent of the rate of change. This assumption is expected to be met for key features of the tropical Pacific Ocean because the upper ocean generally responds rapidly to changes in atmospheric forcing. However, this assumption is unlikely to hold for tuna population dynamics because interannual variability of tuna biomass is driven by demographic processes (recruitment and mortality), which are in turn influenced by environmental variability. Furthermore, due to the slow nature of demographic processes, the repercussions of environmental variability on tuna population dynamics are time lagged. For example, there is a time lag of 8 months between the Southern Oscillation Index and the biomass of young skipjack tuna (aged from 3 to 9 months)17, and a time lag of 12 months between the Southern Oscillation Index and total biomass of skipjack tuna (Supplementary Fig. 13). When combined with the effects of stock–recruitment relationships, and different generation times between tuna species, the speed and duration of climate change processes may have a profound effect on tuna biomass. Therefore, due to the rapidly changing ocean conditions in the RCP 8.5 scenario, the population status of a tuna species in the second and third decade cannot be assumed to be equivalent to that under a scenario with lower emissions by mid-century.
To address the complications associated with the population dynamics of tuna in a changing environment, we generated synthetic RCP 4.5 and RCP 2.6 2011−2050 time series by recycling the years from RCP 8.5 simulations. Note that recycling the ‘equivalent’ years from RCP 8.5 simulations to imitate those projected for the RCP 4.5 and RCP 2.6 scenarios involves re-using the same years multiple times because of their lower rate of change. To avoid looping the forcings over the same year multiple times, we selected several years around the equivalent RCP 8.5 year while enlarging the temporal window with increasing differences in the rates of GHG change between the two scenarios and ensuring that the mean CO2e within this window was equal to those in the target RCP 4.5 or RCP 2.6 scenario. The inverse mapping of the RCP 8.5 curve from arrays of CO2e values to the equivalent years in the RCP 8.5 simulation (Supplementary Fig. 14) provided the selected range of RCP 8.5 years to imitate the RCP 4.5 and RCP 2.6 scenarios. The NEMO–PISCES model variables from those years were then used to compute monthly climatology for each year of the surrogate RCP 4.5 or RCP 2.6 forcing to provide smoothed time series of forcing variables over the complete time range. The temporal evolution of epipelagic ocean temperature is compared for four climate models and three RCP scenarios in Supplementary Fig. 14.
The biomass changes projected for the three tuna species in 2050 under RCP 8.5 and under the lower surrogate emissions scenarios were then computed for all Pacific Island EEZs (Supplementary Fig. 15) following equation (1) (Supplementary Methods). The biomass changes projected under the RCP 4.5 forcing are smaller in magnitude than those for RCP 8.5, demonstrating that the effect of climate change is less pronounced in the simulations under this lower-emissions scenario.
The simulations under the surrogate RCP 2.6 forcing did not follow the expected pattern and were deemed to be too unreliable for use in this study (Supplementary Methods).
Estimating changes in tuna biomass in EEZs and the high seas
For this analysis, we produced reference biomasses for skipjack, yellowfin and bigeye tuna for the period 1979−2010 from quantitative assessment studies using SEAPODYM, which estimates population dynamics, habitats, movements and fisheries parameters with an MLE approach (Supplementary Note 7). The fit between observations and predictions (for catch and catch size frequencies) was used to validate the optimal solutions of the models within and outside the time window for the model parameter estimates. The fit was analysed spatially by fishery to ensure that there were no regional biases. Once the optimal solution was achieved, a final simulation was made with the same set of parameter estimates but without considering any fishing, to obtain the unfished biomass dynamics during both the historical period and the projection for the twenty-first century. The differences in unfished biomass between the historical period (2001−2010) and projections in 2050 (mean of 2046−2050) for each species were used to compute the weighted mean change in total tuna biomass in the EEZs of the ten Pacific SIDS, the high-seas areas shown in Supplementary Fig. 1 and the EEZs of the other Pacific SIDS listed in Supplementary Table 1 for the RCP 8.5 and RCP 4.5 emission scenarios by 2050.
Estimating changes in catch in EEZs and the high seas
To evaluate the impacts of climate change scenarios on purse-seine fisheries, comparisons were restricted to the EEZs of the ten tuna-dependent Pacific SIDS and the high-seas areas, particularly EPO-C (Supplementary Fig. 1).
To estimate the effects of projected changes in biomass of skipjack, yellowfin and bigeye tuna due to RCP 8.5 and RCP 4.5 on purse-seine catches in the EEZs of Pacific SIDS and in high-seas areas by 2050, in the absence of management interventions to reallocate catch entitlements to maintain historical access rights for Pacific SIDS, we assumed that there would be a direct relationship between projected changes in biomass and catch. Because purse-seine catches are composed of different proportions of the three tuna species, and because each species is projected to have a different response to climate change (Fig. 2), changes in purse-seine catches by 2050 were estimated using the weighted mean response of the three tuna species to RCP 8.5 and to RCP 4.5. These estimates were derived from the average relative abundance of each species in purse-seine catches in the EEZs of the ten Pacific SIDS (Supplementary Table 3) and in high-seas areas (Supplementary Table 4) and the projected percentage change in biomass of each species under each emission scenario (Supplementary Tables 17 and 18).
The weighted average percentage changes in biomass of all tuna species combined were then applied to the 10-year average (2009−2018) purse-seine catches from the EEZs of the ten Pacific SIDS and high-seas areas (Supplementary Tables 3 and 4) to estimate the changes in purse-seine catches for these jurisdictions by 2050 under RCP 8.5 and RCP 4.5. In the case of Kiribati, which has three separate EEZ areas (Fig. 1), we estimated the change in catch for each EEZ area and amalgamated the results to produce the overall estimated change in purse-seine catch for the country.
The projected percentage change in total purse-seine catch differs from the percentage change in total tuna biomass due to variation in the relative contributions of the three tuna species to total catch and to total biomass.
Estimating the effects of tuna redistribution on economies
To assess the effects of climate-driven redistribution of tuna on the economies of the 10 Pacific SIDS, we assumed that estimated changes in purse-seine catch within their EEZs due to the redistribution of tuna biomass described above would result in a proportional change in access fees earned from purse-seine fishing and associated operations.
To estimate the effects of RCP 8.5 and RCP 4.5 on the capacity of Pacific Island governments to earn access fees from industrial tuna fishing, and the contributions of these access fees to total government revenue excluding grants (‘government revenue’), we used annual averages of government revenue, tuna-fishing access fees earned by the ten Pacific SIDS and the percentage contribution of access fees to government revenue for the period 2015−2018 (Supplementary Table 2) as a baseline. We applied the projected average percentage changes in total purse-seine catch in each EEZ for RCP 8.5 and RCP 4.5 (summarized in Supplementary Tables 17 and 18) to the average annual access fees received in 2015−2018 by each of the Pacific SIDS to estimate the change in value of their access fees by 2050 under each emissions scenario. The change in value of access fees was used to estimate decreases or increases in government revenue in 2050 relative to 2015–2018 under both emissions scenarios in US$ and percentage terms, assuming that the relative contributions of other sources of government revenue remain the same.
The estimated percentage changes in government revenue for each Pacific SIDS do not account for (1) management responses; (2) variation in the value of access to particular EEZs and the willingness of fleets to pay for this access due to the effects of changes in tuna biomass on catchability of each species, levels of fishing effort/catch rates, the price of tuna or cost of landing tuna; and (3) the impact of tuna redistribution on the degree of control that Pacific SIDS exert over fisheries targeting tuna. The third factor is expected to be particularly important. For example, substantial movement of tuna from the EEZs of PNA countries into high-seas areas would be expected to limit the effectiveness of the VDS69 by reducing the degree of control over the fishery exerted by PNA members.
Overall, it is important to note that the simple approach used to assess the potential effects of tuna redistribution on government revenue is intended only to provide indicative information on the magnitude of these impacts. To obtain robust estimates of climate-driven changes in government revenue, more complex bio-economic analyses will be required, beginning with, for example, a fleet-dynamics analysis to investigate the potential response of purse-seine vessels to redistribution of tuna and the flow-on effects on access fees.
Reporting Summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Source: Ecology - nature.com
