Phenology metrics for ocean waters with application to future climate change in the Northwest Atlantic Ocean

Phenology has evolved from studies of periodic events in biological life cycles to include cycles of physical variables and how changes in the latter affect ecosystem components. With the increased interest in climate change as a driver of shifts in the timing of seasonal events, the number of published papers related to phenology has proliferated over the past three decades. A Web of Science (WoS) search using the term phenology returned about 1,250 papers for the 1990–1994 pentad, growing to almost 10,000 for the 2015–2019 interval (see Section S1 for more details). Notably, ocean-related papers were far fewer, with the search term “phenology and ocean” returning just 5 for 1990–1994, increasing to 423 for 2015–2019.

There have been several review/synthesis studies of phenology in the ocean from which some common themes have emerged (see, e.g., Edwards and Richardson, 2004; Koeller et al., 2009; Rijnsdorp et al., 2009; Ji et al., 2010; Burrows et al., 2011; Mackas et al., 2012; Poloczanska et al., 2013; Poloczanska et al., 2016; Staudinger et al., 2019; Ardyna and Arrigo, 2020; Zhao et al., 2022). One theme is that phenology-related shifts are likely faster in the ocean than on land. Evidence for this difference is summarized by Poloczanska et al. (2013) who, in a meta-analysis of the literature, reported an oceanic seasonal advancement approximately double that reported on land (4.5 vs. 2.3 days decade⁻¹) and a range expansion rate for marine species about 10 times that of terrestrial species (70 vs. 6 km decade⁻¹). Similarly, Asch (2015) reported mean phenology shifts of 4.7 days decade⁻¹ (earlier) for larval fish species in the California Current ecosystem, also about twice the rate as on land. Another recurring theme is that trophic mismatches are expected as climate warms and the phenology of predator/prey species changes. This expectation is highlighted in, for example, Edwards and Richardson (2004) who examined planktonic data from 1958 to 2002 and Mackas et al. (2012) who studied zooplankton phenology in the northern hemisphere ocean.

While many of the reviews of ocean phenology consider the pan-oceanic to global scale, there are important studies focused on the area of interest of this study, the Northwest Atlantic Ocean. Staudinger et al. (2019) synthesized the literature on phenology in the Gulf of Maine (GoM) region, providing estimates for changes in physical metrics (spring and fall thermal transitions, stratification onset), primary and secondary productivity (shifts in phytoplankton bloom and zooplankton abundance metrics), as well as metrics for fish, macroinvertebrate, and seabird changes (see their table 1). Record et al. (2019) investigated phytoplankton phenology in the GoM, finding a notable shift toward later timing for the spring and fall bloom on the century timescale, although Friedland et al. (2023) reported that this shift does not hold when considering only the last 20 years of data. The influence of changing seasons has also been observed at higher trophic levels. For example, Pendleton et al. (2022) documented the relationship between the time of peak habitat use and the (changing) onset of spring conditions for fin and right whales in the Cape Cod region of the GoM, with implications for resource management, and Langan et al. (2021) reported changes in the migration timing of 12 fish species in Narragansett Bay (Rhode island).

Commonly, ocean phenology studies focus primarily on the present climate. This focus is supported by our WoS literature review (see Section S1), which shows that overall about 50% of studies mention “climate change” but less than 3% mention “future climate” (for recent examples of the latter, see Asch et al., 2019; Asch et al., 2022; Yamaguchi et al., 2022; Kléparski et al., 2023). Furthermore, less than 5% mention standard phenology metrics such as the change in the onset of spring or the duration of summer. This analysis invites two questions. First, are there other phenology metrics of relevance to ecosystem changes in the ocean? And second, what are the projected changes in phenology under future climate change?

This article presents an expanded set of phenology metrics, focusing on how they may change in a future climate. The metrics are based on ocean physical variables that can be derived from all future climate ocean model simulations, with each metric expressing the change in the variable of interest between the future and present climates. Our set of phenology change metrics is not meant to be exhaustive. Rather, it is a proposed starting point for use in a preliminary assessment of how the phenology of an area is likely to change and to illustrate the types of ecosystem-related calculations that can be made using future climate ocean model output. Indeed, these metrics can be used broadly over a range of scientific questions, so we refer to them interchangeably as phenology or climate change metrics.

An additional focus of this study was to consider how the future climate change metrics can be applied to predict changes in the ecosystem. To do so, we searched the literature for relationships between phenology variables and ecosystem properties such as species growth rates and thermal habitat. We demonstrate that some of our climate change metrics can be applied directly to predict a phenological change in a species, while some can be used indirectly to predict a change that could lead to a phenological shift. Some metrics can also be used to predict non-phenology-related changes to ecosystem properties. These applications are presented to illustrate the utility of the metrics in predicting possible future climate ecosystem changes.

