Future Physical and Biogeochemical Ocean Conditions under Climate Change along the British Columbia Continental Margin

ABSTRACT

Climate change impacts coastal ecosystems through large scale changes in temperature, stratification, circulation and ocean acidification. Here, the potential response of the British Columbia continental margin to climate change is investigated using a regional ocean circulation-biogeochemical model to downscale climate change projections from the Canadian regional and global climate models (CanRCM4/CanESM2) under two Intergovernmental Panel on Climate Change emission scenarios (RCP 4.5 and RCP 8.5). Projections of future physical and biogeochemical conditions for the 2041–2070 period are compared to the recent past (1981–2010). We found an overall annual average warming of >1.6°C in sea surface temperature, increase in stratification in the upper layer, and decrease in surface pH of as much as 0.21. Increasing stratification and changing winds have a limited impact on nitrate availability, phytoplankton biomass and primary production, whilst ocean warming increases primary production by up to 30% in most of the model domain. Increased atmospheric CO₂ contributes to acidification over the model domain with a decrease in pH and aragonite saturation (Ω_arag) at all depths resulting in an increase of 20 to 32% of the volume of Ω_arag ≤1 in the upper 100 m of the continental shelf depending on the climate scenario. Our projected results, therefore, show that future climate change may alter the amount of food available for higher trophic levels and the habitat of benthic species, since bottom waters on the shelf will be undersaturated with respect to aragonite for 2–3 months in mid-summer. Both climate change scenarios results in a similar pattern of changes but projected changes were stronger and more extensive under RCP 8.5 showing the benefit of mitigation efforts in reducing the effect of climate change on marine ecosystem stressors.

RÉSUMÉ

[Traduit par la rédaction] Le changement climatique a des répercussions sur les écosystèmes côtiers en raison des changements à grande échelle de la température, de la stratification, de la circulation et de l’acidification des océans. Ici, la réponse potentielle de la marge continentale de la Colombie-Britannique au changement climatique est étudiée au moyen d’un modèle régional de circulation océanique et biogéochimique pour réduire les projections de changement climatique des modèles climatiques régionaux et globaux canadiens (CanRCM4/CanESM2) dans le cadre de deux scénarios d’émissions [BDE1] (RCP 4.5 et RCP 8.5) du Groupe d’experts intergouvernemental sur l’évolution du climat (GIEC). Les projections des conditions physiques et biogéochimiques futures pour la période 2041–2070 sont comparées au passé récent (1981–2010). Nous avons observé un réchauffement annuel moyen global de >1,6 °C de la température de surface de la mer, une augmentation de la stratification dans la couche supérieure et une diminution du pH de surface pouvant aller jusqu’à 0,21. La hausse de la stratification et le changement des vents ont une incidence limitée sur la disponibilité des nitrates, la biomasse du phytoplancton et la production primaire, tandis que le réchauffement de l’océan augmente la production primaire jusqu’à 30% dans la majeure partie du domaine du modèle. L’augmentation du CO₂ atmosphérique contribue à l’acidification sur le domaine du modèle avec une diminution du pH et de la saturation en aragonite (W_arag) à toutes les profondeurs entraînant une augmentation de 20 à 32% du volume d’W_arag ≤1 dans les 100 m supérieurs du plateau continental en fonction du scénario climatique. Nos résultats projetés montrent donc que le changement climatique futur pourrait modifier la quantité de nourriture disponible pour les niveaux trophiques supérieurs et l’habitat des espèces benthiques, puisque les eaux de fond du plateau seront sous-saturées en aragonite pendant 2 à 3 mois au milieu de l’été. Les deux scénarios de changement climatique donnent lieu à des changements similaires, mais les changements projetés sont plus importants et plus étendus dans le cadre du RCP 8.5, ce qui révèle l’intérêt des efforts d’atténuation pour réduire l’effet du changement climatique sur les facteurs de stress de l’écosystème marin.

KEYWORDS:

• climate change
• British Columbia coast
• regional climate model
• biogeochemistry
• downscaled projections

1 Introduction

The British Columbia (BC) continental margin along the Canadian Pacific coast is located within a dynamic transition zone, where the Pacific Ocean Current bifurcates into the northeastward-flowing Alaska Current and the southeastward flowing California Current (Cummins & Freeland, 2007). This region is very productive and serves as habitat of many economically important species, including salmon, herring, and Pacific hake (e.g. Ware & Thomson, 2005). Despite its importance, the dynamics of this ecosystem and their susceptibility to climate change and ocean acidification remains poorly understood. Primary productivity is strongly influenced by summer wind-driven coastal upwelling off the Vancouver Island shelf, which injects nutrients into surface waters (Peña et al., 2019; Whitney et al., 2005) but also delivers carbon-rich (Feely et al., 2008) and oxygen-poor waters (Crawford & Peña, 2013) into the shelf making this region particularly vulnerable to future ocean acidification and hypoxia (Feely et al., 2008; Hauri et al., 2013). In addition, the estuarine outflow from the Juan de Fuca Strait feeds the near-shore Vancouver Island Coastal Current preventing denser hypoxic water from reaching the shore and supplies ample amounts of nutrients to the Vancouver Island coast year-round (Bianucci et al., 2011; Crawford, 1991; Crawford & Peña, 2013). In northern BC (Fig. 1), the general circulation is estuarine in summer and coastal winds force strong downwelling in winter and weak upwelling in summer as downwelling relaxes (Hannah et al., 1991). The broad continental shelf of this region is also highly productive as nutrient-rich waters flow onto the shelf at depth in summer and are subsequently entrained into surface waters (e.g. Whitney et al., 2005). All the aforementioned physical processes undergo large interannual to decadal fluctuations due to the El Niño Southern Oscillation (ENSO; Wolter & Timlin, 2011), the Pacific Decadal Oscillation (PDO; Mantua et al., 1997) and the North Pacific Gyre Oscillation (NPGO; Di Lorenzo et al., 2008) that play an important role in regulating year-to-year variability in primary production (Peña et al., 2019; Whitney et al., 2005). The pronounced natural climate variability of this water, makes long-term anthropogenic trends difficult to detect, requiring that time series of observations extend for several decades (Cummins & Masson, 2014).

Fig. 1 Map of the British Columbia continental shelf showing the bathymetry. The red line shows the outer boundary of the British Columbia Continental Margin (BCCM) model, the green line the continental shelf boundary region, the purple and green shaded areas show the coastal regions: northern coast (NBC) and west coast of Vancouver Island (WCVI), respectively. The yellow triangle shows the location of the 6 weather buoys used in Fig. 3, and the magenta lines indicate the location of the transects used in Figs 12 and 13.

