The ecological risk conceptual model was converted into a simulation model that quantifies the toxicological risks to NKI sea otters from potential exposure to remnant EVOS-SSOR. The quantitative model simulates the co-occurrence of sea otter pits and SSOR based on realistic spatial relationships. We used the object-oriented simulation language StellaTM
(ver. 9.0.2; copyright isee systems, inc. [www.iseesystems.com]) to simulate each exposure pathway in the conceptual model. Six primary sources of empirical data were used in the model: (1) our observational data on foraging, diving, and dietary attributes of each class of PWS sea otters, discussed above (), supplemented with other data on sea otter characteristics (); (2) characterization of the spatial distribution of sediment and habitat types of the NKI Study Area () based on our shoreline survey, which was conducted by skiff following the shoreline ~50 m offshore at low tides in June 2008; (3) the most recent NOAA survey data on the spatial distribution of SSOR on the shorelines of NKI (Michel et al. 2006
; Short et al. 2006
); (4) the chemical characterization by NOAA of 41 PAHs in the SSOR on the shorelines of NKI (EVOSTC 2008
); (5) the survey in summer 2006 by Boehm et al. (2007a)
of sea otter pits on shorelines of NKI during particularly low tides; and (6) chemical characterization of 41 PAHs in prey based on mussels collected in 1998–2002 (Boehm et al. 2004
); clams collected in 2002–2004 (Neff et al. 2006
); crabs, whelks, and worms collected in 2002–2004 (Neff et al. 2006
); and seawater collected in 2005 (Boehm et al. 2007b
). Note that in all cases, measured PAH concentrations reported as 0 or non-detected were assigned here the value of ½ Method Detection Limit.
Model Parameterization for Simulating Co-occurrence of Sea Otter Pit and SSOR
For a few parameters in the model, no empirical data exist, requiring those parameterizations to be based on informed expert judgment (see Acknowledgments for list of expert participants in modeling workshops). These parameters included: (a) grooming efficiency for removing particles (75% for HOR); (b) fraction of sediment particles sticking to fur (5% for HOR); and (c) thickness of SSOR-coating on paws and clams (33x thickness of HOR oil film [characterized as dark colors film; see HAZMAT 1996
]). Each of the expert-judgment-based parameters were assigned conservatively and were either subject to parameter-sensitivity analyses (discussed later) or were assigned the most conservative possible number (e.g.,
100% efficiency for ingestion of oil-phase SSOR coating on paws and clams; zero dilution of oil-phase SSOR released from pits as it floats to become surface film, 100% of which coats the sea otter's fur).
The quantitative model predicts daily doses of each of 41 PAH analytes assimilated by an individual sea otter subject to specific activities each hour of the day, utilizing the individual-based modeling (IBM) approach (DeAngelis and Gross 1992
), in which environmental constraints parameterized by empirical data provide the framework in which sea otter behavior occurs. The model predicts exposures and effects on an individual sea otter under the paired situation of: (1) an at-risk sea otter,
which has potential direct PAH exposure from excavating a pit that intersects SSOR, consumption of diet and exposure to sediments in the vicinity of SSOR, and exposure to background levels of PAHs in diet, drinking water, and sediments from areas separate from the SSOR; and (2) a not-at-risk sea otter,
which experiences identical situations except exposure is only to background levels of PAHs in diet, drinking water, and sediments. Many model parameters are age- or gender-specific; consequently, separate models simulate each of seven classes of sea otters: (1) older pups (ages 2–6 months); (2) juveniles (ages 6–12 months); (3) subadults (ages 1–3 years); (4) adult males in the territorial phase (limited to a single bay, BI); (5) adult males in the non-territorial phase (allowed to roam over all three bays); (6) adult females with dependent pups (also limited to a single bay, BI); and (7) adult females without pups. Each adult phase is assumed to extend over a 6-month period. Dependent pups ≤2 months were not modeled separately as they are unable to digpits and fully rely on their mother's milk for food (Payne and Jameson 1984
), which apparently does not transfer significant quantities of PAHs (e.g., Bulder et al. 2006
A major component of the IBM model is a series of stochastic functions that incorporate variability for the sea otters and the environmental conditions in which they actually live. Each distribution, derived in all cases from empirical data, is sampled using a pseudo-random number generator, that is, meeting all criteria for randomness following a particular sequence initiated by a “seed,” thereby allowing reconstruction of specific runs and parameter-sensitivity analyses. The model's stochastic functions assign the sea otter to a specific bay daily, based on the frequency distribution of shoreline lengths (pBAY
) from the SCAT maps of NKI, and assign the following parameters hourly (): (1) assignment to an SSOR-subdivision or a non-SSOR-subdivision, based on the bay-specific frequency distribution of SSORand non-SSOR subdivisions (PSSOR subdivion
); (2) assignment to the ITZ or the STZ, based on class-specific frequency distributions from our sea otter observational data (psea otter
ITZ); (3) assignment to a site with sediments that potentially could have SSOR (bedrock sediments cannot have SSOR), based on the frequency distribution of sediment types from our shoreline survey discussed earlier (pSSOR·potential sediments
); (4) assignment to sediments that potentially support clams and thus could be excavated by sea otters (i.e.,
fine-grained sediments), based on the same shoreline survey (Pclam·potential sediments
); (5) assignment to an SSOR field or not, based on the SSOR frequency distribution from the 2003 NOAA intense-survey sites (here the term field
means an area comparable to the cross-ITZ 12.5 m-wide swath sampled in that survey) (EVOSTC 2008
) (pSSOR field
); (6) probability of a patch of SSOR being within the LITZ and thus potentially containing sufficient clams that sea otters might dig, based on the frequency distribution from the NOAA SSOR surveys (pSSOR LITZ
); (7) assignment of a dive resulting in a pit or not, based on the frequency distributions of infauna prey in the ITZ from our feeding observational database (psea otter ITZ pit
); and (8) probability that the hour is within the feeding period of the day, based on our feeding observational database (Pfeeding
activity). Thus within a particular bay, the effective probability of co-occurrence of a sea otter pit and a patch of SSOR () is a product of the individual probabilities, calculated as:
Effective Probability of Co-occurrence of Sea Otter Pit and a Patch of SSOR
Note that the model does not literally assign this effective probability value (except for some sensitivity analyses, discussed later); rather, co-occurrence emerges from the hour-by-hour assignment of each stochastic parameter based on the data distributions. If a patch of SSOR is intersected during a time-step, then the model assigns the SSOR category (e.g.,
HOR, MOR) based on the frequency distribution of NKI SSOR from the NOAA surveys. The model then assigns a specific TPAH concentration to that particular SSOR category, randomly sampling from a lognormal distribution fit to the SSOR PAH data (EVOSTC 2008
), and the mass of each analyte is calculated based on the relative fraction of TPAH of each analyte. If no co-occurrence occurs during a time-step, the model assigns only PAH concentrations from background sources (e.g., diet, seawater, and STZ and non-SSOR-ITZ sediments). Similarly, for the not-at-risk sea otter, all SSOR cooccurrence algorithms are bypassed.
The diet consumed by a sea otter class is assigned based on our feeding observation data (). The foraging values used to parameterize the model are similar to those reported for sea otters occupying similar habitats with similar histories of occupation elsewhere. However, PWS sea otters feed more in the STZ than the ITZ because clams, their primary food in PWS (Calkins 1978
; Doroff and Bodkin 1994
; Johnson and Garshelis 1995
), are more abundant in the large benthic STZ compared with the narrow ITZ (Dean et al. 2002
). Similar clam-dominated diets have been observed for sea otters in soft-bottomed habitats in Kodiak (Doroff and DeGange 1994
) and Washington State (Laidre and Jameson 2006
). The estimated daily mean number of foraging dives across sea otter classes in PWS ranged from 348 to 684 dives-day–1
* feeding period; ), 2–5x higher than recorded for sea otters in southeast Alaska (Bodkin et al. 2004
). This is as expected because dive frequency is inversely related to water depth, and the feeding areas in PWS are shallower than the foraging areas reported in southeast Alaska, where a substantial proportion of dives were > 60 m deep (Bodkin et al. 2004
). Rates or proportions of intertidal foraging by sea otters have not been reported for other areas. However, we recognize that using land-based observations likely biased our data toward foraging in shallower areas and led to higher estimated rates of diving compared with telemetry-based observations (see Ralls et al. 1995
). Any such bias in the model diving rate parameters would increase the projected exposure of the NKI sea otters to SSOR (because of over-estimated ITZ diving rates), rendering the modeling results more conservative. We used separate estimates of foraging parameters for each sea otter class because numerous studies confirmed our observations of significant gender- and age-specific differences in foraging (Garshelis et al. 1986
; Johnson and Garshelis 1995
; Ralls et al. 1995
; Bodkin et al. 2004
; Laidre and Jameson 2006
) and even individual variation in prey selection (Estes et al. 2003
). As will be shown, variations in rates of ITZ foraging and pit excavation are important to theecological modeling outcomes because these rates directly affect PAH exposures for the at-risk sea otters.
The consumption rate (1019 kJ·kg–1
dry weight of sea otter per day) and the energetic value of each prey species (J·mg–1
dry weight of prey) were taken from Dean et al. (2002)
. The model uses daily dietary energetics needs, prey frequency distribution, and energetic value of each prey species to calculate the required mass of each prey consumed per hour when feeding in order to sustain the sea otter. The prey concentrations of each PAH analyte are assigned to that mass, depending on the source of the prey (ITZ SSOR subdivision, ITZ non-SSOR subdivision, or STZ). Seawater consumption rate is assigned at 23% of total water influx (Costa 1982
), or 62 ml·kg–1.
converted to quantity per hour for the feeding period of each day. For all dietary, seawater, or oil-phase PAHs, the quantity ingested is assumed to be completely assimilated and contribute to the total daily dose. However, for particulate-bound PAHs in the plume of particles from the pit excavation, a PAH-specific partitioning coefficient controls the amount of each PAH released from the sediment particle and assimilated (data taken from Durrell et al
.  for each of the 41 PAHs analyzed).
The model calculates the total hourly intake for each of the 41 modeled PAH analytes from each exposure pathway for the at-risk sea otters, with or without cooccurrence with SSOR, and for the not-at-risk sea otters, aggregated to the PAH-ring class (discussed later), and summed at the end of every 24 h to a daily dose that is normalized to the weight of the sea otter (as mg PAH·kg–1 sea otter·day–1). For each class of sea otter, 1000 yr of activity were simulated, resulting in an initial distribution of daily doses of 2–3-ring PAHs, 4–6-ring PAHs, and TPAH for the paired (at-risk and not-at-risk) sea otters. These daily doses were converted to average daily doses by simulating a sequence of days that could occur in the life of a sea otter in order to derive chronic exposures for comparison with chronic TRVs. To do this, a secondary distribution was created by randomly sampling the initial distribution (with replacement) to select a 180-day sequence (representing the duration of phase of each sea otter class, except that a 120-day sequence was used for the shorter-termed class of older pups) to derive a single average daily dose. This process was repeated 500,000 times, generating the modeled population of NKI sea otters, from which the chronic dose to the 99.9% quantile, maximum-exposed sea otters was calculated.
The model was developed with a consistent bias toward conservative assumptions or data (where conservative
refers to yielding a higher estimate of risk). Some examples, among many, of this conservatism include: (1) if a simulated sea otter was assigned at random to dig a pit that intersected SSOR, then all the pits dug during that hour were assigned to intersecting SSOR; (2) TRVs are based on the geometric 95% quantile lower confidence interval of the toxicity data within each PAH-ring group; (3) the 4–6-ring PAH TRV is based on toxicity data for the most toxic analyte in the group, benzo(a)pyrene, even though the toxicities of the other PAHs are much lower (e.g.,
an order-of-magnitude lower for pyrene and fluoranthene; Di Toro et al. 2007
) and even though benzo(a)pyrene represents ≤0.3% of the composition of the modeled assimilated PAHs. On the other hand, the model balances conservatism with common sense and plausibility. For example, the excavations were limited to the tidal zone where sea otters actually forage (99% of ITZ pits occur in the LITZ; Boehm et al. 2007a
) and in habitats where their infaunal preyexist and where sea otters can physically dig pits (i.e.,
the sediments that constitute clam habitat). This contrasts with the assumptions of Short et al. (2006)
, in which probabilities unrealistically included sea otters digging pits in the MITZ and UITZ and in boulder-cobble fields.
Two types of sensitivity analyses were conducted: (1) structural-sensitivity analyses that examine six alternate model constructs based on different data sources or approaches to assessing co-occurrence; and (2) parameter-sensitivity analyses that examine eleven variations in selected base model parameters.
Two structural-sensitivity analyses directly assigned probabilities of co-occurrence of a sea otter pit with SSOR based on two estimates from Short et al. (2006)
: (a) assuming sea otters can dig pits throughout the LITZ, MITZ, and UITZ up to the +4.8 m MLLW tide level (estimated in Short et al
. , but which does not actually occur); and (b) limiting sea otter pits to the LITZ ≤+0.8 m MLLW (the more realistic case derived from the Short et al
.  data). A third structural-sensitivity analysis assigned the sea otter to foraging only at an SSOR-subdivision throughout the simulation period. The fourth structural-sensitivity analysis used the sea otter diving data for PWS derived from Ballachey and Bodkin (2006)
, instead of our observational data (except we kept our conservative assumption in the model that sea otters excavated a new pit on each dive rather than those authors’ assumption of 2–5 dives per new pit). The final two sets of structural-sensitivity analyses combined the Short et al. (2006)
parameterization (for probabilities) and the Ballachey and Bodkin (2006)
parameterization (for sea otter diving characteristics) as follows: (a) upper-bound (more conservative) parameters, with sea otter pits throughout the ITZ and two dives·pit–1
; and (b) lower-bound (more realistic) parameters, with foraging only in the LITZ and five dives·pit–1
The parameter-sensitivity analyses were conducted on three sea otter classes (older pup, adult male territorial phase, and adult female with pup), chosen to conservatively represent all classes. These sensitivity analyses used identical random-number seeds as the counterpart base models so that any difference between the results is attributable solely to the change in the tested individual parameter. Individual parameters were primarily changed in the direction to increase doses, except for those cases where the parameter was already at its most conservative level in the base model (e.g., 0% dilution for surfaced oil film, and 100% of pits in co-occurrence hour assigned to intercepting SSOR). The parameters examined were: Kow partitioning coefficient for particulate-bound PAHs (all partitioning coefficients were changed to 1.0 to simulate complete bioavailability and assimilation of sediment-bound PAHs); sea otter pit dimensions (width of all pits increased by 10 cm); efficiency of fur grooming to remove particulates (grooming efficiency halved); number of SSOR pits within a co-occurrence hour (reduced by factor of ¼); number of sea otter dives per new pit (changed from 1 to 5); dilution of the oil-phase components from the pit to the water surface where feeding occurs near high tide (oil film thickness reduced by 1/10); increased thickness of the oil-phase coating on clams and/or the sea otter paws (doubled); decreased thickness of the oil-phase coating on the sea otter paws; and fraction of sediment particles attaching to the sea otter's fur (doubled).
The model was subject to extensive quality assurance (QA) to examine model structure, equations, parameters, data sources, documentation, and each simulation output, following a USEPA-approved QA plan developed by the authors for a previous model. All models, parameters, QA runs, and simulation outputs have been archived for further analysis and reproducibility. Altogether 126 scenarios were simulated, resulting in a total modeled population of 63 million individual sea otters.