The suite of phenology metrics derived in this study were computed using present and future climate simulations from the Bedford Institute of Oceanography North Atlantic model (BNAM; Brickman et al., 2016; Wang et al., 2018). The calculations were made on a grid cell basis, producing spatial maps of the phenology change metrics. The horizontal resolution of BNAM is such that it can illustrate the spatial complexity of these calculations on the shelf regions of Atlantic Canada (Figure 1). We considered that using a single ocean model was acceptable because we are interested primarily in the metrics themselves and their applications, as opposed to the variability in their values.

In the following, we present the proposed set of phenology metrics in Section 2 and the details of the ocean model simulation in Section 3. Section 4 contains the results of the future climate projections of phenology change. Examples of how the phenology metrics can be used to predict changes in ecosystem properties in the future are presented in Section 5. We conclude with a summary of our study in Section 6.

In this section we present a set of eight species-independent physically based phenology metrics for ocean waters. Three have been reported in the literature, while the remaining five are new. For each metric we present basic information on how it is calculated (the variable(s) on which it is based and its units) and briefly describe how it relates to ecosystem function. Detailed explanations on how the metrics are calculated are presented in Section S2. Because our interest is in future phenology changes, the phenology metrics and how they are calculated reflect differences between a projected future climate and the present climate.

There are three phenology metrics commonly reported in the (ocean) literature and in the popular press: change in the onset of spring, change in the length of the growing season (or summer), and change in the onset of stratification. These metrics are described in the following sections.

The onset of spring is linked to the emergence and development of numerous terrestrial and marine plants and animals, as well as the arrival times of migratory species (see, e.g., Inouye, 2022 and references therein). Differences in the way species react to changes in the onset of spring can increase the potential for predator-prey mismatches, with negative impacts on individual fitness, population dynamics, and ecosystem function (Staudinger et al., 2019).

The onset of spring in the ocean, when based on physical variables, is typically computed as the arrival time of a particular sea surface temperature (SST) isotherm (see, e.g., Friedland et al., 2015; Wang et al., 2021a), and changes to this arrival time have been calculated from present climate time series (Galbraith and Larouche, 2013; Thomas et al., 2017). Here we calculated when the present climate April SST value occurs in the future and report this result as a time difference (in weeks; see Section S2.2).

There are various calculations for the length of summer in the ocean literature (Galbraith and Larouche, 2013; Pershing et al., 2015; Henderson et al., 2017; Thomas et al., 2017), reflecting the fact that both a calendar-based definition of summer and a longer warm or temperature-based growing season can be of ecological significance. Here we considered the broader temperature-based definition and computed the length of the growing season.

The importance of changes to the length of the growing season has been identified by the Intergovernmental Panel on Climate Change (IPCC) as far back as 2007 (IPCC, 2007). Changes to the length of the growing season have been linked to, for example, an increase in global terrestrial net primary productivity (Piao et al., 2007), cod mortality (Pershing et al., 2015 and see below), and shifts in both distribution and biomass of fish stocks in the GoM (Henderson et al., 2017). Henderson et al. (2017) also noted that a longer growing season “may improve growth conditions, alter maturity schedules and improve reproductive success of some stocks.”

The calculation of the length of the growing season typically uses the departure time of a particular late-year SST isotherm in concert with the onset of spring calculation to derive a duration of the growing season. To accommodate the large range of temperatures in our geographic domain, we defined the length of the growing season based on the present climate mean April and October SSTs, and determined when these two values would occur in the future. The result is the future climate length of the growing season. We report it as a time difference from the present climate length (April–October) in weeks (Section S2.2).

The onset of stratification is often linked to the initiation of the spring phytoplankton bloom, whereby the combination of the development of stratification, which serves to retain nutrients and phytoplankton in the upper mixed layer, and increasing illumination allows accelerated phytoplankton growth to occur (the Sverdrup “Critical depth hypothesis”; Sverdrup, 1953). While this idea has formed the basis for understanding bloom dynamics, further work has indicated that the causal effect of stratification on the initiation of the bloom is not so apparent (Townsend et al., 1992; Behrenfeld, 2010; Lindemann and St. John, 2014).

Here we used the common definition of stratification as the difference in density between a depth of 50 m and the surface. The onset of stratification was computed as the time when a particular value occurs, the latter based on, for example, the annual mean stratification (Staudinger et al., 2019). For the GoM region, spatial patterns in the timing of peak stratification and decadal changes to the onset of stratification have been reported for the present climate (Li et al., 2015; Staudinger et al., 2019). In this study we computed the present climate annual mean stratification, at each grid cell, and determined when this value occurs in the present and future climates. The change in onset is the difference in these two times (in weeks; see Section S2.2).