Global and regional climate model projections can be used to assess the potential effects of climate change on the marine ecosystem. In most eastern boundary upwelling regions, models project that alongshore winds will intensify under climate change (Bakun et al., 2015; Wang et al., 2015; Xiu et al., 2018), but the effect of this change on primary production differ among these regions depending on phytoplankton nutrient limitation levels and the dynamics controlling the upward flux of nutrients to the coast (Lachkar & Gruber, 2013; Pozo Buil et al., 2021; Rykaczewski et al., 2015). Also, future alteration of the properties of the source water advected into the region might play an important role (Bograd et al., 2019; Rykaczewski & Dunne, 2010). In the California Current System (CCS), earth system models (ESMs) project net primary production to increase in response to nutrient enrichment of upwelling source waters due to a decrease in ventilation of North Pacific Ocean (Rykaczewski & Dunne, 2010). A northward displacement of the North Pacific Current, as it has been projected by Toste et al. (2019), would enhance the transport of cold, nutrient-rich water from higher to lower latitudes, positively affecting multiple trophic levels in the CCS (Malick et al., 2017; Sydeman et al., 2011). Off the Vancouver Island shelf, an ensemble of 18 climate models predicted increased upwelling summer winds in the twenty-first century (Merryfield et al., 2009). However, using a regional ocean circulation model to downscale global projection, Foreman et al. (2014) found that slightly stronger summer winds along the western Vancouver Island shelf do not result in appreciable changes to upwelling. Therefore, although global models have been critical for assessing the likely responses of the marine ecosystem under possible future climate, high-resolution climate models are needed in coastal regions where the coarse resolution of global models do not represent coastal processes adequately (Siedlecki et al., 2021; Xiu et al., 2018).

Whereas wind intensification as a result of climate change could benefit marine populations by increasing nutrient input into the euphotic zone, it could also be detrimental for many organisms by increasing the transport of acidic waters into shelf habitats. Moreover, as the ocean absorbs excessive anthropogenic atmospheric CO₂, ocean acidification leads to a decline in the saturation state (Ω) of seawater with respect to calcium carbonate minerals, such as calcite or the less stable form aragonite (Feely et al., 2016). Ocean acidification has the potential to significantly impact individual marine organisms and ecosystems by reducing the formation of skeletons and shells by marine calcifying organisms (Fabry et al., 2008; Stojkovic et al., 2013;) and affecting larval and adult fish (Munday et al., 2009; Ross et al., 2011). On the west coast of the US, a regional circulation-biogeochemical model has been used to assess the effect of increasing atmospheric CO₂ emission on the CCS that show a rapid progression of ocean acidification from 1995 to 2050 leading to increased intensity, magnitude, and duration of acidified water over the continental shelf (Gruber et al., 2012; Hauri et al., 2013). Additionally for the CCS, dynamically downscaled projections (i.e. using coarse global climate model fields to force a high-resolution regional ocean model) from three ESMs were used to produce climate projections and to assess model uncertainty in the projections (Pozo Buil et al., 2021). In Canadian Pacific waters, regional models have been used to downscale ESMs providing useful information on the potential effect of climate changes under future climate scenarios (Foreman et al., 2014; Holdsworth et al., 2021). However, these high-resolution projections either did not include biogeochemical processes or employed atmospheric climatologies to force the ocean model.

To gain a better understanding of the combined effects of climate change and changes in ocean acidification on conditions over the west coast of Canada, we use a high resolution coupled circulation-biogeochemical model developed and parameterized for the BC continental margin to downscale projections from an ESM. Multi-year simulations were carried out under historical (1981–2010) and future (2041–2070) conditions to produce climate projections under two representative concentration pathways (RCP) 4.5 and 8.5 scenarios (CMIP5; Taylor et al., 2012). We analyse changes in a range of ecosystem-relevant variables (i.e. warming, acidification, and changes in ocean productivity) and the plausible mechanisms associated with those changes. The remainder of the manuscript is organized as follows. Section 2 briefly describes the regional model, the ability of the historical simulation to reproduce observations, and the approach used to downscale the ESM. Section 3 presents projected changes in the ecosystem-relevant variables and then discussion of the results is presented in Sections 4.

2 Methods

a The Coupled Physical-Biogeochemical Model

The model we employ is an implementation of the Regional Oceanic Modeling System (ROMS; Shchepetkin & McWilliams, 2005) coupled to a lower trophic level biogeochemical model for the British Columbia continental margin (BCCM). A detailed description of the BCCM model is given in Peña et al. (2019). The model has a horizontal resolution of 3 km and covers the entire Canadian Pacific continental margin, extending from the Alaska Border (∼55 °N) to south of the Columbia River (∼47 °N), and seaward about 400 km from the shore (Fig. 1). It has 42 terrain-following vertical levels stretched to give high resolution in surface and bottom boundary layers with a surface layer thickness ranging from a minimum of 0.3 m in the shallowest areas (15 m deep) to ∼1.6 m at the deepest depth (c.a. 3800 m deep) and on average about 12 levels within the top 50 m, allowing a reasonable representation of the surface boundary layers. Freshwater discharges from 154 rivers affecting the study area, derived as in Morrison et al. (2012), are implemented as point sources at the coast.

