ABSTRACT

Individual-Based Modeling (IBM) techniques offer many advantages for spatially explicit modeling of marine fish early life history. However, computationally efficient methods are needed for incorporating spatially explicit circulation and prey dynamics into IBM's. Models of Nutrient-Phytoplankton-Zooplankton (NPZ) dynamics have traditionally been formulated in an Eulerian (fixed spatial grid) framework, as opposed to the pseudo-Lagrangian (individual-following) framework of some IBMs. We describe our recent linkage of three models for the western Gulf of Alaska: 1) a three-dimensional, eddy-resolving, wind- and runoff-driven circulation model, 2) a probabilistic IBM of growth and mortality for egg and larval stages of walleye pollock (Theragra chalcogramma), and 3) an Eulerian, stage-structured NPZ model which includes production of larval pollock prey items. Individual fish in the IBM are tracked through space using daily velocity fields generated from the hydrodynamic model, along with self-directed vertical migrations of pollock appropriate to each life stage. The NPZ dynamics are driven by the same velocity, temperature and salinity fields as the pollock IBM, and provide spatially and temporally varying prey fields to that model. The resulting prey fields yield greater variance of individual attributes (e.g. length), relative to models with spatially uniform prey. Practical issues addressed include the proper time filtering and storage of circulation model output for subsequent use by biological models, and use of different spatial grids for physical and biological dynamics. We demonstrate the feasibility and computational costs of our coupled approach using specific examples from the western Gulf of Alaska.
 
 

1. Introduction

Modern fisheries research continues to establish important linkages among fish populations, their prey and predators, and characteristics of the immediate physical environment (e.g. circulation, temperature). As interactions among populations and the environment are clarified by field measurement and laboratory experiment, there is a growing need for fisheries models that simultaneously include several trophic levels and specific populations. A single modeling approach may not serve well for all trophic levels and populations, however.

At present the majority of marine ecosystem models are formulated in terms of the aggregate properties of each model component, e.g. the mean biomass of phytoplankton per unit volume or the mean number of fish per unit area. In effect, such models follow the evolution of the "mean individual" of represented populations through time. An increasingly popular class of population models, commonly referred to as "individual-based" models, or IBMs, keep track of a number of distinct individuals within a population, each of which interacts with other individuals and its physical environment based on its present state and possibly its past history. Such models are frequently stochastic in design. Recent examples of this approach include DeAngelis et al. (1991), Rose et al. (1993), Crowder et al. (1993), Werner et al. (1993), Werner et al. (1996), and Hinckley et al. (1996). Continuing developments in computer speed and architecture have rendered the individual-based approach to fisheries modeling increasingly tractable.

Realistic population or ecosystem modeling may require an explicit treatment of spatial history. Many marine species spend their different life stages in substantially different physical environments; for example, eggs may be laid in a near-bottom, near-coastal environment, while larvae develop near the surface in shallow continental shelf waters and juvenile and adult years are spent in the deeper ocean. Even for marine species with no appreciable horizontal transport, systematic vertical migration through light, temperature, and prey concentration gradients can play a significant role in determining growth and survival.

The present work considers the simultaneous inclusion of circulation and prey fields in a spatially explicit, three-dimensional IBM for walleye pollock (Theragra chalcogramma) near Shelikof Strait, Alaska. An overview of the area and the life history of pollock spawning in Shelikof Strait is shown in Fig. 1. The early life history of this population has been studied extensively by the Fisheries Oceanography Coordinated Investigations (FOCI) program; recent summaries of this research are presented in Schumacher and Kendall (1995) and Kendall et al. (1996) . Spawned eggs develop into larvae and juveniles as they are advected to the southwest with the prevailing currents. The models we use to capture this early life history include passive advection and active vertical locomotion of individual fish through their habitat, and feeding on prey items whose concentration varies in space and time.

Specifically, we describe our experience with coupling three complex models: 1) The three-dimensional, eddy resolving, semi-spectral primitive equation circulation model (SPEM) of Haidvogel et al. (1991), adapted to this region by Hermann and Stabeno (1996). This physical model is driven by winds and runoff appropriate to the region, and has been calibrated using current meter and drogued drifter data (Stabeno and Hermann, 1996). 2) The stochastic, spatially explicit IBM of walleye pollock described by Hinckley et al. (1996). Output from an early version of this model has been compared with larval surveys spanning several years (Hermann et al., 1996). The IBM now includes prey selection based on stage and size of pollock, and the effects of turbulence on feeding, as described in Megrey and Hinckley (1998, this volume), and summarized in Fig. 2. 3) The three-dimensional Nutrient-Phytoplankton-Zooplankton (NPZ) model of Hinckley et al. (In prep). The NPZ model includes zooplankton segregated into two species groups: Neocalanus spp., the biomass dominant of the region, and Pseudocalanus spp., which are a major prey item for the larval pollock (Napp et al., 1996; Incze et al., 1997; Kendall and Nakatani, 1992). Pseudocalanus are further segregated into 13 developmental stages, spanning eggs to adult copepodites (Fig. 3). The NPZ model is formulated on a fixed three dimensional grid. Both of the biological models receive input from the circulation model in the form of currents, temperature and salinity (Fig. 4). The IBM also uses wind data directly, to calculate effective turbulence levels at larval depths through the empirical formula of MacKenzie and Leggett (1993).