We propose five new phenology metrics based on physical variables. These metrics relate directly to changes in the timing of particular events in the future or measure a difference in a physical quantity that can result in a phenology change. Two of the five metrics measure the period of time that the future climate spends above a threshold temperature in the present climate as a general thermal stress indicator. The other three are based on degree day (DD) calculations. Degree days (units of °C × days) are defined as the integral of temperature over time: $DD=∫t0t1Tdt$⁠, where T is the temperature, t is the time, t0 is the start time, and t1 is the end time. Two of these DD-based metrics measure when a given present climate number of DDs “arrive” in the future, providing an indicator of the speeding up of ecosystem development times in the future. The last metric is the difference in DD at a given time of year.

We define this metric as the period of time (in weeks) that the future climate spends above the maximum present climate SST or above the maximum present climate bottom temperature (Tbtm). Because the maximum temperature that an organism experiences in the present climate may approach its thermal tolerance limit, the period of time in the future that is spent above this value can indicate potential thermal stress. This concept is related to “thermal safety margin,” which measures the difference between observed temperature and the thermal tolerance limit of a species (Boyce et al., 2022). Details of this calculation are provided in Section S2.1.

The planktonic development index refers to the difference between present and future climate arrival times of April–September SST DDs. DDs are typically related to development times for various life stages (Neuheimer and Taggart, 2007; Neuheimer et al., 2008; Neuheimer and MacKenzie, 2014). For species with floating life stages, DDs calculated for the surface layer can be considered to be a larval development metric. Similarly, this metric could apply to planktonic or pelagic species. For this arrival time calculation, the number of present climate SST DDs from April to September was computed for each grid cell. Then, in the future climate, starting in April, we computed when the same number of DDs “arrives.” The arrival time is the difference, in weeks (Section S2.2).

The demersal development index refers to the difference (in weeks) between present and future climate arrival times of January–December and October–April Tbtm DDs. For bottom-dwelling organisms, Tbtm DDs can be related to growth and to the timing of a developmental stage (e.g., spawning). Generally, bottom DDs are computed for the calendar year, but different periods can be more relevant. Neuheimer and MacKenzie (2014) focused on the overwintering period (roughly October–April) as crucial for the development and release of groundfish eggs. As such, we calculated the Tbtm DD arrival times for two periods, January–December and October–April. The calculation is identical to that for the arrival time of SST-based DDs (Section S2.2).

The calculation of each metric requires computing a given quantity in the present and future climates. Thus simple future-climate-minus-present-climate difference maps are available for all variables (see Section S2.3). As an example, if an organism has an allometric relationship as a function of DD, then the absolute difference can imply a difference in this allometric quantity in the future. Absolute difference maps can be used to infer possible phenological changes. This usage is illustrated for bottom DD and bottom temperature differences in Section 5.

Simulations of present and future ocean-climate states are required to compute the phenology metrics. Typically an ensemble of (usually coarse resolution) ocean model simulations would be used in order to provide an estimate of the variation in these calculations. However, in this study, where our focus is on the metrics themselves and their application, we chose to use simulations from a single high-resolution ocean model, BNAM, which was designed specifically to resolve processes in the Labrador Sea and on the shelf regions of the NW Atlantic from the Newfoundland/Labrador (NL) shelf down to the Florida Strait.

BNAM output (Brickman et al., 2016; Wang et al., 2018) has been used in numerous present and future climate ecosystem and process studies of the NW Atlantic Ocean (see, e.g., Brickman et al., 2018; Greenan et al., 2019; Shackell et al., 2019; Beazley et al., 2021; Brickman et al., 2021; Wang et al., 2021b, and references therein; Czich et al., 2023). We present only basic information here. BNAM is a 1/12-degree resolution model of the North Atlantic Ocean, based on the NEMO-OPA code (Madec, 2008), on a domain that covers 7–75°N latitude and 100°W–25°E longitude. This resolution translates to a grid cell size of about 6 km in the Atlantic Canada region. Output is monthly averaged physical quantities, which are linearly interpolated in time to provide weekly resolution of the phenology metrics. The model is run in present and future climate modes. Its future climate simulations are designed to produce climatologies for four future periods and representative concentration profiles (RCPs): 2046–2065 (2055) RCPs 8.5 and 4.5; 2066–2085 (2075) RCPs 8.5 and 4.5. The production of these climatologies is accomplished by applying forcing appropriate to each period and running the model to equilibrium. The simulations produce mean states (climatologies) centered on the various time periods (including the present) but do not provide data for the intervening years. In this study, for simplicity, we used one future climate simulation, RCP8.5, 2046–2065 (2055). Hereafter the term “future climate” refers to that single simulation, which we consider to represent climatological conditions for the mid-21st century (for RCP8.5), a scenario consistent with the current rate of emissions (Schwalm et al., 2020). As BNAM’s present climatology is centered around the mid-1990s, for calculation purposes the middle of the 21st century is about 60 years in the future.

