DOI: 10.7291/d1jq2k
title: Recommendations for quantifying and reducing uncertainty in climate projections of species distributions

description: Projecting the future distributions of commercially and ecologically
 important species has become a critical approach for ecosystem managers to
 strategically anticipate change, but large uncertainties in projections
 limit climate adaptation planning. Although distribution projections are
 primarily used to understand the scope of potential change - rather than
 accurately predict specific outcomes - it is nonetheless essential to
 understand where and why projections can give implausible results and to
 identify which processes contribute to uncertainty. Here, we use a series
 of simulated species distributions, an ensemble of 252 species
 distribution models, and an ensemble of three regional ocean climate
 projections, to isolate the influences of uncertainty from earth system
 model spread and from ecological modeling. The simulations encompass
 marine species with different functional traits and ecological preferences
 to more broadly address resource manager and fishery stakeholder needs,
 and provide a simulated true-state with which to evaluate projections. We
 present our results relative to the degree of environmental extrapolation
 from historical conditions, which helps facilitate interpretation by
 ecological modelers working in diverse systems. We found uncertainty
 associated with species distribution models can exceed uncertainty
 generated from diverging earth system models (up to 70% of total
 uncertainty by 2100), and that this result was consistent across species
 traits. Species distribution model uncertainty increased through time and
 was primarily related to the degree to which models extrapolated into
 novel environmental conditions but moderated by how well models captured
 the underlying dynamics driving species distributions. The predictive
 power of simulated species distribution models remained relatively high in
 the first 30 years of projections, in alignment with the time period in
 which stakeholders make strategic decisions based on climate information.
 By understanding sources of uncertainty, and how they change at different
 forecast horizons, we provide recommendations for projecting species
 distribution models under global climate change.