In section 2, as justification for the overall modeling approach, we present an overview of two basic categories of spatially explicit models, Lagrangian and Eulerian, and their relative merits. In section 3, we consider applied and theoretical issues relevant to coupling the Lagrangian biological (IBM), Eulerian physical (SPEM), and Eulerian biological (NPZ) models, and describe our present coupling techniques. Relevant coupling issues include: the use of time-filtered versus unfiltered circulation fields for tracking individuals, the choice of a time step for updating individual positions, the use of different spatial grids for biological and physical models, and the choice of spatial boundary conditions for each model. Issues relevant to coupling the IBM with SPEM receive the most extensive treatment, as we have the longest experience with these two models. In Section 4, we describe results of float tracking experiments using SPEM output, which bear directly on coupling issues. We then briefly describe results from the NPZ model driven by SPEM velocities, and their effect when used as prey fields for the IBM. Finally, we summarize computational and storage requirements for the models as currently implemented. Section 5 presents the conclusions.
 

2. Eulerian and Lagrangian Models

Spatially explicit biological models may be constructed using either Eulerian and Lagrangian frames of reference, a distinction drawn by other plankton modelers (e.g. Lande and Lewis, 1989). An Eulerian biophysical model is here defined as one which follows the evolution of some quantity at discrete, fixed physical locations. Typically these locations are the fixed grid points used by a numerical model, although analytical solutions throughout the domain are sometimes possible for simple Eulerian models. Changes in any quantity are due not only to local biological processes, but also to advective and diffusive exchange with adjacent locations, e.g.:

    (1)
 

Here B represents the modeled quantity, u = (u,v,w) is the fluid velocity, k represents horizontal and vertical eddy diffusion coefficients, x = (x,y,z) is the spatial location, t is time, and f_E represents changes due to biological processes (birth, growth , death, etc.) at a fixed location. This can expressed verbally as
 

local rate of change = advective change + diffusive change + nonconservative biological change
 

Note that f_E is a function of space and time, and may be either deterministic or probabilistic. This Eulerian biological model, with f_E deterministic, is the form typically used by the oceanographic community for NPZ dynamics (e.g. Jamart et al., 1977; Jamart et al., 1979; Frost, 1987; Frost, 1993; Wroblewski, 1977; Lewis et al., 1994; Moisan and Hofmann, 1996a; Moisan et al., 1996; Franks and Walstad, 1997).

In the hydrodynamic literature a "Lagrangian reference frame" is one which moves with a discrete parcel of fluid. By extension, we may define a Lagrangian biological model as one which follows an individual organism or group of organisms moving through space as it is advected by the ambient currents. Mathematically we represent changes following an individual (or group) i as:
 

     (2)

where

     (3)

Here B_i (t) represents some property of the individual (or group) at time t, x_i(t)= (x_i[t],y_i[t],z_i[t]) represents the spatial path through time, and f_L represents any nonconservative changes due to biological processes. Such changes may include processes based on past history. The quantity u_i represents each component of the fluid velocity u plus any directed motions, such as vertical migration due to buoyancy or locomotion. Note the potentially "individual-based" nature of the Lagrangian model, though in fact such models need not be IBM's per se. Examples of Lagrangian models which are not individual-based include those of Flierl and Davis (1993), Botsford et al. (1994), and Moisan and Hofmann (1996b).

When a conservative (or linearly decaying) substance is being tracked (f_L =f_E = 0), the density of "individuals" represents concentration of the substance. In that case the Lagrangian and Eulerian approaches are in fact mathematically equivalent, and both have been used in pollutant dispersal studies. The Lagrangian approach is often found to yield better resolution of strong gradients (fronts and streaks) of concentration, due to spurious diffusion introduced by numerical techniques used to solve the Eulerian equations. Conversely, the Eulerian approach is more useful when broad spatial coverage is required (i.e. relatively uniform concentration over a wide area). Some pollutant dispersal models use the Lagrangian approach in the initial stages of dispersion from a point source, then extrapolate particle densities onto an Eulerian grid to model the later stages of dispersion (Dimou and Adams, 1992).

A Lagrangian approach (e.g. an IBM which tracks positions) is especially preferable to the Eulerian approach when past history influences present response to the environment. However, as in pollutant dispersal models, Lagrangian IBM's may require huge numbers of individuals to ensure adequate spatial coverage of the modeled population in its domain. The problem is especially acute when mortality is high and dispersion is rapid. In these situations, some way must be found of reseeding the population to make up for lost individuals; even with such reseeding there may be a steady loss of fine-scale information about the concentration field of remaining individuals. Rose et al. (1993) and Scheffer et al. (1995) have explored ways of reseeding lost individuals in spatially aggregated, high-mortality situations through random sampling of the current surviving population. Perhaps this approach could be applied to the spatially explicit IBM's as well, by treating location as another attribute of individuals. Another approach would be to extrapolate the IBM results onto an Eulerian grid once the population becomes widely dispersed. These issues merit future detailed exploration, but are beyond the scope of the present paper.

In practice, typical oceanic populations of plankton are too large, and too widely dispersed, for an individual-based approach to be feasible in three dimensions. One- and two-dimensional implementations have been successful, however, for both phytoplankton (Woods and Onken, 1982; Woods and Barkmann, 1993; Platt and Gallegos, 1980; Falkowski and Wirick, 1981) and zooplankton (Batchelder and Williams, 1995). Here, we choose to use an IBM for young pollock, an Eulerian NPZ model to generate prey fields, and an Eulerian physical model (SPEM) to generate velocity, salinity and temperature fields for both.
 

3. Methods for coupling Lagrangian and Eulerian models

3.a. IBM-SPEM coupling

At its core, our "Lagrangian" problem is one of accurately tracking an individual in a time- and space-variable, partially unknown environment. When a hindcast of a population for a specific time period is attempted, we must consider both the accuracy of the hindcast velocity field and the accuracy of the scheme used for tracking individuals within that velocity field. Both the physical and biological environment experienced by the tracked individuals - in particular, the prey field they experience - may be strongly affected by these factors.