In this section we present results from the phenology change calculations. Because we used the output from a single ocean circulation model, our results should be taken as providing a general idea of what the possible changes might be. However, the high resolution of the BNAM, compared to lower resolution CMIP6 alternatives (Coupled Model Intercomparison Project phase 6; Eyring et al., 2016), allows shelf processes to be resolved and thus reveal horizontal spatial differences in the climate change metrics. The resolution of shelf processes, coupled with representations of other expected changes in the hydrological cycle incorporated into the BNAM simulation, leads to a broad range of interesting results. For each phenology metric we review the calculation, summarize the resulting projected change by the middle of the century, and discuss the result.

This calculation estimates when the present climate April SST occurs in the future, reporting this difference in weeks. The BNAM result (Figure 2) shows projected changes in the onset of spring of 3–4 weeks in most of the GoM-SS-GSL region, diminishing to 1–3 weeks toward the north and east. On the northern Labrador shelf, the onset can be negative, indicating a delayed spring. This outcome is likely due to the effects of Greenland glacial melt and other changes in freshwater input in the northern region which, when advected southward along the shelf, result in colder April SSTs in the future and thus a delayed spring.

The onset of spring calculation is one for which the BNAM result can be compared to existing present climate estimates for the region. Thomas et al. (2017), from analysis of SST data for the GoM, found an advancement of spring by approximately 0.5 days year⁻¹, similar to Friedland et al. (2023, see their figure 3). For the GoM, BNAM projects the onset of spring to be about 4 weeks earlier in 2050, translating into a rate of approximately 0.5 days year⁻¹ earlier. For the GSL, Galbraith and Larouche (2013) reported, from present climate SST data, an advancement of the 12°C isotherm, representative of early summer, of approximately 0.4 days year⁻¹, with little sensitivity to the choice of isotherm. Although more recent data indicate variability in this calculation (Galbraith et al., 2020), the BNAM projection that spring (averaged over the GSL) will be about 3 weeks earlier by mid-century translates into a similar advancement rate. Thus, the BNAM result indicates that the onset of spring will continue to advance at a rate similar to what is observed in the present climate for both the GoM and GSL.

This calculation determines when the present climate April and October SSTs occur in the future and reports this length of time as a change, in weeks, relative to the 6-month present climate period. The BNAM result (Figure 3) shows projected increases in the length of the growing season of 4–6 weeks in the GoM-SS regions, 3–5 weeks in the GSL, and 1–3 weeks toward the north and east.

The change in the length of the growing season is also a calculation for which the BNAM result can be compared to estimates based on present climate data. For the GoM, Thomas et al. (2017) reported a rate of increase of approximately 1–1.5 days year⁻¹, slightly lower than the 1.6–2.0 days year⁻¹ deducible from figure 3 of Friedland et al. (2023), while the BNAM rate is approximately 0.6 days year⁻¹ (based on a GoM average 5-week increase over 60 years). For the GSL, Galbraith and Larouche (2013) reported a rate of increase of approximately 0.9 days year⁻¹ from present climate SST, while the BNAM rate is approximately 0.5 days year⁻¹ (based on a GSL average 4-week increase over 60 years). Thus the BNAM results indicate that the lengthening of the growing season in the GoM and GSL will continue into the future but at a reduced rate compared to the present climate. The calculated BNAM length of growing season is not identical to that of Thomas et al. (2017), Friedland et al. (2023), or Galbraith and Larouche (2013), so the above comparisons should be taken with this caveat in mind.

In this calculation the annual mean present climate stratification is computed, and the time when this value occurs in both present and future climates is determined. The change in the onset of stratification is the difference in these two times, in weeks. The BNAM result (Figure 4) indicates little difference in the seasonal cycle of stratification between the present and future climates. Most of the shelf region is about 5–10 days earlier in the future, but there are areas where the future climate onset is later than present climate (negative values). Figure 4 provides an example of the possible spatial complexity of stratification changes in shelf seas when simulated by a high-resolution ocean model. The difference in density between the surface and 50 m is affected by seasonal air-sea interactions which drive surface fluxes and mixing, with associated ice formation and convection. In the eastern Canadian shelf region, estuarine circulation (driven by river inputs) and deep water intrusions also affect water column stratification (Brickman et al., 2018; Galbraith et al., 2020). The BNAM includes realistic representations of these processes and their projected changes, as well as projections of changes in glacier melt, which result in the complex spatial pattern seen in Figure 4.

