We present an extension of the previous work by Speirs et al. (2005, 2006), who modeled the population dynamics of C. finmarchicus over the entire North Atlantic. This extension took two key forms. First, we incorporated recent developments in our understanding of Calanus biology. Second, we modified the model of Speirs et al. (2006) so that it could represent the population dynamics of both C. finmarchicus and C. helgolandicus. Full mathematical details of the model, along with relevant parameters, are given in Appendix 1 (Supplementary Material). Here we will summarize the modeling framework of Speirs et al. and then the extensions to it.

The model of Speirs et al. was discrete in time and space. It covered the entire North Atlantic, ranging from 30 to 80°N and 80°W to 90°E. The population of C. finmarchicus was distributed over a regular grid of cells of size 0.5°longitude by 0.25°latitude. They had two update processes. First, the population of each cell was updated to account for development, reproduction and mortality. After these updates, the population is redistributed between cells to account for physical population transport. A separate physical model was used to create the flow-field and temperature drivers for the relevant biological and physical update. The annual cycle of food in each cell was estimated by deriving phytoplankton carbon fields from satellite sea-color observations. 1997 was used as the target year for simulations because this was the year when the Trans-Atlantic Study of Calanus (TASC) collected a large number of time series of C. finmarchicus abundance in the North Atlantic. The framework of Speirs et al. was as follows. Surface developers are made up of eggs (E), naupliar stages (N1 to N6), and copepodite stages (C1 to C5). Finally, there are diapausers (C5d) and adults (C6).

Calanus development follows the equiproportional rule, that is relative stage duration is independent of temperature (Campbell et al., 2001). Development from egg to adult can therefore be divided into a fixed number of steps, with each having identical time duration under identical environmental conditions (Gurney et al., 2001). In total, there were 57 development steps, which cover the 13 stages of Calanus development.

This framework allows the entire population to be updated simultaneously, and for the entire population to be simulated with high computational efficiency (Speirs et al., 2006). However, modeling the populations of C. finmarchicus and C. helgolandicus required one modification.

We began with the hypothesis that differences in the response of growth and development to temperature are sufficient to explain the geographic distributions of both species. In other words, all equations and parameters would be the same, except for those related to growth and development. This could not be satisfactorily achieved in the original framework. Large-scale patterns of fecundity are not only the result of the effects of environmental conditions, but also of body size. Further, the ability of animals to diapause is strongly influenced by size (Wilson et al., 2016). We therefore incorporated body size into the framework. Large-scale patterns of fecundity and diapause duration could therefore be represented as the combined effects of body size and the environment, and did not require the introduction of interspecific differences. The geographic domain used by Speirs et al. covers all regions of high C. helgolandicus abundance (Bonnet et al., 2005), and was therefore maintained.

The following biological processes are represented in our model: development, egg production, diapause and mortality. In each case, we modified the model of Speirs et al. to account for recent developments in the understanding of Calanus biology.

A recent review of the differences between the two species found that the only known relevant difference was the influence of temperature on ingestion, and thus growth, development and fecundity (Wilson et al., 2015). We therefore constrained the model by making a number of assumptions about the differences between the species based on this review. These assumptions were as follows:

Further, we take the following assumptions and simplifications about the biology and ecology of both species.

The key modeled relationships between body size, development time, egg production rate and diapause duration with temperature are shown in Figure 1.

There are no apparent interspecific differences in body size, and large-scale geographic patterns of body size are largely driven by temperature (Wilson et al., 2015). We therefore modeled body size under the simplified assumption that it is determined by temperature experienced at birth for all development classes (Figure 1A). This assumption is derived from the fact that egg size is determined by temperature (Campbell et al., 2001) and that the existence of an exo-skeleton likely greatly constrains size over all development classes. The temperature-prosome length relationship of Campbell et al. (2001) was used with a multiplier, which was fitted based on the relationship between predicted and observed female prosome length. Prosome length reduces linearly with increasing temperature. This approach contrasts with Speirs et al., which did not represent size.