The Navier-Stokes equations which form the basis of hydrodynamic models are an exact representation of the fluid evolution, but can never be solved exactly (that is, at arbitrarily high spatial and temporal resolution) on any finite computer. Nonetheless, hydrodynamic models have achieved some success in resolving both large scale and, to a lesser extent, mesoscale (~10km) velocity features of coastal and open-ocean environments. Increasingly sophisticated data assimilation schemes have allowed for ever more accurate hindcasts and forecasts of velocity fields using these equations when hydrographic and velocity data are available (e.g. Fukumori and Malanotte-Rizzoli, 1995). Even when such data are not available, or when the model is not intended as a hindcast, one may wish to reproduce the dynamics of a particular region in some statistical sense; for example, to achieve the correct scale and frequency of passage of eddies and correct mean flows, as in models of the Gulf Stream.

Ideally we wish to use a three-dimensional hydrodynamic model to generate accurate velocity fields through which individuals can be tracked in a stochastic IBM with arbitrarily high spatial and temporal resolution. Several possibilities exist for coupling a circulation model with a stochastic IBM:

i) Run the hydrodynamic and individual-based models in parallel, updating individual positions at each new time step of the hydrodynamic model. Run the combined model many times, one run for each realization of the IBM.

ii) Run the hydrodynamic and individual-based models in parallel just once, computing all realizations of the IBM simultaneously.

iii) Run the hydrodynamic model once, tracking a number of passive floats whose initial position corresponds to initial organism concentrations in the model domain. Store the resulting Lagrangian series of position, temperature, salinity, etc. for later use by an IBM without any further spatial tracking (Fig. 5a).

iv) Run the hydrodynamic model once, storing all relevant gridded velocity, temperature, and salinity fields at each model gridpoint and time step, for subsequent use by an IBM which performs its own spatial tracking (Fig. 5b).

All approaches have both advantages and drawbacks. Method i) is at present computationally prohibitive for most three-dimensional hydrodynamic models, as it requires a potentially large number of runs. Method ii) reduces computation of the physical model, but requires huge amounts of computer memory to keep track of all the biological realizations simultaneously. Method iii) is efficient and allows for many individuals, but eliminates the possibility of adding individual behaviors (in particular, vertical locomotion) which vary as a function of the individual's unique, partly random, life history. In other words, any realization of the IBM in method iii) cannot feed back on the float tracks, which are determined by the hydrodynamic model alone. Method iv) offers the possibility of feedback, but requires potentially huge amounts of data storage to accommodate the three-dimensional circulation model output. Hence we propose a modification of iv):

iv-a) Run the hydrodynamic model once, storing suitably low-pass filtered (high frequencies removed), decimated time series of all relevant gridded velocity, temperature, and salinity fields at each model gridpoint, for subsequent use by an IBM which performs its own spatial tracking.

Note that the characteristics of the filter used in iv-a), and the interval for subsampling, depend on what portion of the full spectrum of motions is considered "expendable" for purposes of spatial tracking. For example, time-filtering would introduce a bias if locomotion of individuals were somehow correlated with higher frequency current fluctuations (e.g. swimming in response to tides or internal waves), that had been removed. In general, primitive equation hydrodynamic model output contains far less energy at superinertial frequencies than is observed in real oceanographic measurements, and such fluctuations that do exist in the model are poorly correlated with the real data. Nevertheless, subsampling without filtering, while partially adequate for spatial tracking, could lead to serious aliasing errors when much high frequency energy is present. We suggest filtering be done to eliminate this poorly correlated high frequency band, but retain significantly correlated subinertial motions. The significance of tides is context-dependent, and needs to be considered carefully.

Consider how spatial tracking in method iv-a) might be efficiently accomplished with the filtered fields, and whether any of the information lost in the filtering process might be added back in by the tracking algorithm. Neutrally buoyant floats have been tracked in many numerical ocean models. Instantaneous float velocities are generally derived by interpolating between gridpoints to the current location of each float. Generally the floats are advanced in time using the same time step as is used for solving the governing hydrodynamic equations, but this need not be the case, and especially not in our situation where previously stored, filtered velocity fields are used. How, then, to choose a time step for float tracking which is sufficiently short but not computationally prohibitive? One requirement for successful float tracking is that the Lagrangian decorrelation time T_L be much longer than the time step dt used for float tracking:
 

    (4)
 

The property T_L is defined as follows. When moving within a turbulent flow field, particle velocities eventually become decorrelated from their starting velocities. This loss of "memory" can be characterized by the integral of the lagged correlation function:
 
 

    (5)

where

    (6)
 

here u'u' is the mean turbulent kinetic energy of the flow, u'_i(t) is the velocity of the particle i at time t, Tmax is a suitably long time used for averaging and {} denotes an ensemble mean over all particles.

When eq. (4) is satisfied, the moving float (individual) can be reasonably approximated to travel at constant velocity for a period of time which is longer than the time between updates of the velocity field used for tracking. When this not the case, float tracks will diverge from "true" particle paths. (Haidvogel [1982] has also demonstrated that adequate spatial resolution is crucial for accurate float tracking; with overly coarse resolution, nonlinear interactions present in the hydrodynamic equations will be improperly represented and floats will diverge from their true path as a result).

We generally expect time-filtering of velocities to increase T_L. This is helpful in allowing a longer time step for tracking, as well as reducing the size of the velocity file, if the higher-frequency variability can be reasonably sacrificed. If tides are a critical aspect of the individual's life history, there may be no easy way to avoid using unfiltered velocities, and a short tracking time step. Note, however, that some information about tides could be stored in compact form (amplitude and phase) for later use; such reconstructed tidal velocities could be added to filtered (subtidal) velocity fields directly in the tracking algorithm, potentially reducing storage requirements.