This calculation determines the maximum SST in the present climate, for each grid cell, and computes the number of weeks in the future climate that are greater than this value. This metric is relevant to species currently experiencing SSTs (or Tbtms, see Section 4.5) near their thermal tolerance maximum, as it indicates the degree to which these species may be under stress in the future. The BNAM result (Figure 5) shows projected durations above the maximum present climate SST of 4–12 weeks in most of the GoM-SS-GSL-NL region. For a region on the Labrador shelf, the projected duration is shorter, as the future climate is colder than the present climate. This region is associated with the advection of Greenland glacier meltwater around the Labrador Sea.

This calculation determines the maximum Tbtm in the present climate, for each grid cell, and computes the number of weeks in the future climate that are greater than this value. The BNAM result (Figure 6) shows projected durations above the maximum present climate Tbtm of about 52 weeks for most of the GoM-SS-GSL region, with patches of lower values (<20 weeks) in shallow water nearshore and in bank regions. This result means that the bottom temperature over most of the GoM-SS-GSL region will be warmer than the present climate for the entire year. The NL shelves show a more complicated spatial pattern, with values varying from 8 weeks to 52 weeks. The negative regions (purple shading) toward the north are areas where the future climate is colder than the present climate, due to the release of Greenland glacier meltwater.

This calculation determines the time it takes for the present climate number of DD for April–September to be attained in the future, expressed as a difference (in weeks). The BNAM result (Figure 7) projects that the present climate SST-based April–September DD will arrive 1–3 weeks earlier in the future for most of Atlantic Canada and the GoM.

This calculation was performed for the standard calendar year (January–December), as well as the overwintering months (October–April), which are important to the timing of spawning (Neuheimer and Mackenzie, 2014). The BNAM result for January–December (Figure 8a) projects earlier arrival times of 5–15 weeks for the GoM/SS, 5–10 weeks for most of the GSL, and 5–50 weeks for the Grand Banks/NL shelf region (see Section S2.2 for an explanation of these high values on the NL shelf). The spatial pattern for the overwintering months (Figure 8b) is generally similar to that of the calendar year but with a lower range of values. In this case we find earlier arrival times of 6–8 weeks for the GoM/SS, about 5 weeks for most of the GSL (but note the roughly 20-week earlier region in the southern GSL), with the Grand Banks/NL shelves ranging about 0–20 weeks earlier.

An example of an absolute difference calculation for annual bottom DD is shown in Figure 9. The BNAM result projects differences of 400–700 DD for the GoM and western Scotian Shelf and less than 400 DD for the GSL and regions east and north. We use this information in Section 5.

In this section we provide examples of how the climate change metrics can be used to predict future changes in ecosystem components. Applications of phenology metrics fall into two broad categories: direct and indirect. In direct applications, formulae exist that use values of the phenology metric to predict a phenological, or other, change in an organism. Indirect applications include possible, or expected, phenological responses and non-phenology-related changes to ecosystem components. An example of the former would be a prediction of thermal stress on an organism during a phase of its seasonal cycle, from which an expected response would be to change its seasonality. An example of the latter is the use of absolute difference maps such as the absolute difference in bottom DD (Figure 9). These maps can be used to compute phenological responses, but they could also be used to predict non-phenology-related changes to ecosystem properties like habitat suitability. Although these latter applications are not phenological strictly speaking, they are included because they have broad utility in predicting possible future climate ecosystem changes.

To find applications we searched the literature for examples of species with developmental or life history traits dependent on variables contained in the climate change metrics. Examples were also taken from the work of the authors. Implicit to the use of any specific formula or species behavior based on present climate data is that it also applies to the future climate. This “all things being equal” assumption, of course, may not be strictly true, but is necessary to make predictions of ecosystem change.

Mills et al. (2013; 2017) studied the effects of marine heat waves on the GoM inshore lobster population. In this region, as spring advances, lobsters move from offshore waters into shallower coastal areas where they are available to the fishery, resulting in a characteristic rise in landings in midsummer. Mills et al. (2013) found that during the 2012 marine heat wave the onset of spring was about 21 days earlier than climatology and that this earlier onset was accompanied by a shift in the peak of the lobster landing curve also toward 21 days earlier compared to climatology (see figure 2 of Mills et al., 2013).