The biogeochemical model, described in detail in Peña et al. (2019), consists of an NPZD (Nutrient-Phytoplankton-Zooplankton-Detritus) module with three nutrient compartments (nitrate, ammonium, and silicate), two phytoplankton groups (diatoms and flagellates), two zooplankton groups (herbivorous zooplankton and omnivorous zooplankton) and three detrital pools (slow-sinking detritus, large fast-sinking detritus, and silicon detritus). Temperature affects all physiological rate model parameters of phytoplankton, zooplankton, and detritus according to Q₁₀ factors, which are set to 2 for phytoplankton, and 3 for zooplankton and remineralization by bacteria. The model was expanded by adding a carbon cycle module that adds three new compartments; dissolved inorganic carbon (DIC), alkalinity (Alk), and mineral CaCO₃ (D_CaCO3), to the biogeochemical model (see Appendix A). The organic carbon transformations are tied to those of organic nitrogen with a fixed Redfield C:N ratio of 106:16. However, different remineralization rates for slow-sinking carbon and nitrogen detritus (DC₁ and DN₁), and for large fast-sinking carbon and nitrogen detritus (DC₂ and DN₂) allow deviations from Redfield C:N ratios (see Table A1). The exchange of CO₂ across the air–sea interface is parameterized following Wanninkhof (2014) and the carbonate chemistry is simulated following the standard Ocean Carbon Cycle Model Intercomparison Project Phase 2 protocols (available at http://ocmip5.ipsl.jussieu.fr/documentation/OCMIP/phase2/simulations/Abiotic/HOWTO-Abiotic-1.html).

b Model Simulations

To examine the potential effects of future climate changes on the west coast of Canada, we use climate projections to create perturbations that were added to the initial conditions and time-varying forcing fields of a realistic multi-year regional hindcast simulation (described below). This dynamic downscaling approach is known as the pseudo-global-warming method (Hara et al., 2008; Kimura & Kitoh, 2007) or the delta method (Alexander et al., 2020) and has been used previously in similar studies (e.g. Foreman et al., 2014). This technique allows a comparison of the climate in the recent past and future years, where the latter is similar to the control years in terms of interannual variations but includes future climatology (Hara et al., 2008). To allow the model to adjust to its lateral boundary and surface forcing, the model was spin-up for 3 years using forcing fields from 1979 repeatedly before running the model with forcing data for 1980 and following years. Including or removing the first 1–2 years of model output did not affect the climatologies considered in this study and downscaled anomalies do not systematically diverge over the time span of the run (Fig. S1). The same spin-up procedure was applied to the future simulations. Average results from the last 30-years of simulations of the coupled circulation-biogeochemical model are presented for the historical (1981–2010) and the future (2041–2070) periods and apart from the initial and forcing fields, use the same model settings.

1 HINDCAST SIMULATION

The 30-year baseline simulation (1981–2010) was previously used to examine the interannual variability of physical and biogeochemical conditions along the BC coast (Peña et al., 2019). Atmospheric conditions (i.e. 3-hourly wind and daily average air humidity, air temperature, air pressure, downward shortwave and longwave radiation at the surface, and precipitation) were obtained from the North American Regional Reanalysis (NARR; Mesinger et al., 2006) and used to compute atmospheric heat fluxes with standard bulk formulae. Tidal boundary forcing was obtained from the barotropic tidal model of Foreman et al. (2000). Initial and monthly open boundary conditions for the ocean fields (temperature, salinity, u and v) were derived from the Simple Ocean Data Assimilation version 2.1.6 results (Carton & Giese, 2008), whereas monthly climatology of nitrate and silicate, repeated over each year of simulation, were obtained from the World Ocean Atlas 2013 (https://www.nodc.noaa.gov/OC5/WOD13/; Garcia et al., 2014). Initial and boundary conditions for DIC and Alk were obtained from the Global Ocean Data Analysis Project (GLODAPv2, Lauvset et al., 2016; Olsen et al., 2016). A seasonal cycle was added to the upper layer DIC and Alk lateral boundary conditions using the monthly climatology of surface ocean pCO₂ from Takahashi et al. (2014) and an observationally based parameterization of alkalinity with regard to sea surface temperature and salinity (Lee et al., 2006), following Hauri et al. (2013). The seasonal cycle of DIC and TA variation were propagated throughout the upper water column by assuming that its seasonal amplitude attenuates proportional to the seasonal amplitude of temperature. Initial and boundary conditions for all other biogeochemical variables were set to low constant values (≤ 0.05 mmol m⁻³). River runoff was defined according to Morrison et al. (2012) and seasonally varying riverine inputs of nitrate and silicate were derived from available observations at the main rivers, as described in Sutton et al. (2013), whereas DIC and Alk were set to 1000 mmol m⁻³ (Ianson et al., 2016). Atmospheric pCO₂ was increased from 342 ppm in 1981 to 391 ppm in 2010 following observations at Cold Bay, Alaska (https://gml.noaa.gov/dv/data/).

Extensive model validation of the hindcast simulation were presented in Masson and Fine (2012) and Peña et al. (2019), which showed that the model can reproduce with reasonable accuracy the important oceanographic features of the BC shelf region. This includes the major seasonal currents, elevation time series from coastal tide gauges and offshore satellite altimetry along-track sites, eddy kinetic energy patterns, the 1997–98 El Niño event, the vertical and cross-shore gradient of density and nutrients, and the mean distributions and seasonal variability of upper layer nutrients and primary production. That the model is able to reproduce seasonal and interannual variability in the past serve as a validation of the models to be used to project impacts of future climate change.

In order to evaluate the ability of the model to reproduce key carbonate chemistry variables, we compare the model results to in situ observations of DIC and Alk, as well as temperature and salinity, from GLODAPv2 and Institute of Ocean Science (Fisheries and Oceans Canada, DFO) surveys at several locations within the model domain and during the study period. The location of in situ observations are shown in Fig. S2. We first assess the ability of the model to reproduce observations within the coastal region (see Fig. 1) and offshore region binned by decade for temperature and salinity but not for DIC and Alk since most observations are from the last decade (>80%). As shown in Fig. 2a, modelled temperature, salinity, DIC, and Alk are well correlated with observations (r ≥ 0.9) in the offshore region, whereas correlations are lower for DIC and Alk (r ∼0.7) and slightly lower for temperature and salinity (r∼0.85) in the coastal region. The model is also able to reproduce the observed variability among regions, with most normalized standard deviation in the rage of 0.8–1.2. Second, we assess the capability of the model to reproduce the vertical structure. In the upper 100 m of the water column, correlations are lower for DIC and Alk (0.63 and 0.58, respectively) than for temperature and salinity (>0.8), but simulated variables are well correlated with observations (r ≥ 0.9) farther below (Fig. 2b). At the same time, the model reproduces the observed variability equally well in both vertical layers with normalized standard deviation in the range of 0.8 to 1.2.

Fig. 2 Taylor diagram of model simulated temperature, salinity, DIC, and TA compared to observations for: (a) shelf (filled symbols) and offshore (open symbols) regions, and (b) 0–100 m (filled symbols) and below 100 m depths (open symbols). The data has been binned by decade: 1981–1990 (red), 1991–2000 (blue), 2001–2010 (magenta), and 1981–2010 (black). Radial distance from the origin is the normalized standard deviation of the modelled parameters. The angle from the vertical represents the correlation between the model and the observations. The distance between the model point and the observation point (x on the abscissa) indicates the normalized root mean square (RMS) misfit between model and observational estimates.

Then, we compared model values to observations from the North American Carbon Program (NACP) West Coast Cruise in May 2007 (Feely et al., 2008) and found a couple of shortcomings. First, the model underestimates Alk in the upper 100 m by 10–30 µmol kg⁻¹ along a transect in the northern shelf and off Vancouver Island (Fig. S3 & Fig. S4). Second, the uplift of isopleths of salinity, DIC, and Alk along the transect off Vancouver Island are not confined to the nearshore in the model but extend farther offshore (Fig. S4), possibly reflecting a poor representation of the cross-shore gradient in alongshore winds due to the coarse resolution of the atmospheric forcing.

2 REGIONAL CLIMATE CHANGE SIMULATIONS

To dynamically downscale climate projections, climatological perturbations were computed as the differences between mean values for the future period (2041–2070) and the historical period (1981–2010) under two climate scenarios: a moderate mitigation (RCP 4.5) and a high emission (RCP 8.5) scenario (Moss et al., 2010). The atmospheric fields to force the regional climate change simulations were constructed by adding daily climatological perturbations (future-historical) from the Canadian regional, atmospheric-only, model (CanRCM4) simulations to the hindcast atmospheric fields during each year of the hindcast simulation. Monthly climatological perturbations (future-historical) from the Canadian Earth System model (CanESM2) simulations were added to the hindcast initial and ocean boundary conditions, repeated over each year of the hindcast simulation, to produce initial and along the boundaries conditions for the regional climate change simulations. In addition, we forced the model with increasing atmospheric pCO₂, as specified by the climate scenario. Like the hindcast simulation, the climate change simulations are 34 year long and the average of the last 30 year of the simulations are presented here. The choice of moving our historical simulation ahead by sixty years is somewhat arbitrary, but it facilitates comparison with earlier studies using a similar lead time (e.g. Foreman et al., 2014; Holdsworth et al., 2021) and it is far enough in the future that most changes can be attributed to changes in green-house forcing rather than internal model variability. At the regional scale, internal variability is still important in global models.

The CanESM2 (Arora et al., 2011) provides ocean initial and monthly lateral open boundary conditions at a horizontal resolution of ∼1.41° longitude and 0.94° latitude and 40 vertical levels, whereas the regional atmospheric only model (CanRCM4; Scinocca et al., 2016), which was driven by the CanESM2 global model, provides all the required daily atmospheric forcing fields at a horizontal resolution of ≃ 22 km to run the model. This model combination was selected because it was the only regional climate model (RCM) driven by an earth system model (ESM) in our model domain available to us and because three-dimensional monthly fields of ocean biogeochemical variables (NO₃, DIC, and Alk) were available from this ESM rather than only annual means as it is obtainable from phase 5 of the Coupled Model Intercomparison Project (CMIP5) archive. The CanESM2 model includes the nitrogen and carbon cycles but not oxygen, so changes in hypoxic conditions are not considered here. The increased resolution of the RCM is able to better capture spatial variations in oceanic winds and heat fluxes in the coastal region (e.g. Morrison et al., 2014) and having three-dimensional monthly temporal resolution of boundary conditions for biogeochemical, as well as physical, variables allows for a better representation of the marked seasonality of biogeochemical properties of the region (Peña et al., 2018). However, even though the CanRCM4 model has a similar spatial resolution than the NARR dataset, as mentioned above, the CanRCM4 projections were not used directly in our simulations because a comparison of CanRCM4 winds over the historical period of 1981 to 2010 with observations at weather buoys along the BC shelf (Fig. 3 & Fig. S5) revealed systematic biases in the magnitude and timing of the spring and fall transitions and, in particular, of the upwelling season. In comparison, NARR reproduces the observed timing and magnitude of summer winds with reasonable accuracy. Thus, employing CanRCM4 projections to directly force our future simulation would produce incorrect upwelling conditions and erroneous conclusions about their impact on future oceanic conditions.

Fig. 3 Mean monthly winds at the buoys location shown in Fig. 1: observed (black), NARR (blue), and CanRCM4 (green) means over the period 1981–2010, and CanRCM4 RCP 8.5 (magenta) mean over the future period 2041–2070.

To illustrate the perturbations applied to the model initial and boundary fields of temperature, salinity, nitrate, and DIC, the monthly mean and latitudinal/longitudinal variable perturbations were average into annual and horizontal mean values (Fig. 4). Under RCP 8.5, the mean perturbations at the boundaries provide warmer and fresher waters that are maximum at the surface by 2.2–2.8 °C and up to −0.4 in salinity and extend from the surface down to about 1000 and 500 m, respectively and are on average similar among the boundaries. In comparison, DIC increases in the upper 1500 m with maximum mean perturbations of 40 to 55 mmol m⁻³ in the 100 to 500 m range and enhanced supply of nutrients between around 50 to 3000 m with maximum mean perturbation at the southern boundary of ∼3.7 mmol m⁻³ at around 500 m. The perturbations under RCP 4.5 have a similar pattern but are less pronounced, reflecting the lower emission scenario (Fig. S6).

Fig. 4 Annual mean profiles for: (a) temperature, (b) salinity, (c) nitrate, and (d) DIC perturbations from CanESM2 under the RCP 8.5 scenario. Red, blue, and magenta lines show the average lateral profile at the southern, northern and western boundaries, respectively. Solid black line is the domain average initial profile, with grey shading representing the range of values.

3 Results

We present projected changes in future conditions as mean changes in model simulated ocean and biogeochemical conditions for the future period (2041–2070) relative to mean conditions of the recent past (1981–2010). Changes between 30-year means are considered to help filter out climate variability in favour of the long-term climate change signal of primary interest herein. We focus on changes in the spatial pattern and seasonal cycle of ecosystem-relevant variables (temperature, nutrient availability, acidity, and phytoplankton biomass and production) on the continental shelf under moderate mitigation (RCP 4.5) and high emission (RCP 8.5) scenarios.

a Changes in Upper Layer Oceanic Conditions

Annual mean sea surface temperature (SST) increases by 1.6 and 2.5°C compared to the historical period over most of the model domain under RCP 4.5 and RCP 8.5, respectively (Fig. 5). In comparison, sea surface salinity (SSS) decreases less nearshore than elsewhere where it decreases by about 0.2 and 0.3 under RCP 4.5 and RCP 8.5, respectively. Changes in SST and SSS result in an increase in upper layer stratification, as indicated by changes in sea water density differences (50 m minus surface) of >0.17 kg m⁻³ everywhere, except in the Strait of Georgia. However, the model resolution is too coarse to allow for an adequate representation of this region (Peña et al., 2016) so results from the Salish Sea and the many narrow channels and inlets that are not adequately resolved by the model will not be considered further. Despite the increase in upper layer stratification, surface nitrate concentration increases in most regions, except on the continental shelf where lower concentrations are found, specially under RCP 8.5.

Fig. 5 Annual mean (1981–2010) historical run (left column) and future changes (2041–2070 relative to historical) under RCP 4.5 (middle column) and RCP 8.5 (right column) for sea surface temperature (SST), sea surface salinity (SSS), σ_t differences (50 m minus surface) as proxies for changes in upper ocean stratification, and sea surface nitrate (NO₃).

Future changes in ocean condition along the open continental shelf seasonally vary (Fig. 6). For example, average SST is projected to increase by 2 to 3 °C during the year under RCP 8.5, with the highest increases occurring in August to October. In comparison, surface waters will become more stratified in the future during spring and summer (May to October) but less stratified during winter months. However, increases in summer stratification along the coast are not strong enough to overcome the effect of upwelling favourable winds, as shown by the increase in surface and subsurface nitrate concentrations in summer and fall. Contrary to expectation, nitrate concentrations decrease in winter when stratification is lower. This is due to a decrease in winter nitrate concentration in the upper 60 m in the future compared to the historical period (Fig. 7), whereas summer nitrate concentration increases below the upper layer in the future.

Fig. 6 Hovmöller diagram of average sea surface temperature (SST), σ_t differences (50 m minus surface), surface nitrate, and nitrate at 50 m depth, as a function of latitude and month in the shelf region for the historical run (top row) and future changes (2041–2070 relative to historical) under RCP 4.5 (middle row) and RCP 8.5 (bottom row).

Fig. 7 Historical and future average vertical profiles of σ_t and nitrate in February (top) and September (bottom) at about 49° N on the west coast of Vancouver Island shelf.

Future changes in surface pH in response to increasing atmospheric CO₂ are smoother and more uniform than SST changes. Annual mean surface pH decreases on average by about 0.14 pH unit under RCP 4.5 and by 0.21 pH units under RCP 8.5 (Fig. 8). The projected increase in annual mean acidity under both scenarios is greater (lower pH) than the minimum pH values simulated for the historical period over most of the model domain.

Fig. 8 Annual mean (top), and minimum (bottom) surface pH for the historical simulation (left column) and future simulations under RCP 4.5 (middle column) and RCP 8.5 (right column). The minimum values were computed from monthly mean values of the 30-year long simulation.

b Changes in Phytoplankton and Primary Production

Projected changes in annual mean total phytoplankton biomass (diatoms + flagellates) are generally small (<15%), spatially heterogeneous, and similar under both climate change scenarios (Fig. 9). Along the continental shelf, chlorophyll concentration decreases on the west coast of Vancouver Island and on the shallow bank on the east coast of Haida Gwaii and increases in Queen Charlotte Sound and east coast of Hecate Strait. The decrease is stronger in the southwest coast off Vancouver Island, the most productive portion of the BC shelf (Ware & Thomson, 2005). Declines in phytoplankton biomass are often linked to nutrient depletion in surface waters, since higher water temperatures generally lead to more stable and stratified water conditions resulting in less nutrient supply from deeper waters. In our case, however, changes in chlorophyll do not coincide with change in surface nitrate concentration at most locations (Fig. 5), particularly in the southwest coast off Vancouver Island. Also, temporal variation of surface and subsurface nitrate along the continental shelf (Fig. 6) shows a decrease in concentrations in winter, when stratification decreases, and an increase in subsurface nitrate in summer during increased stratification. However, the increase in stratification and changes in nutrient environment of the future scenario might be responsible for changes in the type of phytoplankton in favour of diatoms on the shelf and flagellates offshore.

Fig. 9 Annual mean (1981–2010) historical run (left column) and future changes (2041–2070 relative to historical) under RCP 4.5 (middle column) and RCP 8.5 (right column) for surface chlorophyll, fraction of phytoplankton biomass due to diatoms, depth integrated primary production, and fraction of primary production due to diatoms.

Changes in phytoplankton biomass are not necessarily driven by changes in nutrient availability and upper ocean stratification since changes in top-down grazing control and horizontal transport could also be important. In particular, the growth and retention of phytoplankton in the southwest Vancouver Island coast is known to be strongly influenced by the outflow of nutrient-rich water from Juan de Fuca Strait and also by the Juan de Fuca Eddy (Hickey & Banas, 2008). This quasi-stationary cold cyclonic eddy often dominates the summer and autumn circulation pattern of this region (Foreman et al., 2008; Freeland & Denman, 1982), as is evident in the historical simulation panel (Fig. 10). Therefore, changes in its magnitude might be responsible for the projected future decrease in phytoplankton. However, as shown in the seasonal average future and anomaly panels, the eddy circulation is not projected to change in the future, except for a small increase in the current magnitude around the periphery of the eddy in the fall period. This is consistent with the results of Foreman et al. (2014), which do not find any significant changes to future Juan de Fuca Eddies. In comparison, a weaker buoyancy-driven Vancouver Island Coastal Current (Thomson et al., 1989) is projected in winter (January-March) and spring (March-May) periods, but it becomes stronger in the summer (June-August) period.

Fig. 10 Historical (1981–2010), future (2041–2070), and differences (2041–2070 minus 1981–2010) in seasonally averaged currents at 10 m depth in the southwest Vancouver Island shelf under RCP 8.5.

In contrast to chlorophyll, annual mean primary production is projected to increase by as much as 0.12 gC m⁻² d⁻¹, though in the nearshore waters off the east coast of Haida Gwaii a small decrease is projected (Fig. 9). The relative increase is higher (up to 30%) in the shelf-break where primary production is lower than on the continental shelf, where primary production is higher. Seasonally, primary production is projected to increase throughout the year along the continental shelf, except in early spring, with most of the increase occurring during the summer (Fig. 11). The future increase in primary production is mostly due to an increase in regenerated production, although new production also contributes to the increase in summer. Changes in new production are consistent with those of chlorophyll between 50 to 52 °N showing a slight decrease in early spring and an increase in late summer, but not farther south where chlorophyll decreases the most.

Fig. 11 Hovmöller diagrams of average surface chlorophyll and depth integrated primary production, regenerated production, and new production, as a function of latitude and month in the shelf region for the historical run (top row) and future changes (2041–2070 relative to historical) under RCP 4.5 (middle row) and RCP 8.5 (bottom row).

The fraction of primary production that is regenerated within the surface layer depends on food web processes, such as grazing pressure. In the model, both the maximum rate of phytoplankton growth and of zooplankton grazing are temperature dependent. This suggests that temperature and its impacts on the relative growth rates of phytoplankton and their grazers are the main cause of the projected increase in primary production and regenerated production in the model. Warming has been shown to enhance phytoplankton growth rates; however, this is strongly modulated by nutrient availability, suggesting that under more stratified conditions warming may have a negative effect on phytoplankton growth (e.g. Marañon et al., 2014). Previous studies have shown that coastal waters off Vancouver Island are rich in nutrients (Hickey & Banas, 2008; Mackas et al., 1980; Peña et al., 2019), so phytoplankton growth in this region is less limited by nitrate than farther south in the CCS. Since light rather than nutrient mostly limits phytoplankton growth in this region, the future increase in stratification would allow phytoplankton to remain near the surface where favourable light conditions sustain phytoplankton growth.

c Changes in the Vertical Structure Across the Continental Shelf

To evaluate changes in the vertical structure of physical and biogeochemical tracers, mean winter and spring cross sections at two locations, Queen Charlotte Sound (QCS) and west coast of Vancouver Island (WCVI), are presented (see Fig. 1). We also present changes in the distribution of the aragonite saturation state (Ω_arag) as this property influences habitability for calcification by marine organisms. The values of Ω_arag are calculated from model dissolved inorganic carbon (DIC) and total alkalinity (Alk) using carbonate system equilibrium equations (CO2SYS) (https://www.nodc.noaa.gov/ocads/oceans/CO2SYS/co2rprt.html).

The QCS transect is characterized by deepening of the thermocline near the coast by downwelling in winter and levelling of the thermocline in summer due to relaxation of downwelling (Fig. 12). The deepening of isopleths of nitrate, pH, and Ω_arag follows that of the isopycnal surfaces. Upper layer nitrate concentrations are higher in winter than in summer while pH and Ω_arag are lower in winter than summer. In the future, similar changes occur under both climate change scenarios but are stronger under RCP 8.5 (Fig. 12) than RCP 4.5 (Fig. S7). There is an overall increase in temperature of ∼ 0.5 °C at 300 m that is surface intensified and stronger in winter than in summer, increasing upper layer stratification and reducing the flux of nitrate to the upper layer, especially in winter. Below the upper layer, nitrate concentration increases in winter and summer. pH and Ω_arag decrease at all depths in the upper 300 m with changes being more pronounced in summer than winter. Changes in Ω_arag are higher near the surface (−0.4 and −0.6 in winter and summer, respectively), whereas pH decreases the most just below the surface at ∼ 30 m and by up to 0.18 and 0.23 in winter and summer, respectively.

Fig. 12 Average winter (Dec-Jan-Feb) and summer (Jun-Jul-Aug) contours of density, nitrate, pH, and aragonite saturation state along a transect across Queen Charlotte Sound (Fig. 1) for the historical run (left columns) and future changes (2041–2070 relative to historical) under RCP 8.5 (right columns).

Future changes in all properties, except temperature, are stronger in the southern shelf transect (Fig. 13) than in the norther shelf one. The southern transect is influenced by summer coastal upwelling and the near-shore Vancouver Island Coastal Current (VICC), which provides a source of nutrients and forms a barrier that prevents upwelled waters rich in C from penetrating the inner shelf (Bianucci et al., 2011). Upwelling regions have naturally a lower pH, which cause a decline in the saturation state of seawater with respect to CaCO3 minerals, such as calcite and aragonite. This is because the upwelled waters are enriched in CO₂ from the remineralization of organic matter in the ocean interior, and thus have low pH and Ω. Along this transect, changes in nitrate, pH and Ω_arag are more pronounced just below the surface between 15–40 m in summer. Under RCP 8.5, the summer aragonite saturation state is predicted to fall below saturation (Ω_arag < 1) in waters deeper than 40 m. In winter, nitrate concentrations decrease and changes in pH and are less pronounced in the nearshore region within 20 km of the coast than farther offshore under both climate scenarios (Fig. 13 and Fig. S8).

Fig. 13 Average winter (Dec-Jan-Feb) and summer (Jun-Jul-Aug) contours of density, nitrate, pH, and aragonite saturation state along a transect off the mid-Vancouver Island shelf (Fig. 1) for the historical run (left columns) and future changes (2041–2070 relative to historical) under RCP 8.5 (right columns).

The BC coast is characterized by shallow depth horizons for carbonate saturation, so undersaturated waters are present well within the continental shelf (Fig. 14). In particular, carbonate concentrations are commonly undersaturated (Ω_arag <1) below 100 m depth, and thus inhibit shell formation by calcifying organisms. The Ω_arag saturation horizon is projected to rise up in the water column by approximately 10–50 m and 20–70 m under RCP 4.5 and RCP 8.5, respectively (Fig. 14). The isopleth shoals more in the troughs of the northern region, where the saturation horizon is deeper, than off the west coast of Vancouver Island. The projected increased shoaling of the aragonite saturation horizon would have a detrimental impact on calcifying organisms such as pteropods, though biological sensitivities may begin at higher thresholds.

Fig. 14 Average depth of the threshold aragonite saturation state (Ω_arag= 1) in summer along the continental shelf for the historical run (left plot) and future changes (2041–2070 relative to historical) under RCP 4.5 (center plot) and RCP 8.5 (right plot). Negative values indicate a shoaling of the aragonite saturation horizon.

As a metric of variability in these habitat constraints, we computed the volume of water subject to a variety of omega saturation levels. Water volumes are computed as the sum of grid cell volumes with a given range of Ω_arag that are on the BC continental shelf region (Fig. 1). The volume of Ω_arag varies seasonally, particularly in the upper layer due to the strong influence of phytoplankton productivity in spring-summer that depresses the partial pressure of CO₂, resulting in these waters being Ω_arag supersaturated at present (Fig. 15). In the future, the Ω_arag is projected to drop rapidly, with much of the shelf developing summer-long undersaturation in the 50–100 m depth. While today the water masses in the top 50 m of the continental shelf are nearly always supersaturated, on average 6 and 12% of these waters will become undersaturated under RCP 4.5 and 8.5, respectively. Under both climate scenario, waters with Ω_arag above 1.5 will disappear below 50 m, water with Ω_arag above 2 will mostly disappear in the top 50 m, and the bottom layer of the shelf will be below the saturation state of aragonite for 2–3 months in mid-summer.

Fig. 15 Monthly mean volume of seawater with a particular Ω_arag in the continental shelf region from the historical (left column), RCP 4.5 (center column), and RCP 8.5 (right column) runs in the upper 50 m (top row), in between 50 and 100 m (middle row), and in the bottom layer of the model (bottom row). Error bars ( = 1std) represent the interannual variability for each month.

4 Discussion

We use a regional, high-resolution coupled circulation-biogeochemical model (BCCM) specifically developed and parameterized for the BC continental margin, to provide climate change projections of physical and biogeochemical conditions under two future climate scenarios. To do this, we dynamically downscale projections from the CanRCM4 and CanESM2 models under RCP 4.5 and RCP 8.5 using the pseudo-global-warming method, which has been used previously in similar studies (e.g. Foreman et al., 2014). The RCPs are labelled according to the additional radiative forcing level in 2100 with RC P4.5 and RC P8.5 CO₂ concentrations reaching 538 and 936 ppm and representing a moderate mitigation and a high emission climate scenario, respectively (Taylor et al., 2012). We found similar patterns of changes in the future (2041–2070) relative to the recent past (1981–2010) under both climate change scenarios considered, but projected changes were stronger and more extensive under RCP 8.5. For example, changes in surface pH and temperature are substantially higher (∼70%) under RCP 8.5 than RCP 4.5 showing the benefit of mitigation efforts in reducing the effect of climate change. Analogous new scenarios and ESMs simulations are now available via CMIP6 but are unlikely to substantially modify the main results of this study.

We found that sea surface temperature increases by 1.6 and 2.5°C over most of the model domain under RCP 4.5 and RCP 8.5, respectively. These changes are similar to the warming of 1.8–2.6 °C projected by CanESM2 model over the same future period, and of 1.8–2.4°C projected by the regional NEP36-CanOE model of Holdsworth et al. (2021; available at https://open.canada.ca/data/en/dataset/a203a06d-9c1f-4bb1-a908-fc52912ff658.) in the future period of 2046–2065 under RCP 4.5 and RCP 8.5, respectively (Fig. 16). These changes are ∼20% higher than the range of warming of 1.4–2 °C projected by Foreman et al. (2014; their Fig. 2) for the 2065–2078 period under Special Report on Emissions Scenarios (SRES) A2 emission scenario. Also, the direction of future changes agrees between these models (Fig. 16), projecting sea surface freshening, and an increase in upper layer stratification, primary production and acidification under both climate scenarios. However, the magnitude of projected changes between the regional models are in general more similar and lower than those from the CanESM2, except for surface nitrate. These differences between the ESM and regional models are more marked in the coastal regions, especially on the NBC region. Despite the similarities in the mean values between the regional models, there are significant differences in the spatial patterns of change between these two models (Fig. S9).

Fig. 16 Annual mean anomalies (future – historical) under RCP 8.5 for: (a) seas surface temperature (SST), (b) sea surface salinity (SSS), (c) upper ocean stratification (σ_t at 50 m minus surface), (d) sea surface nitrate (NO₃), (e) depth integrated primary production (PP), and (f) surface pH over the model domain (MD), north BC coast (NBC) and west coast of Vancouver Island (WCVI). Future changes are presented for the CanESM2 (2041–2070 minus 1981–2010), BCCM (2041–2070 minus 1981–2010) and NEP36-CanOE (2046–2065 minus 1986–2005) models. For comparison between models, the ESM2 and NEP36-CanOE values were interpolated to the BCCM model grid prior to averaging the fields. Error bars ( = 1std) represent the spatial variability in each region.

Warming stimulated phytoplankton growth in our model due to the direct effect of temperature on the maximum growth rate of phytoplankton. However, as shown by this and previous studies (e.g. Taucher & Oschlies, 2011), a better understanding and model representation of direct temperature effects on metabolic processes is important for reliably projecting future changes in productivity. In particular, differences in temperature sensitivity between phytoplankton growth and zooplankton grazing should be studied further given their important role controlling phytoplankton biomass under warming conditions (e.g. O'Connor et al., 2009). We found that projected increases in temperatures and enhanced stratification result in increased primary production by up to ∼25% in the continental shelf consistent with previous studies that found light rather than nutrient is more limiting to phytoplankton growth in this subarctic coastal region (Peña et al., 2019) and also consistent with observations that showed surface nitrate on the BC shelf is rarely depleted to undetectable levels (Mackas et al., 1980). The increase in simulated primary production was not in response to increase in upwelling or nitrate concentration of the source water with climate change, as it has been found farther south in the CCS (Bakun et al., 2015; Rykaczewski & Dunne, 2010), but due to enhanced recycling of nutrients with warming as indicated by the increase in regenerated production. There was also a moderate decrease (< 20%) in phytoplankton biomass at most of the coast in the future, likely as a result of increase in grazing pressure with warming since temperature affects all metabolic rate parameters in the model. Similar moderate declines in coastal chlorophyll in the CCS under RCP 8.5 scenario were reported by Pozo Buil et al. (2021) in two of their three downscaled global models analysed. The decrease in chlorophyll was explained by enhanced stratification in combination with weakened upwelling-favorable winds. In comparison, Xiu et al. (2018) found an increase in nutrient fluxes to the euphotic zone with climate change in the CCS, but the phytoplankton response to this increased nutrient input was complex and nonlinear due to zooplankton grazing being an important factor controlling phytoplankton biomass.

The implications of future increases in primary production for living marine resources remain uncertain. Populations of pelagic fish and zooplankton feeding near the surface may benefit from increased primary production. However, mid-water and benthic organisms might be unable to tolerate more frequent low-oxygen and corrosive events driven by local aerobic respiration triggered by the increase in productivity, which consumes oxygen and produces CO₂ in approximate stoichiometric equivalence (Anderson & Sarmiento, 1994). Understanding and projecting future patterns of primary production is, therefore, important not only to estimate potential impacts of climate change on marine ecosystem and fishery yields but also on C and O₂ content.

The projected increase in primary production along the BC continental margin is matched by other responses with potentially deleterious consequences. In particular, surface pH decreased by 0.14 and by 0.21 over most of the model domain under RCP 4.5 and RCP 8.5, respectively. Moreover, significant declines in the mean saturation state of aragonite in the future period were found as a result of the increased concentration of CO₂ in the atmosphere. Since the preindustrial period, uptake of increased CO₂ has caused the Ω_arag horizon to shoal by approximated 30–50 m along the west coast of North America (Feely et al., 2016) and the percent of the year spent above corrosive conditions has dropped by ∼40% in the Northern Salish Sea (Evans et al., 2019). In our simulations, the Ω_arag horizon is projected to rise in most of the shelf by at least 20 m in 60 years, significantly increasing the volume of corrosive waters on the shelf. Similar up lift of Ω_arag and increase of corrosive waters in the shelf region has been projected in the CCS (Hauri et al., 2013). Declines in Ω_arag will pose increasing physiological challenges to calcifying invertebrates such as pteropods, bivalves and echinoderms (Stojkovic et al., 2013). Present day surface pH of BC coastal waters exhibits considerable spatial and temporal variability (Feely et al., 2008). This variability could provide a certain degree of adaptability by marine organisms in the regions to future ocean acidification, but only as long as future changes do not move substantially out of this natural range. However, the annual mean acidity projected by the 2050s under both scenarios was found to be greater (lower pH) than the minimum pH values observed over the historical period.

Local drivers in coastal systems, such as river input, eutrophication, and increased organic matter degradation, will also cause more intense acidification effects. Because future changes in these drivers were not available, notably river discharge data for the future climate simulations, river flow and its nutrient load were unchanged from the historical simulation. Therefore, changes in surface freshwater flux are only due to changes in evaporation and precipitation. Foreman et al. (2014) found that under SRES A2 emission scenario, more precipitation over the watersheds emptying into coastal waters produces a stronger estuarine flow in Juan de Fuca Strait and a stronger VICC. In contrast, we project a decrease in the VICC in winter and spring and an increase in fall under unchanged estuarine outflow. Given this limitation, the projected changes in VICC is considerably uncertain.

The effects of climate change on the marine ecosystem are difficult to predict due to the complex interactions among marine organisms, and between them and abiotic forces. The future regional projections presented here are representative of possible changes in conditions under the RCP 4.5 and RCP 8.5 emission scenario, but selecting a single model among the entire IPCC model ensemble does offer a somewhat limited viewpoint on regional climate change. Frölicher et al. (2016) evaluated the relative importance of the sources of uncertainties in future projections of potential ocean ecosystem stressors from ESMs. They found that scenario uncertainty is the dominant sources of uncertainty at the global-mean scale for pH and that model uncertainty is comparable to or larger than scenario uncertainty for most other variables considered in the mid-century, whereas internally generated climate variability is the dominant source of uncertainty in most coastal marine ecosystems over the next few decades. Less is known about the relative importance of different uncertainty sources of dynamically downscaled future projections. Cheng et al. (2021) dynamically downscaled projections from three ESMs under two climate scenarios in the Eastern Bering Sea shelf and found that ESM model uncertainty was larger than scenario uncertainty. Pozo Buil et al. (2021) downscaled projections from three ESMs under RCP 8.5 for the CCS and found that while trends in the downscaled models reflect those in the ESMs that forced them, the downscaled biogeochemical variables differ more to the ESM than physical variables. In addition, several methods are being used to dynamically downscale projections and the sensitivity of the results to these different methods have not been assessed yet. For example, several studies have used the delta method (e.g. Siedlecki et al., 2021), while Pozo Buil et al. (2021) use a time-varying delta method, Holdsworth et al. (2021) use atmospheric climatologies with augmented winds, Cheng et al. (2021) use direct output of the ESM to drive the regional simulation, and Hauri et al. (2013) forced their regional model with increasing atmospheric pCO₂ and increasing DIC at the lateral boundaries as projected by an ESM and the physical part of the model remined the same. Our study used the delta method and explored scenario uncertainty (RCP 4.5 and RCP 8.5), but due to the computational cost of dynamical downscaling, it was not possible to downscale the output from a range of ESMs in order to reduce any biases intrinsic to individual models. However, in contrast to downscaled projections from other studies, we use atmospheric forcing fields from a regional climate model and monthly-mean boundary conditions of biogeochemical and physical variables. Given the importance of wind forcing in upwelling regions and of adequate representation of changes in the upwelling source water, these improvements likely provide more accurate projections of future conditions.

Projections of how physical and biogeochemical conditions of the ocean ecosystem might change under possible future climate are needed to manage marine ecosystems and plan potential strategies for adaptation. Dynamic downscaling of ESMs plays an important role in addressing this need by providing regional projections of change in multiple ecosystem stressors, such as acidification, warming, deoxygenation, and changes in ocean productivity. These regional model projections have provided the physical and biogeochemical foundation for projections of higher trophic levels (Thompson et al., 2023a; Thompson et al., 2023b) and of changes on ecological connectivity among marine protected areas (Friesen et al., 2021) in BC waters demonstrating the utility of high-resolution projections for regional marine resource applications. Although downscaling exercises are being conducted with more than one ESM to reduce model uncertainty, other sources of uncertainty remain unexplored, such as the influence of downscaling methods on biogeochemical variables. However, running multiple simulations is computationally expensive, so efforts to better understand and identify key drivers of change are also important. Our results of the importance of the direct effect of temperature on metabolic processes for projecting future changes in productivity serve to emphasize the continued importance of improved process representation and improved model parameterization. More work is also needed to determine the capacity for adaption to these potential future changes. New observations and advances in bringing together models and observations can be expected to deliver significant progress in these areas.

Acknowledgements

This work was funded by the Aquatic Climate Change Adaptation Services Program of Fisheries and Oceans Canada (DFO). We are grateful to all the staff at DFO who collected and maintained field measurements, the Canadian Centre for Climate Modelling and Analysis for CanESM and CanRCM data, and the providers of the public-domain datasets used in this study. We thank Wendy Callendar for assisting preparing boundary conditions and Samantha Huntington for compiling biogeochemical data on the BC coast. We also thank two anonymous reviewers for their careful attention and constructive comments that helped improve this manuscript and internal reviewers Mike Foreman, Jennifer Jackson, and Bill Crawford for their useful suggestions.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplemental data

Supplemental data for this article can be accessed at https://doi.org/10.1080/07055900.2023.2239186.

Additional information

Funding

This work was supported by funding from Fisheries and Oceans Canada (DFO).

Appendix – Carbon biogeochemistry module

The inorganic ocean carbon cycle adds 3 new compartments to the biogeochemical model: dissolved inorganic carbon (DIC), Alkalinity (Alk), and mineral CaCO₃ (D_CaCO3), which are subject to physical transport and mixing.

Local changes in DIC occur due to CaCO₃ formation and dissolution, net primary production and respiratory processes, which are coupled to the ecosystem model, and gas exchange at the air-sea interface. $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }DIC}{\mathrm{\partial }t}=-r_{CN}({\mu }_S{P}_S+{\mu }_L{P}_L)-r_{CN}{f}_{CaC{O}_3}{\mu }_S{P}_S}\\ & +r_{CN}(r_1{Z}_1+r_2{Z}_2)+{\delta }_{C1}D{C}_1+{\delta }_{C2}D{C}_2+r_{CN}{\delta }_S{S}_D\\ & +k_{D_{CaC{O}_3}}{D}_{CaC{O}_3}+k_{S_{CaC{O}_3}}{S}_{CaC{O}_3}+F_{C{O}_2}\end{array}$where µ_SP_S and µ_LP_L are the net primary production of phytoplankton (flagellates: P_S and diatoms: P_L) controlled by light and nutrients (nitrate and ammonium for all phytoplankton and also by silicate for diatoms); $f_{CaC{O}_3}$ is the fraction of P_S production due to calcifiers; r₁Z₁ and r₂Z₂ are the respiration of microzooplankton (Z₁) and mesozooplankton (Z₂); δ_C1DC₁ and δ_C2DC₂ are the water column remineralization of slow-sinking carbon detritus (DC₁), which behave similarly to dissolved organic carbon, and large fast-sinking carbon detritus (DC₂); δ_SS_D is the sediment detritus (S_D) remineralization; $k_{D_{CaC{O}_3}}{D}_{CaC{O}_3}$ and $k_{S_{CaC{O}_3}}{S}_{CaC{O}_3}$ are the dissolution of CaCO₃ detritus in the water column (D_CaCO3) and sediments (S_CaCO3); and F_CO2 is the gas exchange flux of CO₂. The carbon fluxes are tied to those of nitrogen through the constant C:N Redfield ratio (r_CN). However, different remineralization rates for slow-sinking carbon and nitrogen detritus (DC₁ and DN₁), and for fast-sinking carbon and nitrogen detritus (DC₂ and DN₂), where nitrogen is recycled faster than carbon, allow deviations from Redfield C:N ratios (see Table A1).

The bottom sediment layer is parameterized as a simple one-layer model added to the bottom layer of the model domain. For a detailed description of the biogeochemical model the reader is referred to Peña et al. (2019).

The sea-to-air CO₂ flux, $F_{C{O}_2}$, is computed using the following equation: $F_{C{O}_2}=k{p}_{C{O}_2}$where α is the solubility of CO₂ in sea water computed from Weiss (1974), k is the gas transfer velocity based on Wanninkhof (2014), and ΔpCO₂ is the difference between the partial pressures of sea surface and atmospheric CO₂. The direction of the air-sea flux is determined by the gradient of the partial pressure, where a negative flux indicates CO₂ uptake by the ocean, and a positive flux outgassing. Oceanic pCO2 is calculated from DIC, Alk, temperature, and salinity using the standard OCMIP carbonate chemistry routines (available at http://ocmip5.ipsl.jussieu.fr/documentation/OCMIP/phase2/simulations/Abiotic/HOWTO-Abiotic-1.html).

Alkalinity is unaffected by air-sea exchange of CO₂ but changes due to new production and nitrification, and by the production and dissolution of CaCO₃ detritus in the water column and sediment, such as: ${\displaystyle \frac{\mathrm{\partial }Alk}{\mathrm{\partial }t}=r_{CN}({\mu }_S^{new}{P}_S+{\mu }_L^{new}{P}_L)+2k_{D_{CaC{O}_3}}{D}_{CaC{O}_3}+2k_{S_{CaC{O}_3}}{S}_{CaC{O}_3}-2r_{CN}{f}_{CaC{O}_3}{\mu }_S{P}_S-\phi N{H}_4}$Where ${\mu }_S^{new}{P}_S$ and ${\mu }_L^{new}{P}_L$ are the production of phytoplankton due to nitrate uptake (i.e. new production) and φNH₄ is the formation of nitrate (i.e. nitrification), which is light limited.

The concentration of mineral CaCO₃ (mmol C m⁻³) depends on net primary production of small phytoplankton and dissolution: ${\displaystyle \frac{\mathrm{\partial }{D}_{CaC{O}_3}}{\mathrm{\partial }t}=r_{CN}{f}_{CaC{O}_3}{\mu }_S{P}_S-k_{D_{CaC{O}_3}}{D}_{CaC{O}_3}-w_C{\displaystyle \frac{\mathrm{\partial }{D}_{CaC{O}_3}}{\mathrm{\partial }z}}}$The formation of D_CaCO3 is parameterized as 15% of net primary production of small phytoplankton such that for each mol of organic carbon formed by this phytoplankton group, 0.15 mol of D_CaCO3 are formed.

Table A1. Summary of the parameters used in the carbon module of the biogeochemical model.

Symbol	Description	Value	Units
rCN	Carbon to nitrogen ratio of organic matter	6.625	mol C: (mol N)^−1
$f_{CaC{O}_3}$	Fraction of PS production due to calcifiers	0.15	–
r1	Basal metabolism rate of Z1	0.075	d^−1
r2	Basal metabolism rate of Z2	0.05	d^−1
δC1	Remineralization rate of DC1	0.03	d^−1
δN1	Remineralization rate of DN1	0.05	d^−1
δC2	Remineralization rate of DC2	0.01	d^−1
δN2	Remineralization rate of DN2	0.03	d^−1
δS	Carbon remineralization rate of DS	0.003	d^−1
$k_{D_{CaC{O}_3}}$	CaCO3 dissolution rate	0.0057	d^−1
$k_{S_{CaC{O}_3}}$	Sediment dissolution of CaCO3	0.002	d^−1
wC	Sinking velocity of $D_{CaC{O}_3}$	20.0	m d^−1