Egg-adult development time was assumed to be influenced purely by temperature and food concentration. The relationship between egg-adult development time and temperature under food-saturated conditions is assumed to follow that derived by the model of Wilson et al. (2015). Development time saturates at high food levels, and we use the relationship between food concentration and development time of Campbell et al. (2001). There is a U-shaped response of development time to temperature (Figure 1B), which contrasts with the monotonically decreasing form used by Speirs et al. The computational approach is that of Gurney et al. (2001) and uses dynamic time-step constraints. This is the same approach as in Speirs et al. (2005, 2006) and it is effective in minimizing numerical diffusion (Gurney et al., 2001; Record and Pershing, 2008).

Fecundity was related to temperature, food concentration and body size. We assumed that egg production and growth are equivalent (McLaren and Leonard, 1995). Egg laying females have stopped growing and we therefore assume that carbon previously directed to growth will be used to make eggs. The growth rate equation of Wilson et al. (2015) forms the basis of our egg production rate (EPR) model for both species, with the food saturation component taken from Hirche et al. (1997). EPR therefore has a dome-shaped response to temperature (Figure 1C). Further, EPR has a saturating response to food concentration and we use a conventional allometric relationship between EPR and carbon weight, i.e., EPR ~ carbon weight^0.75. This contrasts with Speirs et al., who represented EPR as a monotonically increasing function of temperature, but using the same food response as we have assumed. We assume that 50% of adults are female.

A recent modeling study, which synthesized empirical findings, showed that maximum potential diapause duration is largely determined by prosome length and overwintering temperature (Wilson et al., 2016). We therefore modeled diapause duration using the maximum potential diapause duration equation from that study (Figure 1D). Diapause duration declines at higher temperature because of increased metabolic rates, and is shorter at smaller prosome lengths because of lower relative lipid levels and higher relative metabolic costs. We assumed that a fraction of the C5 population enters diapause at the end of the C5 stage. This fraction is dependent on growth rate, with it increasing at lower growth rates, so that more animals diapause when development conditions are poor. In the model, animals exit diapause at the end of their potential diapause duration. This differs from Speirs et al., who assumed that diapause exit was triggered by a photoperiod cue.

Mortality is modeled using a stage-dependent background rate, alongside a starvation and density dependent term. Field studies indicate that mortality in both species is stage-dependent (Eiane et al., 2002; Ohman et al., 2004; Hirst et al., 2007). These estimates of stage-dependent mortality include all sources of mortality. However, we need to distinguish between different sources of mortality to properly represent population dynamics. We therefore used a fraction of the stage-specific mortality rates calculated by Eiane et al. (2002) as the background mortality rate, with additional temperature, starvation and density dependent terms. Starvation dependent mortality was modeled in the same way for both species by assuming that it relates to growth rate; with starvation mortality only occurring below a threshold growth rate and increasing as growth rate decreases. Background mortality is temperature dependent, with mortality increasing with temperature and the relationship taking the form mortality ~ (T/8)^z. Density dependent mortality is assumed to be proportional to total biomass. Mortality was represented the same way as in Speirs et al., with the exception of starvation-dependence. Speirs et al. represented this purely as a function of food concentration. However, the differences in ingestion rate between the two species (Møller et al., 2012) show that C. helgolandicus is likely to face much greater starvation levels at temperatures below approximately 11°C. We therefore viewed growth rate as a better indicator of starvation than food concentration.

Seasonal cycles in food concentration, temperature and oceanic circulation drive the model. The only data with sufficient spatial and temporal coverage of food concentration are satellite estimates of sea surface color. SeaWIFS satellite estimates of chlorophyll were therefore used to derive food fields.