The BNAM change in the onset of spring calculation (Figure 2) indicates that spring would arrive about 3 weeks earlier in the GoM in 2055. Thus the lobster landing onset/peak would be expected to be 3 weeks earlier by mid-century. In a related study, Mills et al. (2017) found a relationship between the start of the lobster season and temperature anomalies at 50 m in April (see their figure 7). Based on this relationship, the BNAM 50 m April temperature change (not shown) supports a 15–20 day earlier start to the lobster season.

Pershing et al. (2015, see equation 10 in their supplementary material) derived a relationship between changes in the length of summer and GoM cod mortality: Mortality ∼ (1 + D/d); where D is half the change in the length of summer, and d is the half-length of the (summer) predation season. For the 2012 marine heat wave they estimated a 44% increase in cod mortality.

For 2055, the BNAM change in the length of the growing season for the GoM is about 5 weeks (Figure 3). Using the same value for d as Pershing et al. (2015), and assuming that all other factors remain the same, this 5-week increase in the length of the growing season means that GoM cod mortality would be about 40% higher in the future than in the present climate. The Pershing et al. (2015) relationship is the only mathematical equation we were able to find that directly relates a change in a phenology metric to an ecosystem response.

Boyce et al. (2021) investigated the decline of herring in the western SS and Bay of Fundy (WSS-BoF) region of eastern Canada in the context of the 2012 marine heat wave. Their figure S7 illustrated the effect of thermal phenology on herring population variability by showing the herring spawning season superimposed on the seasonal cycle of SST for low, average, and high (i.e., 2012) SST years, along with the temperatures above which peak cardiac function and growth rates for herring larvae have been reported to decline. Those temperatures are 1°C and 2°C above the climatological maximum SST, respectively. During the 2012 marine heat wave the SST during the herring spawning season exceeded the peak cardiac function temperature and approached the growth rate decline temperature for parts of the summer, indicating possible deleterious effects on that larval year class.

From Figure 5, the BNAM duration above the climatological maximum SST in this region is about 12 weeks, with the future climate minus present climate difference in SST maximum (not shown) roughly 1°C. Thus herring larvae in the WSS-BoF released during the 3-month spawning season would spend most of their life stage in waters above the present climate SST maximum. Because their future environment, on average, would approach but not exceed a physiological threshold, the BNAM result suggests that thermal stress is only marginally likely to occur. Thus the indirect phenological response of a shift in the onset of spawning is not expected for WSS-BoF herring by the middle of the 21st century.

Snow crab on the eastern SS and southern Gulf of St. Lawrence (sGSL) are mostly located in bottom temperatures less than 3°C (Choi and Zisserson, 2011; Zisserson and Cook, 2017; Hébert et al., 2021); thus their prime habitat is well approximated by the 3°C Tbtm isotherm. Maps of the duration that the future climate spends above the maximum present climate temperature can be used as indicators of potential habitat changes in the future.

A subdomain of the map of future climate duration above maximum present climate Tbtm (Figure 6) is replicated below (Figure 10a) with the BNAM present climate 3°C Tbtm isotherm superimposed. The BNAM isotherm corresponds well with the data sources listed above. Roughly 40% of the overall present climate isotherm area, and virtually all of the eastern SS, is predicted to be above the maximum Tbtm for the entire year in 2055. This result indicates that much of the snow crab habitat may be lost in the future. Because this analysis is based on a model simulation, the 2055 3°C isotherm can be plotted (Figure 10b). This map shows that habitat loss coincides almost perfectly with the area where all 52 weeks in the future climate are predicted to be above the maximum present climate Tbtm.

Neuheimer and MacKenzie (2014) explored the factors that affected spawning time for 21 populations of Atlantic cod (Gadus morhua) across the North Atlantic Ocean. They derived spawning time relationships using the physiologically relevant metric of bottom temperature degree days, which they referenced to the beginning of the over-wintering period (taken as October). Their figure 3 provides information on the day of year (DOY) of spawning, the number of DD, and the geographic locations of the various cod populations.

Using this information and the predictions for changes in Tbtm DD arrival time in the future (Figure 8b), and assuming all other factors remain the same, we estimated spawning times for northwest Atlantic cod populations in 2055. For example, for Georges Bank cod, BNAM predicts about an 8-week earlier arrival time for the same number of DDs, which implies a spawning time shift from about the end of February to early January. For sGSL cod (which overwinter on the eastern SS; DFO, 2000), BNAM predicts about a 9-week earlier arrival time for the same number of DDs, which implies a spawning time shift from about the middle of June to early April.