The model used in the present study is a prognostic, rigid-lid, hydrostatic, three-dimensional, primitive equation model of velocity and salinity fields in the northern Gulf of Alaska. Details of model construction and sensitivity analyses may be found in Hermann and Stabeno (1996). Stabeno and Hermann (1996) demonstrated that the model yields realistic mean flows and eddy statistics (frequency of passage of mesoscale eddies) in the vicinity of Shelikof Strait. There are no explicit tides in this model, but some high-frequency (near-inertial) oscillations are included in the dynamics. Boundary conditions were provided using a telescoped grid with a toroidal, re-entrant channel, as discussed in Hermann and Stabeno (1996). The finely resolved portion of the model domain (4 km mean spacing) is shown in Fig. 1; the finely resolved area extends from the entrance of Shelikof Strait to slightly westward of the Shumagin Islands, and southward to just beyond the Alaskan Stream.
 

3.b. NPZ-SPEM coupling

The characteristic length and time scales of biological properties are strongly influenced by the corresponding scales of their physical environment, but are not identical to them. Ideally we would like to resolve both fields so finely in space and time that all relevant scales are included, but this is rarely possible due the cost involved. As a compromise, we chose different horizontal and vertical scales for the Eulerian physical and biological models. For NPZ dynamics, we generally require a finer vertical scale than for the physical model due to the sensitivity of primary production to changes in mixed layer depth. This entails a sacrifice of horizontal resolution. The character of mesoscale eddies in the vicinity of Shelikof Strait suggests a horizontal spacing of ~12 km as a desirable goal; note, however, that the useful resolution for our NPZ model is fundamentally dependent on the characteristic scales of the modeled circulation, rather than the characteristic (and typically finer) scales of the circulation itself. As a first attempt, we chose to use a model with 1 meter resolution in the vertical, and 20 km resolution (5x the physical model) in the horizontal (Fig. 6). Vertical coverage is from 0 to 100 m depth. Background horizontal and vertical diffusivities in the NPZ model were set at 10 and 10^-5 m² s^-1, respectively.

As with the IBM, we use time-filtered, stored output from the physical model as input for the NPZ model. Time-filtered barotropic streamfunction output from the physical model is subsampled to the grid of the NPZ model. The horizontal gradients of the barotropic streamfunction, divided by the water column depth (or by 100 m in regions of deeper bathymetry), yields depth averaged velocities at each horizontal location. These are used to calculate the horizontal divergence and, subsequently, the vertical velocities. This approximation distorts the true vertical shear and vertical velocities, but avoids spurious (and potentially very large) convergence which typically results from subsampling a full three-dimensional velocity field. Errors in the computed vertical velocity field are likely not too significant for the present case; vertical upwelling of nutrients strongly affects production in some parts of the ocean, but in the Gulf of Alaska downwelling predominates during most of the year. Horizontal advection and mixed layer deepening are hence the more significant sources of nutrients to feed production in this area.

The vertical structure of the NPZ model assumes a mixed layer of variable depth with uniform properties, and a stratified region beneath (down to 100m), where mixing is based on the spatially constant background diffusivity. Mixed layer depth (MLD) time series were generated for the NPZ model using the circulation model's density profiles at the horizontal locations of the NPZ grid. A simple algorithm is presently employed, where MLD is defined as that depth where the salinity field (the primary determinant of density in this region during spring) assumes a value 0.05 psu greater than the surface value. Clearly this is a crude parameterization, hence future circulation modeling in this region will include explicit mixed layer physics.

Since we have used stored, time-filtered and decimated streamfunction fields to drive the NPZ model, we must consider the effects of that filtering on the prey dynamics. Float tracking experiments suggest an acceptable loss of temporal information for our rigid lid circulation model with no tides (see Results section). The effect of varying spatial resolution on the NPZ fields depends on the magnitude of the advection terms, relative to the local nonconservative (biological) dynamics. If biological dynamics predominate over advective change, then coarse resolution of the advective terms is a permissible compromise. We may ultimately be able to invoke spatially dependent diffusivities, to partially compensate for the loss of spatial resolution. In Shelikof Strait biological components with the fastest turnover times (phytoplankton in our NPZ model) will be the least succeptible to aliasing from the choice of a coarse spatial grid, while those with longer turnover times will be most affected. Napp et al. (1996) have demonstrated that advection plays a large role in determining the spatial patterns of phytoplankton and zooplankton ovserved in Shelikof Strait. Future iterations of our coupled SPEM-NPZ models will examine the consequences of decreasing the spatial scale.

The NPZ model domain is a subset of the physical model, and needs boundary conditions for each of its biological variables. In the interior, the balance of eq (1) holds; separating out the spatial terms yields:

     (7)

where B is an NPZ model variable (e.g. Phytoplankton), Lh is the horizontal gradient operator, u^h ,v^h ,k^h are the horizontal velocities and eddy diffusivity, w is the vertical velocity, k^z is the vertical eddy diffusivity and f_biol represents all nonconservative biological dynamics. Perfect knowledge of upstream values would permit direct specification of B at the upstream boundary; instead, we are limited to several springtime measurements per year taken 7 km apart along a transect accross the lower exit region of Shelikof Strait (Napp et al., 1996; Incze et al., 1997; Incze and Ainaire, 1995). Hence, we provide values at the upstream end by running the model in one-dimensional mode (no advection or horizontal diffusion), with assimilation of observed values at those locations.

    (8)