description: <em>Summary</em> We used a
 combination of regional ocean climate projections and simulated species
 distributions (Leroy et al., 2016) to quantify sources of uncertainty in
 projections of spatially-explicit biomass for three species archetypes in
 the CCS (1985-2100; Fig. 1). Species archetypes were simplified
 representations of three general groups of marine finfish found in the CCS
 that comprise ecologically and/or economically important fisheries and
 that might be expected to show variable patterns of redistribution under
 climate change based on their habitat preferences, population dynamics,
 and mobility characteristics: 1) a highly migratory species (HMS) that was
 designed to resemble north Pacific albacore; 2) a coastal pelagic species
 (CPS) that was designed to resemble northern anchovy (CPS); and 3) a
 groundfish species (GFS) that was designed to resemble sablefish. SDMs
 (n=15; Figure 1) were then fitted to simulated biomass data for each
 archetype (training period 1985-2010) and projected from 2011-2100 using
 each of the three regional ocean climate models. Our framework resulted in
 252 SDMs (15 SDM types, three species archetypes, three ESMs, and two
 environmental parameter simulations; Figure 1). To address our study goal
 of assessing SDM performance and understanding sources of uncertainty in
 species distribution projections, we compared the output of SDM
 projections against simulated “observations” for 2011-2100 and quantified
 the uncertainty introduced by the climate projection (ESM uncertainty)
 versus the uncertainty introduced by the SDM structure (SDM
 uncertainty). <em>Environmental Covariates from
 Regional Ocean Projections</em> Environmental
 covariates used in species distribution simulations were obtained from
 regional ocean projections (Pozo Buil et al., 2021) forced by three ESMs
 from phase 5 of the Coupled Model Intercomparison Project (CMIP5) archive:
 Geophysical Fluid Dynamics Laboratory (GFDL) ESM2M, Hadley Center
 HadGEM2-ES (HAD), and Institut Pierre Simon Laplace (IPSL) CM5A-MR. These
 ESMs, hereafter referred to as GFDL, HAD, and IPSL, span the approximate
 range of potential changes in physical and biogeochemical conditions
 across all CMIP5 models (Pozo Buil et al., 2021). ESMs were downscaled
 using the Regional Ocean Modelling System (ROMS) coupled with a
 biogeochemical model (NEMUCSC) (Fiechter et al., 2018, 2021) based on the
 North Pacific Ecosystem Model for Understanding Regional Oceanography
 (NEMURO) (Kishi et al., 2007). The ROMS domain spans the CCS from 30-48°N
 and from the coast to 134°W at 0.1° horizontal resolution with 42
 terrain-following vertical layers (Figure 2). Each downscaled ESM used the
 Representative Concentration Pathway (RCP) 8.5 climate change scenario.
 While we only examined RCP 8.5, it should be noted that using RCPs 2.6 and
 4.5 would result in only minor differences in the spread of future
 environmental change for the variables and ESMs examined here.
 Specifically, uncertainty in biogeochemical change among the chosen ESMs
 in RCP8.5 envelops the uncertainty among RCPs 2.6 and 4.5; while for
 temperature GFDL and HAD represent opposite ends of the spectrum for the
 projected magnitude of warming in the CMIP5 ensemble (Drenkard et al.,
 2021; Pozo Buil et al., 2021). As such, we do not explore scenario
 uncertainty. Environmental covariates used in species distribution
 simulations were sea surface temperature (SST; C), bottom temperature (BT;
 C), bottom oxygen (BO; mmol m-3), mixed layer depth (MLD; m), surface
 chlorophyll-a (Chl-a; mg m-3), and zooplankton concentration integrated
 over 50 m (zoo_50; mmol N m-2) and 200 m (zoo_200; mmol N m-2). These
 environmental covariates were averaged over spring months (March-May)
 annually (1985-2100) to encompass the seasonal period when ocean
 productivity is most influential on the long-term population dynamics of
 most marine fishes in the CCS. <em>Operating
 Models: Simulated Species Biomass</em> Biomass
 distributions for three species archetypes were simulated on the ROMS grid
 for each year and each ESM from 1985-2100. Simulations were run using the
 ‘virtualspecies’ R package (Leroy et al., 2016) that is specifically
 designed to reflect real-world ecological properties and
 species-environment relationships (Meynard et al., 2019). We refer to
 these simulated species distributions as ‘operating models’.  Species
 simulations used a two-step process. First, habitat suitability was
 calculated based on environmental data and specified species’ habitat
 preferences (Table S1). Environmental preferences used to force species
 distributions varied among species archetypes based on representative life
 histories (see Supplementary Material). The domain for the HMS archetype
 was set to the entire CCS, whereas the CPS and GFS archetypes were reduced
 to inshore waters to reflect the CPS archetype’s preference for pelagic
 waters over the continental shelf and slope, and the GFS archetype’s
 preference for demersal shelf and slope habitats (Leeuwis et al., 2019;
 Stierhoff et al., 2020). Second, total habitat
 suitability was calculated, and converted to presence-absence using a
 logistic function (which specifies at what suitability value the species
 becomes present). When species were present, biomass was estimated from a
 log-normal distribution, and when species were absent biomass was set to
 zero. Biomass at each grid cell was multiplied by habitat suitability of
 that same grid cell to provide habitat-informed biomass. For CPS and GFS
 archetypes, an additional biomass multiplier was used to encompass
 population-level dynamics (Figure S1; see supplementary methods) (Punt et
 al., 2016). Specifically, CPS biomass was made to reflect boom-bust
 population dynamics that are common in CPS species in the CCS, while GFS
 biomass integrated a 20-year phase shift between low and high recruitment,
 as has been observed for sablefish (Haltuch et al., 2019). Simulated data
 were generated for each grid cell (HMS = 21912 grid cells; CPS &
 GFC = 4012 grid cells) once per year for 116 years (1985-2100). Detailed
 methods for the simulation are provided in the Supplementary material, and
 R code is provided on GitHub (https://github.com/stephbrodie1/Projecting_SDMs). <em>Estimation Models: Species Distribution Models</em> We parameterized a series of SDMs to estimate the relationship between simulated species biomass and covariates (Figure 1). Because these are fitted to data from an operating model, we refer to these SDMs as ‘estimation models’. Multiple approaches were tested to explore how decisions about model type and parameterization influence model accuracy and predictive performance (Brodie et al., 2020). We used four types of SDMs: generalized additive models (GAM), generalized linear mixed models (GLMM), boosted regression trees (BRT), and multilayer perceptron models (MLP; a type of artificial neural network model) (Table S2). Parameterization options included various combinations of environmental (E), spatial (S), and temporal (T) covariates (Figure 1; see supplementary methods). Spatial and temporal covariates can act as proxies for unobserved or unmeasured processes that drive species distributions, and were included here given their common use in SDMs (typically called spatiotemporal models) (Brodie et al., 2020). We expect spatiotemporal SDMs with no environmental covariates to perform poorly over the projection period. We constructed all SDMs as delta (hurdle) models, where the probability of occurrence (binomial) and positive biomass (log-normal) were estimated as separate processes. All SDMs were trained on data from 1985-2010, where only 500 random samples per year (2% of available data) were used for fitting (n=13 000). Random samples included both presence and absences sampled across the entire domain. No SDM validation or model selection was required as our simulation experiment is designed to explore a range of model parameterizations.   Fitted SDMs were then used to predict species biomass on projected environmental data, for every year and grid cell in the domain. Only 500 randomly sampled grid cells per year (2011-2100) were used for testing purposes (n=45 000), to match the resolution of samples used to train models. Importantly, not all environmental covariates used to simulate species biomass (see 2.3 above) were included in the fitted SDMs. Specifically, we used chlorophyll-a as a proxy for prey fields (zooplankton) to approximate real-world conditions where imperfect information is available for estimating species’ habitat preferences. In addition to the 15 SDM parameterizations listed in Figure 1, we examined SDMs that only contained a single covariate of temperature (either surface or bottom temperature depending on the archetype). This experiment was done to test how under-parameterized models that miss key environmental drivers of species distributions performs, and the degree to which this approach decreases model fit and increases projection uncertainty. We refer to these SDMs as ‘temperature-only’ models (Figure 1). 
description: .rds and raster files can be opened in R statistical
 software.