These estimates use the arrival time difference from October to April (Figure 8b). A more precise calculation would take into account the Neuheimer and MacKenzie (2014) spawning DOY for each population and compute the Tbtm DD arrival time from October to that DOY (which are roughly DOY 60 for Georges Bank cod and DOY 165 for sGSL cod). Performing that calculation (not shown) results in a reduction in the Georges Bank arrival time of about 1 week (8–7 weeks earlier) and a 1-week increase (9–10 weeks earlier) in the arrival time for the sGSL population. This result suggests that a sufficient answer is obtained from the arrival time map for a generic overwintering period (October–April).

If an organism has an allometric relationship as a function of DD, then the change in DD in the future can be translated into a projected change in this allometric attribute. For example, Neuheimer and Taggart (2007) derived a relationship between North Atlantic cod length-at-day and Tbtm DD (see their figure 5): N Atlantic cod length-at-day = 280 mm + (0.034 mm/DD) × DD. The suite of stocks in the NW Atlantic to which this formula applies include those resident on the NL shelves, GSL-ESS, and the WSS-GoM. Figure 9 shows projected increases of about 150, 350, and 600 DDs, respectively, for these three broad regions which translates into growth increments of about 0.5, 1.2, and 2.0 cm year⁻¹ in 2055. This analysis is an example of an indirect application of a climate change metric because increased growth rates in fish can affect maturation timing, prey demand, and predator-prey relations (Neuheimer and Taggart, 2007) leading to phenological changes in the ecosystem.

Bottom DD difference maps can also be used as indicators of potential future climate changes. Czich et al. (2023) developed a species distribution model for halibut in the Atlantic Canadian shelf region, including a future climate prediction of habitat change. Figure 11 is a zoom of Figure 9 DD difference map with the blue contour delineating the strong gradient in the bottom DD difference between the WSS-GoM and the ESS-GSL. This contour is also a reasonable bounding curve for highly suitable present climate halibut habitat (Czich et al., 2023, their figure 5). Thus halibut are currently distributed in a region with strong projected future climate change, suggesting that this region may become unsuitable for halibut. This possibility is confirmed by the Czich et al. (2023) prediction for the change in halibut distribution (their figure 6) which follows the pattern in DD change, with losses in the high DD change region in the WSS-GoM sharply delineated by gains in the low DD change region in ESS-GSL—effectively a west-to-east shift in favored halibut distribution. This application illustrates that maps such as Figure 9 can be used in conjunction with known species distribution information to infer potential future climate changes for a given species.

Koeller et al. (2009) examined the phenology of shrimp and phytoplankton for shrimp stocks across the North Atlantic Ocean. As part of that study they derived a relationship between the summer-to-spring egg development time and bottom temperature (their figure 1) with a slope of approximately −25 days °C⁻¹. Projected changes in Tbtm can be used to predict changes in hatching times for the four NW Atlantic stocks they considered: NL shelf + Flemish Cap, GSL, SS, and GoM (see their figure S4). First note that for a given time period, DD and mean temperature are related: $〈T〉=DD/(Ndt)$ where N is the number of discrete time periods (dt) used in the DD calculation. Applying this relationship to Figure 9 gives approximate Tbtm changes for the four regions of 0.28°C, 0.75°C, 1.0°C, and 1.8°C, respectively. These changes translate to a phenological shift to earlier egg hatching times of about 7, 19, 25, and 38 days for the four NW Atlantic shrimp stocks mentioned above.

In this article, we have presented eight phenology change metrics, based on physical variables from an ocean model, and assessed how they might change by the mid-21st century. We then used examples from the literature focused on present-day climate to illustrate the application of the phenology change metrics to future climate forecasting. We showed that the suite of phenology change metrics can provide insight into expected future impacts. On their own, there is value to using these metrics to gauge future seasonality changes, but their inclusion in climate risk or climate vulnerability assessments is also possible and may help to reduce the uncertainty of assessments, where various vulnerability metrics are usually based on changes in annual average values as opposed to changes in the timing of events.

The phenology change metrics include two that are commonly mentioned in the terrestrial and oceanographic literature (changes in the onset of spring and the length of the growing season), one unique to the ocean (changes in the onset of stratification), and five that are presented as new metrics (based on SST and Tbtm). Of the latter, two measure the time that the future climate spends above the maximum temperature in the present climate—a thermal stress indicator. Two are based on DD, measuring when a given present climate number of DDs arrives in the future—an indicator of speeding up ecosystem development times in the future. The last metric is the difference in a scalar quantity (e.g., temperature or DD) between the future and present climates, which can be applied to predict changes in habitat or allometric properties of organisms and seasonality changes. The use of DD in phenology studies is rare, with about 3% of terrestrial publications and 1% of the oceanographic literature mentioning that term (see Table S1). An example from the oceanographic literature is Asch et al. (2019), who used DD to study mismatches between phytoplankton blooms and fish spawning phenology.