where _nudge is the nudging constant and B_data is the available data. Nudging is currently implemented for Pseudocalanus nauplii in the mixed layer; we use _nudge = 0.5 d^-1 during periods when data were available (mid-April through mid-May). Boundary values computed using eq. (8) are advected into the interior, driven by the velocity field and the horizontal gradient term of eq. (7). The assimilation, through nudging, of the boundary data is a way of compensating for our ignorance about the true boundary values; eq (8) may be thought of as a dynamically consistent interpolator of boundary values. It is similar in spirit to a class of physical boundary conditions which combine simplified physics (e.g. wave dynamics) with nudging towards observed values.

At the downstream boundary, a zero gradient condition was imposed for all biological variables.
 

     (9)
 

3.c. IBM-NPZ coupling

Some of the issues discussed above for the proper application of stored circulation fields to the IBM hold for the application of Eulerian biological fields to the IBM as well. Ideally we could interpolate perfectly accurate Eulerian prey values onto the location and time of the tracked individuals of the Lagrangian IBM. In practice we must settle for statistically correct prey fields, driven in part by the statistically correct velocity fields from the physical model. The mapping of prey values to individual fish locations is straightforward, and is performed in the IBM using simple linear interpolation of NPZ model output in time and space.

Boundary conditions for the IBM are presently as follows. 1) When an individual is located outside the NPZ model domain, we apply the prey values of the NPZ boundary point nearest that individual's present location. 2) When an individual is advected outside the finely resolved circulation model domain, it is removed from the population. We also trap individuals at the downstream (southwestern) end of the physical domain, to prevent their re-entrance at the upstream end.
 

4. Results

4.a. Float tracking experiments

To help address the accuracy of float tracking with pre-stored, filtered and decimated velocity fields, versus "direct" tracking with the unfiltered field, we present results from a simulation of 1987 currents and float tracks in the Gulf of Alaska. Details of model forcing for this year are available in Hermann et al. (1996). In Fig. 7a we illustrate the results of float tracking by direct updating (case a) in the hydrodynamic model; that is, ten floats are tracked using the instantaneous velocity produced by the model at each time step with a 4th-order Runge-Kutta scheme. The time step for this model is .0375 hr. Interpolation of velocities between model gridpoints is performed as described in Hoffmann et al. (1989); briefly, this entails linear interpolation in the horizontal and spectral (here, Chebyshev) expansion in the vertical. Initial positions span a northwest-southeast line between the Alaska Peninsula and the northwestern corner of Kodiak Island, at 40 m depth. Many of the floats circulate around and through a cyclonic-anticyclonic eddy pair, just downstream from the exit of Shelikof Strait. Note how some of the floats, released in deep regions, cross bathymetric contours to the shallower continental shelf.

In Figs. 7b-7d we illustrate the same set of floats, this time tracked with pre-stored, filtered and decimated velocity fields. These fields are produced by applying a cosine-Lanczos 30-hour low-pass filter to the original model time series during execution of the run, and storing the filtered result once per day. This filtering process eliminates all of the near-inertial, internal wave oscillations. We wish to see if accurate float tracking can be achieved with time steps longer than the .0375 hr of the "direct" unfiltered method, and whether the resulting tracks compare favorably with those derived using the unfiltered method. Float tracking is attempted with three different values of the time step: b) dt =1 hour; c) dt = 0.2 hour; d) dt = .04 hour. Linear interpolation in time is used to provide velocity values in between the stored daily values; spatial interpolation is performed as in case a. In each case floats are seeded at the same locations and depths as in case a. In case b (Fig. 7b), the broadly defined paths of the float tracks from case a (Fig. 7a) are reasonably well reproduced, but the details of the paths differ as the floats are advected downstream. Can a shorter time step do better? Case c (Fig. 7c) differs from case b primarily in the improved fit with features in the vicinity of the Shumagin Islands (55° N, -160° E). Case d (Fig. 7d), with highest temporal resolution, exhibits no substantial difference from case c.

In Fig. 8 we plot the root-mean-squared (rms) horizontal displacement of the floats in case b from their corresponding locations in case a, the displacement of floats between cases b and c, and the displacement between cases c and d. During the first four days of simulation, case b generates float tracks within ~2 km of their locations in case a. Thereafter case b results diverge from case a, to a mean displacement of ~40 km (the typical width of the deep strait and sea valley in this area). Near day 30, the mean displacement increases once more, ultimately achieving values of ~140 km. Cases b and c diverge significantly from each other after day 30 while cases c and d are essentially identical through day 60.

Apparently, by shortening the time unit from 1 hour to 0.2 hour, we sufficiently converge on the "true" solution for the filtered velocity field, that is, we have achieved a time step sufficiently shorter than the Lagrangian decorrelation time for the filtered velocities to yield an accurate result. This justifies the use of a time step of 0.2 hour in the float tracking algorithm described in Hinckley et al (1996); the same time step and spatial interpolation method is used to generate the results shown in Section 4.b.

Fewer of the float tracks cross into shallow shelf areas when filtered velocities are used, compared to the unfiltered technique. In future studies, this deficiency could be addressed to some extent by adding a suitable random walk component to the tracked floats, an equivalence first noted by Taylor (1921). We caution, however, that the diffusivity lost by temporal filtering may be spatially dependent, and that spatial variation in the magnitude of a random walk can pose special difficulties (Hunter et al., 1993; Holloway, 1994).

The importance of vertical position in determining the paths of individuals was noted by Hinckley et al. (1996). In Fig. 9, we underscore this point by plotting float paths as in case c, but with all floats constrained to remain at 40 m depth, rather than being free to follow the flow field in three dimensions. The influence of the eddy field is markedly reduced by this constraint; floats do not circulate around the eddy pair as they did in other cases.
 

4.b. Coupled run output, timing and storage