Insufficient observations are available for 1997. We therefore used a climatological 8 day mean of chlorophyll concentration from 1998 to 2000. There is a poor relationship between time series derived from SeaWIFS and field estimates of chlorophyll (Speirs et al., 2005; Clarke et al., 2006). We used the estimates of Clarke et al. (2006), who developed a statistical methodology, where thin plate regression splines modeled local estimates of chlorophyll concentration in relation to SeaWIFS estimates, bathymetry and time of year. Field estimates of chlorophyll concentration in the top 5 m were used, assuming they reflect chlorophyll concentration throughout the vertical distribution of Calanus. However, it is possible that this does not fully capture deep-water chlorophyll concentrations. Phytoplankton abundance was calculated assuming that 1 mg m⁻³ of Chl a is equivalent to 40 mg Cm⁻³ (the approximate median of the values reported by Parsons et al. (1984). Estimates of food extend to regions covered by sea ice, where we masked food levels to zero. This mask was derived from 1997 satellite percentage ice cover from the Defense Meteorological Satellite Program's (DMSP) spatial sensor microwave/imager (SSM/I) (Comiso, 1997). The approach taken to food was the same as in Speirs et al.

Temperature and velocity fields come from the Nucleus for European Modeling of the Ocean (NEMO) Ocean General Circulation Model (OCGM) (version 3.2) (Madec, 2012). The forcings and model implementation are described in Yool et al. (2011). NEMO is resolved at 64 vertical levels, and it resolves the primitive equations on a C-type Arawkawa grid. Ocean surface forcing comes from the DFS4.1 fields produced by the European DRAKKAR collaboration. This differs from Speirs et al., who used the OCCAM model to derive temperatures and flow fields. Computation of the NEMO model was performed using the free Java tool Ichthyop version 3.2 (Lett et al., 2008).

We assumed that surface developers experience the temperatures and velocities which occur at a depth of 20 m. Diapause depth varies in space. We therefore derived a map of diapause from the data reported by Heath et al. (2004). A loess smooth was used to estimate the median diapause depth in regions close to where Heath et al. (2004) reported data. Where the smoothed estimate exceeded bathymetry, we used a depth 10 m shallower than the bathymetry at a location. In other regions we assumed that if bathymetry was greater than 800 m that diapause depth was 800 m. For locations where bathymetry was shallower than 800 m we used the predictions of a general additive model which related median diapause depth with bathymetry using the data of Heath et al. (2004). Transport updates occurred every 7 days. At the start of each time step, 100 seeds were placed at the center of each model cell. Particle trajectories over a 7-day period were then calculated, and transition matrices were calculated to show the proportion of particles which move to each nearby cell. The approach outlined above was in agreement with Speirs et al.

Our underlying goal was to reproduce the biogeography of both species displayed by the CPR. We therefore carried out an extensive set of simulations to assess how well different parameter sets could reproduce the geographic distributions of both species.

As discussed in Section 2.2, laboratory and field data were used to derive the following traits: development time, growth, fecundity, diapause duration, background mortality and body size. The remaining free, i.e., unknown, parameters related to the equations for diapause entry and starvation and biomass dependent mortality. We initially sought a single parameter set for mortality and diapause entry that would result in credible predictions of geographic distributions for both species. However, a large number of exploratory runs showed that this was not possible. We therefore sought parameter sets that reproduce the geographic distributions of both species while minimizing the differences between the model parameters of both species. A suite of runs showed that this was only achievable by assuming that mortality responded differently to temperature in both species.

Model parameters were derived by simultaneously altering the terms for mortality and diapause entry for both species and recording each parameterization's fit to CPR abundance data. First, CPR data was split into cells of dimension 2°E and 1°N, and we then removed cells without a CPR abundance record for each month of the year. Annual mean abundance was then calculated by averaging the mean abundance of the mean monthly abundance for C5 and adults in each cell.

This resulted in 333 cells for model comparisons. Each CPR abundance record represents approximately 3 m³ of filtered seawater (Richardson et al., 2006). Therefore, CPR data must be divided by 3 to get estimates of abundance per m³. This must then be multiplied by a further conversion factor of 20 (Speirs et al., 2006) to provide estimates of abundance (m⁻²) over the top 100 m of the water column.

Simulations began by seeding a large number of eggs over the entire North Atlantic and in the eastern North Atlantic for C. finmarchicus and C. helgolandicus, respectively. The model was then run to a quasi-stable state and we then calculated the correlation coefficient (r) between predicted annual surface abundance (m⁻²) and CPR abundance (m⁻²).

We report two sensitivity experiments. First, we show the geographic distributions of both species when there are no interspecific differences in free parameters, i.e., only differences in growth, development and fecundity are assumed. In this case we are using the diapause entry and starvation and temperature dependent mortality parameters for C. helgolandicus for both species.

Our initial model of diapause duration used a model of maximum potential diapause duration (Wilson et al., 2016), which possibly results in diapause durations which are unrealistically long. We therefore carried out a sensitivity analysis which relates the ability to reproduce the geographic distributions of both species to the assumptions for diapause duration and temperature dependent mortality. Temperature dependent mortality is proportional to (T/8)^z for temperature T (°C). The parameterization assumed different values of z for each species.

Figures 2, 3 compare the model predictions and CPR estimates of bimonthly abundance for C. finmarchicus and C. helgolandicus respectively. Table 1 shows the correlation coefficients between monthly modeled and CPR abundance for both species. The large-scale geographic pattern of C. finmarchicus abundance was successfully reproduced in comparison with CPR. The correlation coefficient between simulated mean annual abundance and CPR abundance over the 2°E by 1°N cells is 0.75. Bimonthly comparisons between C. finmarchicus predictions and the CPR abundance are shown in Figure 2. Importantly, we reproduced the relatively high abundance of C. finmarchicus in the West Atlantic in autumn. In addition, the model predicts a year round surface population in coastal waters in the West Atlantic, in accordance with CPR. However, it perhaps over-predicted abundance in November and December.

A comparison of bimonthly predictions of C. helgolandicus abundance with the CPR abundance is shown in Figure 3. The correlation coefficient between predicted mean annual abundance and CPR abundance over the 2°E by 1°N cells was 0.76. Importantly, C. helgolandicus was restricted to the continental shelf. The autumn bloom of C. helgolandicus in the North Sea was also reproduced. However, predicted abundance in November and December in the region to the south west of the British Isles appears too high.

Figure 4 shows simulated combined abundance for stage C5 and adult C. finmarchicus compared with those from the time series. Predicted peak abundances are within a factor of 2 of those recorded in the time series, with the exception of the Westmann Islands. OWS I is notable for getting the scale of the first generation very accurate, but we predicted a much larger second generation than is apparent in the time series. We failed to show the apparent sharp increase in C5 and adult at OWS M before day 100. Additionally, the second peak in C5 and adult abundance at OWS M appears to be time shifted by approximately 50 d.

We compare predictions for C. helgolandicus with field time series and time series derived from CPR in Figure 5. The timing of the autumn peak of C. helgolandicus abundance at Stonehaven was successfully reproduced. However, we failed to reproduce the small spring bloom. Predicted time and the magnitude of peak abundance was close to that in the L4 time series. However, abundance appeared to be over-predicted during winter.

Predicted EPR is compared with field time series at OWS M and L4 for C. finmarchicus and C. helgolandicus respectively in Figure 6. Predicted C. helgolandicus EPR is lower in the first half of the year of the time series, and is slightly time shifted compared with the time series. Predictions depart significantly from the times series in the second half of the year, with EPR being significantly higher than in the time series. The C. finmarchicus EPR time series at OWS M is of short duration. We can therefore only make a limited comparison. However, the predicted EPR is approximately the same as the median EPR in the time series.

In the results shown in Section 3.1, the only differences between the species are the relationship between growth, development and fecundity and temperature, and a parameterized difference in the response of mortality to temperature. Figure 7 shows the predicted geographic distribution of C. finmarchicus when the temperature-dependent mortality parameter for C. helgolandicus was used. The geographic distribution in the west Atlantic is successfully reproduced. However, the geographic distribution in the east Atlantic is too southerly, with a large population predicted to exist in the Celtic Sea.

Exploratory simulations showed that the C. helgolandicus predictions were sensitive to diapause assumptions. First, the model performed well if C. helgolandicus was assumed to remain at the surface year round and to never diapause. In fact, this simplified model arguably performed better than the original. The key features of the distribution of C. helgolandicus were largely reproduced, with the correlation coefficient (0.78) of model performance compared with CPR actually improving in comparison with our original model.

Further exploratory simulations showed that the state of populations of C. helgolandicus is sensitive to diapause duration. A sensitivity analysis showed that small changes to diapause or mortality assumptions can result in C. helgolandicus becoming an oceanic species. Figure 8 shows the correlation coefficient between predictions and CPR abundance of C. helgolandicus under varying assumptions for diapause duration and the scaling of mortality with temperature. A small reduction in how steeply mortality scales with temperature results in a reduction in model performance, with C. helgolandicus becoming an oceanic species. Likewise, an increase in diapause duration can result in C. helgolandicus becoming an oceanic species. Notably, the high sensitivity to changes in temperature dependent mortality was not evident diapause duration is reduced by 60%, which is potentially a more biologically realistic assumption for diapause duration.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.