In general there are no unique accepted methods to compute the (known) phenology metrics. The literature, for example, includes at least four ways to compute the onset of spring in the ocean (see Galbraith and Larouche, 2013; Friedland et al., 2015; Thomas et al., 2017; Wang et al., 2021a). We chose methods designed to accommodate our spatially large mid-latitude study region, which includes temperatures that range from below zero to above 20°C. The new phenology metrics we have presented are all adaptable to different computation methods that may be demanded by the region of interest and/or the nature of the ocean model output.

The objective of this study was to introduce new phenology change metrics relevant to ecosystem application, rather than to focus on the actual values computed (Figures 2–9). We used the BNAM high-resolution downscaled future climate ocean model for our calculations as it is able to show the possible spatial complexity of the calculations. Such complexity is manifested by regions that show counterintuitive results; for example, later onsets or shorter summers in the future (negative regions in Figures 2–6). Counterintuitive results are particularly apparent in the onset of stratification map (Figure 4), which shows both advancements and delays in onset within the GSL, likely due to changes in river inputs and deep water intrusions into the Gulf. For two metrics, the onset of spring and the length of the growing season, calculations based on data were available for the GoM and GSL. Our projections indicated continual advancement of spring, at about the same rate as the present climate data, and a lengthening of summer but at a reduced rate.

We have presented eight applications of the metrics (Section 5) to provide examples of how they could be used to predict possible species-specific responses to climate change. The applications illustrate the wide range of possible uses of model-derived climate change metrics. They also point out a dichotomy in searching for phenology applications. While some metrics may be considered to be phenological measures (e.g., changes in the length of the growing season), their application may not predict a seasonal change (e.g., GoM cod mortality; Pershing et al., 2015). Conversely, metrics not representing phenological measures (e.g., absolute difference metrics) can be used to predict changes in phenology (e.g., shrimp egg hatching times; Koeller et al., 2009).

In addition to the specific applications presented, the spatial maps could be used by managers and scientists to gain insight on how the ecosystem will change in the future. For example, arrival time maps (Figures 2, 4, 7, and 8) indicate possible earlier starts or shorter durations for events in the future and provide a sense of the potential spatial variability in these quantities. If the maps are used in this way, then using an ensemble of future climate model simulations becomes more important as that would allow the creation of mean and variability maps. In this study we have used one high-resolution future climate ocean model, BNAM, as it can resolve processes that affect bottom temperature, a key variable in the climate change metrics. As the poor performance of the low-resolution CMIP models in simulating present climate bottom temperature in shelf regions is increasingly documented (see Loder et al., 2015; Wang et al., 2024), the need to use high-resolution models becomes more apparent. Fulfilling this need will become possible as more high-resolution regional downscaled model simulations on appropriate spatial domains become available.

This study has presented climate change metrics computed from ocean model physical variables. Part of the motivation for this choice is that these metrics can be calculated from any future climate ocean model simulation. However, output from Earth system models (which include biogeochemical components) would allow access to a set of biological and chemical variables to formulate and estimate new metrics. The expansion of the set of metrics presented herein using biogeochemical variables would be a logical next step in developing useful ecosystem change metrics from future climate model simulations.

The model data and programs used to create the figures are available as separate files of Supplemental materials.

The supplemental files for this article can be found as follows:

Phenology metrics for ocean waters with application to future climate change in the Northwest Atlantic Ocean: Supplemental materials.docx.

The comments of Dr D. Boyce, Dr R. Stanley, and two external reviewers of this manuscript greatly improved its presentation and readability.

The authors have no competing interests to declare.

Contributed to conception and design: DB, NLS.

Contributed to acquisition of data: DB.

Contributed to analysis and interpretation of data: DB, NLS.

Drafted and/or revised the article: DB, NLS.

Approved the submitted version for publication: DB, NLS.

How to cite this article: Brickman, D, Shackell, NL. 2024. Phenology metrics for ocean waters with application to future climate change in the Northwest Atlantic Ocean. Elementa: Science of the Anthropocene 12(1). DOI: https://doi.org/10.1525/elementa.2024.00001

Domain Editor-in-Chief: Jody W. Deming, University of Washington, Seattle, WA, USA

Associate Editor: Maxime Geoffroy, Fisheries and Marine Institute of Memorial University of Newfoundland, St. John’s, NL, Canada

Knowledge Domain: Ocean Science