Some illustrative results from the coupled models for mid-June, 1987 are shown in Figs. 10a-c. A total of 1600 individuals are used in the IBM, with mortalities set to zero. A cyclonic eddy is located in the deepest part of the sea valley, near 56.5 ° N, -156 ° E (Fig. 10a). Total numbers of all copepodite stages of Pseudocalanus are shown for zooplankton (Fig. 10b). Positions of individual fish are shaded according to length (Fig. 10c). The individual lengths vary, due to the different prey fields encountered during their life history; larvae greater than 8 mm ingest a mixture of nauplii and copepodites. A cluster of medium-sized larvae (8-16 mm) is trapped within the eddy. A tongue of low zooplankton concentration is being advected around the southern rim of the eddy, and shorter fish are associated with that feature. A tongue of higher zooplankton concentration, advecting around the northern rim of the eddy, has fewer but generally longer fish associated with it. Regions downstream and shoreward of this feature also exhibit greater lengths.

While some patchiness of fish size can be produced through purely kinematic effects related to the time of spawning (e.g. a cluster of late-spawned, smaller individuals trapped in an eddy), clearly the spatially variable prey field contributes to the patchiness of fish attributes in this model. The use of a dynamic prey field also results in different mean growth and survival for groups of individuals spawned at different places and times (Hinckley et al., 1988, this volume).

All model results described here were generated using a CRAY Y-MP vector processing supercomputer. Statistics of runs and output files are listed in Table 1. As shown, the physical model consumes the largest amount of computer time. Indeed, this is the primary motivation for decoupling the physical and biological runs; by decoupling we allow for much more testing and sensitivity analysis of the biological models. The output files from the physical model, while large, are not prohibitively so, and can be transferred to local workstations for detailed analysis and plotting. Statistics of the IBM are dependent on the number of individuals followed.
 

5. Conclusions and future work

Spatially explicit biophysical models, spanning multiple trophic levels, can be a powerful tool for probing hypotheses about the early life history of marine fish. Throughout our ongoing work with coupled hydrodynamic, population (IBM), and trophodynamic (NPZ) models, we have attempted to balance desirable features with practical limits on computing power. Our experience thus far has led to the following conclusions. 1) Storage of pre-filtered and decimated model velocity and scalar fields, for later use by Lagrangian biological models, is both feasible and efficient using presently available computers. Supercomputing platforms with reasonable large (several gbyte) storage capacity are necessary for much of this work, however. 2) The results are reasonably accurate given a sufficiently small time step for spatial tracking; however, care must be taken to consider the information lost through filtering. Practical techniques exist to recover some of the lost information, such as tides. 3) Spatially varying prey fields, derived using a traditional Eulerian approach, can be applied in a straightforward manner to spatially explicit (Lagrangian) IBM's, and should typically enhance the variance of individual attributes. 4) Boundary conditions for Eulerian prey models should benefit from a careful use of data assimilation, when field data are available.

As our models evolve to greater complexity, there is a continuing need for calibration and sensitivity analyses of the coupling methods which link them. For example, there is a strong need for further sensitivity analyses of grid spacing in the NPZ model. Equally important is the need for further comparison of three-dimensional NPZ model output with chlorophyll and zooplankton data. For chlorophyll, satellite-derived estimates may be useful for comparing at least broad spatial patterns, and possibly higher-order spatial statistics, with those produced from the model. These and other issues will be explored in future reports.
 

Acknowledgements

This is contribution number 1920 from NOAA/ Pacific Marine Environmental Laboratory, contribution number 318 to NOAA's Fisheries Oceanography Coordinated Investigations, and contribution number 475 from the Joint Institute for the Study of the Atmosphere and Oceans under cooperative agreement NA90RAH00073.
 

Batchelder, H. P. and Williams, R. (1995). Individual-based modeling of the population dynamics of Metridia lucens in the North Atlantic. ICES J. Mar. Sci. 52: 469-482.
 

Botsford, L. W., Moloney, C. L., Hastings, A., Largier, J. L., Powell, T. M., Higgins, K., and Quinn, J. F. (1994). The influence of spatially and temporally varying oceanographic conditions on meroplanktonic metapopulations. Deep-Sea Res. 41 107-145.
 

Crowder, L. W., Rice, J. A., Miller, T. J., and Marschall, E. A. (1993). Empirical and theoretical approaches to size-based interactions and recruitment variability in fishes. In DeAngelis, D. L., Gross, L. J. (eds.) Individual-Based Models and Approaches in Ecology, Chapman and Hall, pp. 237-255.
 

DeAngelis, D. L., Godbout, L., and Shuter, B. J. (1991). An individual-based approach to predicting density-dependent dynamics in smallmouth bass populations. Ecol. Model. 57 91-115.
 

Dimou, K. N. and Adams. P E. (1991). Representation of sources in a 3-D Eulerian-Lagrangian mass transport model. In Wrobel, L. C., Brebbia, C. A. (eds.) Water Pollution Modeling, Measurement and Prediction 1991, pp. 251-264.
 

Falkowski, P. G., and Wirick, C. D. (1981). A simulation model of the effects of vertical mixing on primary productivity. Mar. Biol. 65 69-75.
 

Flierl, G. R., and Davis, C. S. (1993). Biological effects of Gulf Stream meandering. J. Mar. Res. 51 529-560.
 

Franks, P. J. S. and Walstad, L. J. (1997). Phytoplankton patches at fronts: a model of formation and response to wind events. J. Mar. Res. 55: 1-29.
 

Frost, B. W. (1987). Grazing control of phytoplankton stock in the open subarctic Pacific Ocean a model assessing the role of mesozooplankton, particularly the large calanoid copepods Neocalanus spp. Mar. Ecol. Prog. Ser. 39 49-68.
 

Frost, B. W. (1993). A modelling study of processes regulating plankton standing stock and production in the open subarctic Pacific Ocean. Prog. Oceanogr. 32 17-56.
 

