49 Spatial implementation of multi-stanza and IBM

Ecospace biomass dynamics

For functional groups not represented by multi-stanza population dynamics accounting, Ecospace represents biomass (B) dynamics over a set of spatial cells (k) with the spatially discretized rate formulation

[latex]dB_{ik}/dt=g_iQ_{ik}-Z_{ik}B_{ik}-(\sum\limits_km_{ikk'})B_{ik}+\sum\limits_{k'}m_{ik'k}B_{ik'} \tag{1}[/latex]

where Bik is the biomass of functional group i in spatial cell k; g_i is conversion efficiency of food intake by group i into net production; Q_ik is total food consumption rate by group i in spatial cell k; Z_ik is total mortality rate of group i biomass due to predation, fishing, etc.; m_ikk΄ is instantaneous movement rate of group i biomass from cell k to cell k΄; and m_ik΄k is movement rate of group i biomass from cell k΄ to cell k.

All of the terms on the right hand side of Eq. 1, except g_i, are treated as dynamically variable over time so as to reflect changes in food availability (Q_ik), fishing effort and predation risk (Z_ik), and seasonal changes in movement patterns (m_ikk΄). Food consumption rates Q_ik are calculated as sums over prey types j (i.e., Q_ik = ΣjQ_jik). Likewise, total mortality rates are calculated as sums over predator types and fishing fleets f: Z_ik = M_oi + Σ_f F_ifk + Σ_jQ_ijk/B_ik , where M_oi is unexplained mortality rate, the fishing rate components F_ifk by fleets f are predicted from spatial distributions of fishing effort for each “fleet” f over the grid cells k, and the Q_ijk/B_ik ratios represent predation rate components of M (i.e., M_ijk = Q_ijk/B_ik) calculated from predator j consumption rates Q_ijk.

The Ecospace grid cells are arranged as a rectangular grid with rows r and columns c, so that each cell k exchanges biomass directly only with those cells k΄ that are in adjacent rows and columns. If cell k represents row r, column c, then k΄ is restricted to cells (r – 1,c), (r + 1,c), (c – 1,r), and (c + 1,r). Exchanges at the map perimeter are set to zero, except for groups that are assumed to be advected across the map, in which case biomasses at the map boundary are set to constant (Ecopath base estimate) values.

Original representation of multi-stanzas in Ecospace

When the multi-stanza option was originally developed for Ecosim, it was not incorporated directly into Ecospace. Instead, each stanza was treated as its own higher-order functional group for Ecospace biomass-dynamics calculations without accounting for age structure within the stanza. Rather, the age-structure of each stanza was assumed to be in equilibrium. Body weight was computed grid-wide (not cell-specifically) for each stanza. Feeding rates were assumed proportional to a relative search-area index P_s calculated from a prediction of the numerical abundance of the stanza N_s as P_s = N_s̅P_s, where ̅P_s is the initial (t = 0) per-capita mean of the relative area-searched index P_s,t, i.e.,

[latex]\bar P_s=\sum\limits_aN_{a,0}W_{a,0}/\sum\limits_aN_{a,0}\tag{5}[/latex]

Dynamics of the numbers in each stanza N_s were computed for each cell by the differential equation

[latex]dN_{sk}/dt=R_{sk}-Z_{sk}N_{sk}-(\sum\limits_{k'}m_{ik'k}N_{sk})\tag{6}[/latex]

where R_sk is an approximate difference between recruitment (incoming) rates and exit (to next stanza) rates for stanza s in spatial cell k. If the age structure within the stanza is assumed to remain near equilibrium, the R_sk term in Eq. 6 can be approximated as

[latex]R_{1k} = E_{tk} (1-exp(Z_{sk}(a_2(1))/12)) \; \text{for} \; s = 1 \tag{7a}[/latex]

[latex]R_{sk} = N_{s-1,k} Z_{s-1,k} / (\exp(Z_{s-1,k}(a_2(s-1)-a_1(s-1))/12)-1) \; \text{for} \; s>1 \tag{7b}[/latex]

Eq. 7a represents egg production rate minus survival rate to the age at exit from stanza s = 1. Egg production is assumed to be approximately proportional to biomass B_sk of the oldest (adult) stanza s in cell k. Eq. 7b is derived from the equilibrium of the delay-differential equation for N_s that results from assuming spatial gain and loss rates to be approximately balanced, so that the dominant effects on N_s are gains from individuals progressing from the previous stanza and from losses of individuals as they progress to the next stanza and mortality within the stanza.

The equilibrium assumption needed for derivation of Eq. 7 can lead to inaccurate predictions because it can result in incorrect size distributions if incoming and outgoing numbers are not in balance, and size then affects the predation-rate parameters (areas searched, maximum prey-consumption rates). Eq. 7 is a relatively poor approximation for both egg production and net rates of numbers gained through graduation from younger stanzas and loss to older stanzas, so this early version of Ecospace tended to predict incorrect absolute values for cell-specific numbers N_sk relative to Ecospace-predicted cell-specific biomasses B_sk, but the predicted spatial distributions of abundances were at least qualitatively reasonable. In past applications, Ecospace generally predicted that N_sk was relatively high in cells with high egg production, in cells with favourable habitat, and in cells near seasonally varying optimum migration positions for migratory stanzas.

Predicting multi-stanza spatial distribution from continuous mixing-rate models

In an effort to avoid the large computer-memory requirements and massive accounting calculations (for typical models, on the order of 10³ more calculations per time step) required for replicating the full age-structure accounting for multi-stanza groups for every grid cell of large Ecospace models, a simple approach based on combining the overall multi-stanza population accounting (as described in the Ecosim Age-structured dynamics chapter, eq. 2-4 above) with the relatively simple Eq. 7 diffusion model for predicting relative spatial abundances by stanza. This approach depends on two key assumptions:

1. the Ecosim multi-stanza accounting system can be applied for each multi-stanza population as a whole (summed over all Ecospace grid cells), given reasonable estimates of mean food consumption rates q_s,t and mortality rates Z_s,t averaged over the grid cells (a basic assumption that is made anyway in the non-spatial Ecosim representation of any large area) and
2. the diffusion model, Eqs. 6-7, gives reasonable predictions of the relative distribution of the biomass of each stanza over grid cells whether or not the absolute numbers N_sk are predicted correctly, hence preserving effects of complex spatial-overlap patterns among stanzas.

We then perform the Ecospace time solution on monthly time steps using the following four-step procedure:

1. We use the results from integration of Eqs. 6-7 to apportion the spatial distribution of total stanza biomass B_st  over spatial cells k to give B_sk relative cell biomasses using B_sk = B_stN_sk/Σ_kN_sk.
2. The spatial B_sk biomasses (and relative predator-search areas P_sk = P_stN_sk/Σ_kN_sk) are then used in the Ecosim foraging arena and fishing rate calculations for each cell k to predict food-consumption rates Q_sk and mortality rates Z_sk.
3. Biomass-weighted average food-consumption rates ̅q_s,t, and mortality rates ̅Z_s,t, for the whole population are calculated as,[latex]\bar q_{s,t}=\sum\limits_k B_{sk}q_{sk}/B_{st} \ \ \text{and} \ \ \bar Z_{s,t}= \sum\limits_k B_{sk}Z{sk}/B_{st} \tag{4a,b}[/latex]
4. The system-scale multi-stanza accounting is done by means of Eqs. 2-4 with the biomass-weighted averages (̅q_s,t, Z_s,t) to give predicted total population age and size structure and total stanza biomasses B_s,t+1 at the start of the next month.

This procedure retains some information about predicted changes in spatial abundance patterns due to mixing processes and spatial variation in mortality rates Z_sk because Z_sk is included in the prediction of relative numbers N_sk by cell from Eqs. 6-7, but it discards information about spatial variation in growth rates qso in favour of using a single system-scale prediction of body growth (Eq. 3 with consumption rate q_s,t represented by ̅q_s,t).

Further, it fails to account for the cumulative divergence that can take place in both age and size structure for relatively sedentary species resident in spatial cells that are protected from fishing. That is, for “adult” stanzas containing many age classes, it fails to represent the potential accumulation of older, more fecund animals in protected areas, considered by some to be a key benefit of marine protected areas. One possible solution to allowing accumulation of large adults in specific areas is to split the oldest stanza group into a number of stanzas, but doing so is an approximate fix rather than a solution.

For resident species, the mixing-model approach also fails to account for regional variation in growth rates associated with spatial cells that have higher basic (primary and lower-trophic-level) productivity or reduced intraspecific competition due to limited recruitment. For these reasons, the mixing-model approach is best suited to analyses of pelagic systems, where relatively high mobility results in averaging of feeding and mortality rates over substantial areas. When used for systems with many resident or sedentary species, the approach is potentially misleading and should be used only to provide computationally “quick-and-dirty” policy screening for options such as size and spacing of MPAs, to be followed by more careful screening according to the more detailed individual-based approach described below.

Individual-based approach for predicting spatial patterns in growth, survival, and distribution

Most regional populations exhibit at least some degree of localized or cell-scale variation in recruitment, body growth, and survival rates, and erosion of this local structure has serious implications for maintenance of both biodiversity and overall productivity. The original and mixing-model approaches described above cannot adequately capture such local structure, which can result from the cumulative effects of the development of a fishery or from MPAs. Ecospace therefore includes a much more detailed and realistic approach to the representation of localized trophic-interaction effects based on concepts of individual-based modeling (IBM).

In the IBM approach, we retain the spatial biomass-dynamics accounting for non-multi-stanza species represented by Eq. 1 and the multi-stanza population dynamics accounting of Eqs. 2–4, but rather than solving Eqs. 2–4 once for each stanza using spatially averaged (grid-wide) food-consumption and mortality rates, we divide the age-0 recruits for each multi-stanza population (N_0,t) into a large number np of packets (cohorts). Each packet is assumed to represent some number of identical individuals of the population, and all packets from the monthly recruitments start out with the same individual biomass and numbers at recruitment (N_p,0,t = N_0,t/n_p). Each packet is then followed independently as it moves among spatial cells on the grid. This approach is similar to that recommended by Rose et al.^[2] and Scheffer et al.^[3]

The growth-survival Eqs. 2-3 are then solved for each packet, yielding its predicted age and size dynamics (N_p,a,t and W_p,a,t). Packets are discarded from the overall population when they reach a maximum age (denoted amax) beyond which N_p,a,t is negligible. Each packet p is assigned an initial spatial position X_p,0,t,Y_p,0,t, and movements of the packet over time are predicted from both random (diffusive) and oriented (migratory) changes in position. At each simulation time step, the ecological conditions (food intake rates, mortality rates) for the spatial cell in which each packet is located are used in Eqs. 2–3. The overall accounting for B_s,k and P_sk needed for trophic-interaction predictions (impacts from and on biomasses of non-stanza species in each cell k) then involves simply summing B_p,k,t and P_p,k,t packet biomasses and predation search areas over those packets present in each cell k, before foraging arena predictions of Q_sk, Z_sk are performed for that cell.

The obvious advantage of the IBM approach is that it retains the cumulative history of each packet’s space-use pattern, in the form of the packet’s numerical (worth) and body size (weight) states. For sedentary species, local differentiation in growth and accumulation of older animals is represented by how packets in different local areas (cells) fare over time. Further, through use of restricted movement rules, collections of packets can easily be made to form distinctive local populations, presumably key units of local adaptation and biodiversity, including as a consequence of local environmental conditions.

A disadvantage of the approach is that it requires massive computation, both for the survival-growth calculations and for movement of a sufficient number of packets over the simulated grid to permit realistic spatial distributions and variation. This number must be determined by trial and error; the number of packets must be increased until results stop changing. Most of the computational effort (typically about 90%) ends up being in the simulation of movement as changes in the locations X_p,a,t, Y_p,a,t.

Monthly survival-bioenergetics updates for each packet are based on food intake and mortality rates predicted for the spatial cell where the packet is located at the start of the month. No attempt is made to integrate q or Z rates over times within the month spent in different cells; doing so would be prohibitively computationally intensive. This omission amounts to assuming either that cell sizes are set large enough that most movements over any month occur within a single cell or that spatial correlation in productivity and predation risk among nearby cells there is reasonably high, so movements over such cells would result in the same predicted food intake and mortality rates obtained from the initial cell. Effects of violating this assumption could be tested if the model were run with varying grid cell sizes.

The initial or spawning position for each packet (X_p,0,t, Y_p,0,t) is set to the centre of a cell k, where the probability of recruiting to cell k is set equal to E_kt/Σ_kE_kt and E_kt is the predicted total egg production in cell k for month t summed over all packets that are in cell k at the start of the month. This procedure allows spawning to occur well away from locations of larval settlement or juvenile growth because larval dispersal and juvenile migration can be explicitly represented, through either different or similar movement-simulation rules as used for packets of older fish. In particular, the IBM approach “encourages” formation of local stock structure; recruitment tends to occur near centres of egg production. In the context of MPAs, lower mortality rates Z in designated cells can result in the accumulation of older, more fecund fish, and those cells can thus become local areas of high reproduction.

Monthly movements by each packet are simulated as a set of ns increments to the X,Y values that determine location on the grid. The user specifies an average annual movement distance, which implies an average monthly movement distance. The number of moves ns each month is then set so that the distance per increment cannot exceed the width of one cell. Each movement is made only in a cardinal direction (N,S,E,W), so that only X or Y (not both) changes for each move. The probability of choosing each of the four directions, k΄, is set to m_skk΄/Σ_k΄m_skk΄, where m_skk΄ is the instantaneous movement rate from cell k to k΄ calculated for the continuous biomass model (see Eq. 1). As noted above, the m_skk΄ can be set equal for all k΄, to represent purely diffusive movement, or biased to represent avoidance of cells with unsuitable habitat, movement toward preferred habitats, or seasonal migration patterns. This method for choosing movement directions makes it possible to use the same user interface for entering assumptions about movement distances and orientation for multi-stanza populations as for groups represented only by biomasses, and it ensures that the multi-stanza movement patterns are broadly comparable with predictions from the computationally faster continuous mixing-model version of Ecospace.