Fukumori, I. and Malanotte-Rizzoli, P. (1995). An approximate Kalman filter of ocean data assimilation: an example with an idealized Gulf Stream model. J. Geophys. Res. 100(C4): 6777-6793.
 

Haidvogel, D. B., Wilkin, J. L. and Young, R. (1991). A semi-spectral primitive equation ocean circulation model using vertical sigma and orthogonal curvilinear coordinates. J. Comput. Phys. 94: 151-185.
 

Haidvogel, D. (1982). On the feasibility of particle tracking in Eulerian ocean models. Ocean Modelling 45 4-9.
 

Hermann, A. J. and Stabeno, P. J. (1996). An eddy resolving model of circulation on the western Gulf of Alaska shelf. I. Model development and sensitivity analyses. J. Geophys. Res.101(C1): 1129-1149.
 

Hinckley, S., Hermann, A. J., and Megrey, B. A. (1996). Development of a spatially explicit, individual-based model of marine fish early life history. Mar. Ecol. Prog. Ser.139: 47-68.
 

Hinckley, S., Hermann, A. J., Meir, K. L. and Megrey, B. A. (1998). The importance of spawning location and timing to successful transport to nursery areas: a simulation modeling study of Gulf of Alaska walleye pollock. This volume.
 

Hoffmann, E. E., Hedstrom, K. S., Moisan, J. R., Haidvogel, D. B., and Mackas, D. L. (1991). The use of simulated drifter tracks to investigate general transport patterns and residence times in the coastal transition zone. J. Geophys. Res. 96 15,041-15,052.
 

Holloway, G. (1994). On modeling vertical trajectories of phytoplankton in a mixed layer. Deep-Sea Res. II 41 957-959.
 

Hunter, J. R., Craig, P. D., and Phillips, H. E. (1993). On the use of random walk models with spatially variable diffusivity. J. Comput. Phys. 106 366-376.
 

Incze, L. S., Siefert, D. W., and Napp, J. M. (1996). Mesozooplankton of Shelikof Strait, Alaska: abundance and community composition. Cont. Shelf Res. 17: 287-305.
 

Jamart, B. M., Winter, D. F., Banze, K., Anderson, G. C., and Lam, R. K. (1977). A theoretical study of phytoplankton growth and nutrient distribution in the Pacific Ocean off the northwestern U.S. coast. Deep-Sea Res. 24 753-773.
 

Jamart, B. M., Winter, D. F., and Banze, K. (1979). Sensitivity analysis of a mathematical model of phytoplankton growth and nutrient distribution in the Pacific Ocean off the northwestern U. S. coast. J. Plank. Res. 1 267-290.
 

Kendall, A. W., Jr. and Nakatani, T. (1992). Comparisons of early life history characteristics of walleye pollock Thergra chalcogramma in Shelikof Strait, Alaska and Funka Bay, Hokkaido, Japan. Fish. Bull. 90: 129-138.
 

Kendall, A. W., Jr., Schumacher, J. D. and Kim, S. (1996). Walleye pollock recruitment in Shelikof Strait: applied fisheries oceanography. Fish. Oceanogr. 5 (suppl.): 4-18.
 

Lande, R. and Lewis, M. R. (1989). Models of photoadaptation and photosynthesis by algal cells in a turbulent mixed layer. Deep Sea Res. 36: 1161-1175.
 

Lewis, C. V. W., Davis, C. S., and Gawarkiewicz, G. (1994). Wind forced biological-physical interactions on an isolated offshore bank. Deep-Sea Res. 41 51-73.
 

MacKenzie, B. R. and Leggett, W. C. (1993) .Wind-based models for estimating the dissipation rates of turbulent energy in aquatic environments: Empirical comparisons. Mar. Ecol. Prog. Ser. 94: 207-216.
 

Megrey, B. A. and Hinckley, S. (1998). The effect of turbulence on feeding of larval fishes: a sensitivity analysis using an individual based model. This volume.
 

Moisan, J. R. and Hofmann, E. E. (1996a). Modeling nutrient and plankton processes in the Coastal Transition Zone. Part 1: A time- and depth-dependent model. J. Geophys. Res. 101(C10): 22647-22676.
 

Moisan, J. R. and Hofmann, E. E. (1996b). Modeling nutrient and plankton processes in the Coastal Transition Zone. Part 3: Lagrangian drifters. J. Geophys. Res. 101(C10): 22693-22704.
 

Moisan, J. R., Hofmann, E. E. and Haidvogel, D. B. (1996). Modeling nutrient and plankton processes in the Coastal Transition Zone. Part 2: A three-dimensional physical-bio-optical model. J. Geophys. Res. 101(C10): 22677-22691.
 

Napp, J. M., Incze, L. S., Ortner, P. B., Siefiert, D. L. W., and Britt, L. (1996). The plankton of Shelikof Strait, Alaska: standing stock, production, mesoscale variability and their relevance to larval fish survival. Fish. Oceanogr. 5 (suppl. 1): 19-38.
 

Platt, T., and Gallegos, C. L. (1980). Modeling primary production. In Falkowski, P. G. (ed.) Primary Productivity in the Sea, , Plenum Press, New York, 531 pp.
 

Rose, K. A., Christensen, S. W., and DeAngelis, D. L. (1993). Individual-based modeling of populations with high mortality: a new method for following a fixed number of model individuals. Ecol. Model. 68 273-292.
 

Scheffer, M., Baveco, J. M., DeAngelis, D. L., Rose, K. A., and Van Nes, E. H. (1995). Super-individuals: a simple solution for modelling large populations on an individual basis. Ecol . Model. 80: 161-170.
 

Schumacher, J., and Kendall, A. W. Jr. (1995). An example of fisheries oceanography: Walleye pollock in Alaskan waters. U. S. Nat. Rep., Int. Union Geod. Geophys. 1991-1994, Rev. Geophys. 33: 1153-1163.
 

Stabeno, P. J., and Hermann, A. J. (1996). An eddy resolving model of circulation on the western Gulf of Alaska shelf. 2. Comparison of results to oceanographic data. J. Geophys. Res.101(C1): 1151-1161.
 

Taylor, G. I. (1921). Diffusion by continuous movements. Proc. London Math. Soc. (ser. 2) 10: 196-212.
 

Werner, F. E., Page, F. H., Lynch, D. R., Loder, J. W., Lough, R. G., Perry, R. I., Greenberg, D. A. and Sinclair, M. M. (1993). Influences of mean advection and simple behavior on the distribution of cod and haddock early life stages on Georges Bank. Fish. Oceanogr. 2: 43-64.
 

Werner, F. E., Perry, R. I., Lough, R. G. and Naimie, C. E.. (1996). Trophodynamic and advective influences on Georges Bank larval cod and haddock. Deep-Sea Res. 2 (Top. Stud. Oceanogr.) 43: 1793-1822.
 

Woods, J., and Barkmann, W. (1993). Diatom demography in winter - simulated by the Lagrangian Ensemble method. Fish. Oceanogr. 2 202-222.
 

Woods, J. D., and Onken, R. (1982). Diurnal variation and primary production in the ocean - preliminary results of a Lagrangian ensemble method. J. Plankton Res. 4 735-756.
 

Wroblewski, J. S. (1977). A model of phytoplankton plume formation during variable Oregon upwelling. J. Mar. Res. 35 357-394.

Figure 1. Summary of early life history for walleye pollock spawning at the exit of Shelikof Strait. Primary regional currents are shown in inset. Subsequent to spawning, hatched larvae from eggs drift southwest with the Alaska Coastal Current; some are transported up onto the shallow, near-coastal shelf area while others are transported out into the deeper waters of the Alaskan Stream. Dashed line indicates the boundary of the finely resolved physical model domain (~4 km mean spacing).
 

Figure 2. Summary of the individual-based model structure.
 

Figure 3. Summary of the NPZ model structure.
 

Figure 4. Summary of coupling among the physical and biological models.
 

Figure 5. Two schemes for implementing a Probabilistic-Lagrangian biophysical model of the early life stages of marine fish. Schematic a) represents scheme iii) described in the text; b) represents scheme iv) described in the text. "SPEM" refers to the Semispectral Primitive Equation Model of Haidvogel et al. (1991), used for generating ocean current, salinity and temperature fields. "IBM" refers to the Individual-Based Model of marine fish. [.]_i(t) refers to the value of some physical or biological quantity associated with individual i.
 

Fig. 6. Overlay of SPEM (+) and NPZ (o) model grids. Horizontal grid spacing of the NPZ grid is 5x larger than the SPEM grid. Triangles () indicate upstream locations where Pseuodocalanus values are nudged to available data.
 

Figure 7. a) Float tracks generating by directly updating float positions at each time step in a hydrodynamic model (SPEM) simulation of 1987. Ten floats, released along a cross-shelf line at 40m depth on DOY 135 for 1989, are shown. Tracks are coded by time interval: DOY 135-165 (black); DOY 165-195 (red); DOY 195-225 (green); DOY 225-255 (dark blue); DOY 255-270 (light blue). Model bathymetry (meters) is shown in greyscale.b) Float tracks generated by updating float positions using previously stored, filtered, daily velocity fields, with dt=1.0 hour for the float tracking algorithm. c) Results as in (b) using dt=0.2 hour for the float tracking algorithm. d) Results as in (b) using dt=0.04 hour for the float tracking algorithm.
 

Figure 8. Root-mean-square relative displacement of floats tracks among the various tracking experiments. Black line - filtered velocities with dt=1.0 hour (case b) versus "direct" tracking (case a). Red line - filtered velocities with dt=1.0 hour (case b) versus dt=0.2 hour (case c). Green line - filtered velocities with dt=0.2 hour (case c) versus dt=0.04 hour (case d). Mean relative displacement in each case is plotted as a function of time since the beginning of the experiment (DOY 105, 1987).
 

Fig. 9. Float tracks generating by directly updating float positions at each time step in a hydrodynamic model (SPEM) simulation of 1987, where all floats are constrained to remain at 40m depth. Initial float positions and coding of tracks are as in Fig. 6.
 

Fig. 10. Results of coupled models for mid-June, 1987. a) Surface velocities (m s^-1) from the circulation model. b) NPZ model output. Shading represents total concentration of all copepodite stages of Pseudocalanus spp. in the mixed layer (# m^-3). Bathymetry is shown by solid lines (m) and barotropic streamfunction is shown by dashed lines (contour interval is 0.1 x 10^6 m^3/s coastward of the shelf break and 1.0 x 10^6 seaward. c) IBM output. Positions of individual pollock are shaded by length (mm).
 

Table 1. Timing and Storage statistics for the three coupled models.
 

Model type	Circulation  (SPEM)	Prey  (NPZ)	Pollock  (IBM)
Model size	257x97x9 gridpoints	20x20x100  gridpoints	1600  individuals
Resolution	4 km horiz  5 m vert	20 km horiz  1 m vert	4 km horiz  5 m vert
Daily output	velocity  salinity  temperature	Nitrate concentration  Phytoplankton biomass  Zooplankton biomass  (14 stages/species)	length  weight  position
Length of simulation (days)	270	95	210
Output file size (gbytes)	1.0	0.5	.02
CRAY YMP CPU (hours)	125	5	5