Using individual-based bioenergetic models to predict the aggregate effects of disturbance on populations: A case study with beaked whales and Navy sonar

Anthropogenic activities can lead to changes in animal behavior. Predicting population consequences of these behavioral changes requires integrating short-term individual responses into models that forecast population dynamics across multiple generations. This is especially challenging for long-lived animals, because of the different time scales involved. Beaked whales are a group of deep-diving odontocete whales that respond behaviorally when exposed to military mid-frequency active sonar (MFAS), but the effect of these nonlethal responses on beaked whale populations is unknown. Population consequences of aggregate exposure to MFAS was assessed for two beaked whale populations that are regularly present on U.S. Navy training ranges where MFAS is frequently used. Our approach integrates a wide range of data sources, including telemetry data, information on spatial variation in habitat quality, passive acoustic data on the temporal pattern of sonar use and its relationship to beaked whale foraging activity, into an individual-based model with a dynamic bioenergetic module that governs individual life history. The predicted effect of disturbance from MFAS on population abundance ranged between population extinction to a slight increase in population abundance. These effects were driven by the interaction between the temporal pattern of MFAS use, baseline movement patterns, the spatial distribution of prey, the nature of beaked whale behavioral response to MFAS and the top-down impact of whale foraging on prey abundance. Based on these findings, we provide recommendations for monitoring of marine mammal populations and highlight key uncertainties to help guide future directions for assessing population impacts of nonlethal disturbance for these and other long-lived animals.

Budget (DEB) formulation of Kooijman (2009).Demand-driven growth implies that rates of energy expenditure are independent of the energy acquired from feeding, because growth rate does not depend on energy intake.Instead, individuals regulate their body condition (fraction of reserves to total mass) by adjusting their feeding effort (see below).For fetuses, structural length is a linear function of time since conception (  ).Structural mass is directly related to structural length through a length-weight relationship and is therefore a fixed function of age too.Total body mass equals the sum of structural mass and reserve mass and includes structural mass of any fetus.Maintenance body mass is a weighted measure of body mass that discounts the mass-specific contribution of reserve mass to metabolic rate through parameter   (Table S1).
Energy acquisition from prey feeding (IR) and milk suckling (IL) are directly used to cover all energetic costs.All individuals spend energy on maintenance costs (CM) and growth (CG).Gestation (Cp) and lactation (CL) costs are paid by, respectively, pregnant and lactating females only.Females can be simultaneously pregnant and lactating.If energy acquisition is insufficient to cover all energy expenditures, additional energy is made available by catabolizing reserves.Vice versa, if energy intake exceeds energy expenditure, the surplus energy is stored as reserves (anabolism).Consequently, dynamics of reserve mass (dF/da) depend on the rates of energy intake and expenditure, which in turn depend on the reproductive status of an individual: With the reserve mass at birth as initial condition F(0) = Fb.In eq. ( 1), the anabolic   =  + if   > 0 and catabolic   =  -if   < 0 conversion parameters account for the efficiency of storing and utilizing reserves.
The function describing energy assimilation rate from prey feeding (IR in MJ day -1 ) consists of three different components.The first component is a type II functional response of prey density R, with a maximum consumption rate that depends on structural surface-area S 2/3 and the handling time parameter h.The attack rate parameter of the functional response (   ) varies between areas (j) (Table S2).The second component represents the feeding effort ) and describes the fraction of time spent foraging.The feeding effort is modelled as a decreasing sigmoidal function of body condition (F/W) that equals 0.5 at the target body condition  = 0.3.The last component of prey foraging models the prey feeding efficiency (capture success), which is lower for younger whales and gradually increases with age according to a hyperbolic function ( The lactation rate is proportional to calf structural surface-area (   2/3 ) and subject to the decrease in feeding effort of the calf in response to calf body condition.The lactation rate is modified by a hyperbolic function of maternal body condition (Fm / Wm) that equals zero if maternal body condition falls below   and one at   /  = .Furthermore, beyond age TN, milk demand by the calf declines as the calf approaches weaning age (TL).
Costs for structural growth (CG) are obtained from the derivative of structural mass with respect to age multiplied by costs of growth parameter   .Metabolic costs follow a quarter power relationship of maintenance body mass (WM) following Kleiber (1975) and Lockyer (2007).Pregnancy costs (CP) are proportional to fetal growth in structural mass.Fetal survival is assumed to be 0.8 over the whole of the gestation period and probability of abortion is independent of female state.Lactation costs (CL) for the female are proportional to the lactation rate (IL) of the calf, accounting for the efficiency of milk production and consumption through parameter (  ).

Mortality and survival
Calves and weaned females experience age-dependent mortality rate (Db(a)) according to a Siler model with three components: i) constant background mortality, ii) early life mortality that declines with age and iii) senescence mortality that increases with age (Barlow and Boveng, 1991;Siler, 1979).Weaned males experience constant background mortality only, with rate   .In addition, individuals experience additional starvation-induced mortality if their body condition falls below the starvation body condition threshold   = 0.15.The starvation mortality rate is a hyperbolic function that increases with decreasing body condition (Table S2)., which is a decreasing function of age with H(0) = 1.0.In natural populations, there is variation in age at death even in the absence of starvation mortality.To capture this variation, each individual was assigned a survival threshold: a uniform random number between 0 and 1.If H(a) falls below this threshold the individual was considered death and removed from the population.

Reproduction
Whales are assumed to initiate pregnancy when the absolute amount of reserves exceeds the 'pregnancy threshold', defined as   =    +   , where    is the reserve mass threshold below which starvation mortality occurs.Once the female becomes receptive she is assigned the 'waiting' status, as she awaits the onset of pregnancy.Each time a female enters the waiting period, its duration is randomly determined from an exponential distribution with a mean of 445 days.Gestation last TG days and during gestation fetuses experience a daily mortality rate of   , which is independent of the body condition of their mother.Hence, termination of pregnancy by the female was not considered.Fetal growth in structural length is a linear function of time since conception (Table S1).Fetuses do not contain any reserve mass, as they are directly linked to the reserves of their mother.Upon parturition, the minimum amount of reserves to offset starvation mortality (Fb) is transferred from the mother to the neonate (Table S2).Lactation lasts TL days and females can reinitiate the next pregnancy (enter the waiting period) during the last TG days of the lactation period if their reserve mass exceeds Fpreg.This restriction is to prevent females from having two calves simultaneously.Waiting females can become pregnant during the last phase of lactation if the waiting period happens to be shorter than the gestation period (TG).

Prey dynamics
Prey dynamics (dR/dt) are modelled using a semi-chemostat approach, with a fixed prey productivity term   and a first-order prey removal term .In absence of predation by whales, prey growth rate is a linear decreasing function of prey density and prey density equilibrates at  =   .During analysis, the maximum prey density   was varied to obtain a beaked whale population abundance that matched observed densities (around 100 Zc and 60 Md individuals).Total consumption of prey by all whales decreases prey biomass.For sake of simplicity, prey biomass (MJ m -3 ) is expressed in MJ assimilated energy and therefore includes any overheads and efficiency of prey assimilation by whales.The volume scalar V converts the volume unit of prey density (m 3 ) to the volume unit of whale density (V m 3 ).The ordinary differential equation of prey biomass equals: 2. Parameterization Little is known about the biology of beaked whales compared to other odontocetes and several parameter values for Globicephala melas (hereafter Gm) from Hin et al. (2021Hin et al. ( , 2019) ) are therefore retained for Zc and Md.Model parameters and their sources are summarized in Table S3.Derivation of several parameters is described below.

Asymptotic length
Length estimates for beaked whales are given in Mead (1984) and, more recently, MacLeod (2006).According to Mead (1984), the maximum length of females Zc is 754 cm.However, Heyning (1989) notes that reported lengths > 7 m for Zc probably represent misidentified individuals, as these are all derived from individuals that co-occurred with larger beaked whale species (Berardius spp.and Hyperoodon).Indeed, the vast majority (93%) of Zc individuals reported in MacLeod ( 2006) are under 7 m in length and according to Heyning (1989), the largest accurately identified Zc measured 693 cm in length.Data from catches of Zc in the North Pacific indicate a maximum length for females of 670 cm (MacLeod, 2006;Nishiwaki and Oguro, 1972).Because there is no significant different in total length between sexes of Zc, we adopt 670 cm as our asymptotic length estimate for Zc (Table S3).For Md we adopt an asymptotic length of 475 cm, based on the observation by MacLeod ( 2006) that Md rarely reaches over 4.8 m in length.

Size at birth
For Zc, we adopt the mean length at birth estimate of Mead (1984), which is 270 cm.Following (New et al., 2013), length at birth for Md was set to 198 cm.This value also falls in between the length of longest fetus (190 cm) and length of shortest calf (261 cm), as reported by Mead (1984).

Von Bertalanffy growth constant
The Von Bertalanffy growth constant was calculated by rewriting the growth equation and using estimates of mean age (am) and length at sexual maturity (lm).
For Zc, Mead (1984) reports mean length at sexual maturity of 580 cm and age at sexual maturity was estimated at 2250 days, which give k = 0.000663 for Zc.Applying the same procedure for Md with mean length and age at sexual maturity of, respectively, 378 cm and 3285 days (9 yrs; Claridge, 2013) gives k = 0.000319.

Structural lengthmass relationship
To estimate the scaling constant and exponent of the lengthmass relationship of structure, we collect length-weight data from stranded Zc and Md individuals.Data from emaciated or starved individuals were excluded from analysis, as well as those reported to be in a moderate or advanced state of decomposition.Because of the limited number of data points, we pooled observations from both species and fitted a single equation.This approach is supported by the general observation that body shape and proportions are extremely conserved in Ziphiids (Mead et al., 1982).The data and the fitted equation are shown in Figure S1.We assume that 27% of total weight consisted of reserve mass and multiplied the scaling constant by (1 -0.27) to obtain  1 = 2.93 ⋅ 10 −5 .Applying these estimates to newborns that are assumed to receive minimum amount of reserves (Table S2) leads to Fb = 13.9 for Md and Fb = 33.2 for Zc (Table S3).

Energetic parameters
We set the metabolic rate constant   to 0.6, which equals 2 times the resting metabolic rate as estimated by Kleiber (1975).This low proportionality constant reflects the 'cheap' body structure and low energy expenditure that is associated with extreme deep-diving lifestyle of beaked whales (Pabst et al., 2016).The parameters  + and  − represent the anabolic and catabolic conversion efficiencies of build-up and use of reserves.For Gm, Hin et al., (2019) argued that the catabolic conversion efficiency resembled the energy density of lipids (reported as 40 MJ kg - 1 by Lockyer (1993)) and assumed a roughly 90% efficiency of catabolic conversion.This led to the value ε − = 35.Hin et al. ( 2019) considered anabolic conversion less efficient and set ε + = 55.However, lipids in the blubber of Gm are mostly (>90%) triacylglycerols, while blubber lipids of beaked whales are dominated by wax esters (70-100%; Koopman, 2007Koopman, , 2018)).Although wax-rich blubber has a higher energy density than blubber consisting of triacylglycerols, it is currently unknown whether wax esters in the blubber of beaked whales (and other odontocetes) are used as an energy source during periods of energy deficiency (Koopman, 2018).Moreover, the blubber of beaked whales displays low density of microvasculature, suggesting that blubber of beaked whales is relatively inert and mainly serves other functions than that of an energy store (Pabst et al., 2016).If indeed beaked blubber plays a marginal role as energy store, beaked whales must mobilize other types of lipids or different biochemical compounds (such as proteins) to overcome periods of energy deficiency.To account for the possibility that beaked whales have a lower potential to store energy, we take a precautionary approach and considerably decrease the values of the anabolic and catabolic conversion efficiency parameters compared to Gm, and adopt  − = 20 and  + = 30 by default.The parameter   denote the energetics cost of growing one kg of structural mass.For beaked whales, no empirical estimates are available for this parameter, but for Gm it was derived to be around 30 MJ kg -1 by Hin et al. (2019).In line with the notion that beaked whales are build 'cheap', we deviate from this value and adopt a lower value of   = 20 MJ kg -1 .
The lactation rate scalar parameter   was recalculated for Zc and Md using the approach described in Hin et al. (2019).This approach assumes that milk supply of a female at target body condition (   =   ) is sufficient to cover metabolic and growth expenses of her calf during its first year of life.For Gm this led to a value of   = 2.7 (Hin et al., 2019) and values for Md and Zc are similar (Table S3).
The maximum consumption rate of the type II functional response of prey density was estimated based on dive statistics of Zc.The average foraging dive cycle (time between start of consecutive deep foraging dives) of Zc is on average 121 min (Tyack et al., 2006), which allows Zc to perform 11.9 deep dives per day.On each foraging dive there are around 30 prey capture attempts (Tyack et al., 2006), resulting in a maximum consumption of 357 prey items if all attempts are successful.A calorific content of 850 kCal per prey (3.56 MJ; Southall et al., 2019) leads to a maximum consumption rate of 1270 MJ per day.Using this value, the handling time parameter h of the type II functional response of a fully-grown Zc individual with structural mass of  ∞ =  1  ∞  2 = 2398 kg, becomes h = 2398 2/3 / 1270 = 0.141.We assume that there are no differences in maximum consumption between Zc and Md other than those related to size.In other words, we do not have any reason to believe that Zc and Md individuals of identical size have a different maximum consumption rate.With the derived value of h = 0.141, the maximum consumption rate of a fully-grown Md becomes

Pregnancy threshold
The female age at first reproduction is largely determined by the age at which female reserve mass exceeds the pregnancy threshold for the first time ('age at first receptive').In addition to the amount of reserves below which individuals experience starvation mortality (  ), the pregnancy threshold includes a constant value (Fneo) that can be used to control the modelled age at first receptive.For Gm, Hin et al. (2019) assumed that Fneo consisted of the amount of stored reserves that are needed to produce a neonate, which can be derived from the structural size at birth (  =  1    2 ) and amounts to: Here, the first term represents the amount of female reserve mass required to build the structural mass of a newborn and the second term equals the amount of reserves transferred from the mother to the newborn at birth (Fb).For Gm this approach resulted in a modelled age at first reproduction that corresponded well with its observed counterpart (Hin et al., 2021).To infer values of Fneo that result in good correspondence between observed age at maturity and the modelled age at first receptive we ran the full population model for a range of Fneovalues.We should note that the observed age at maturity for Md (9 y) was assumed one year earlier than the observed age at first parturition (Claridge, 2013) and does not include a receptive period before the onset of pregnancy, as is included in our model.Accounting for this 'waiting' period, the modelled age at first receptive should correspond to 9 -445 / 365 = 7.8 y.For Zc we do not know how the estimate for the age at maturity (6.2 y) was derived and the corresponding modelled age at first receptive should therefore be in the range of 5 -6.2 yrs.The values of Fneo reported in Table S3 lead to good correspondence between modelled and observed age at first maturity.Cuvier's beaked whales (Ziphius cavirostris) at SOAR that were used to estimate transition rate between discrete geographical areas used in the model.Red lines demarcate the boundaries that were used to allocate beaked whale tracks to different model areas.The parts of the tracks that overlap with gray shadings are considered to be on the hydrophone range and all other parts of the tracks (no shading) are off-range areas.In the right panel, these boundaries coincide with the hydrophone array at the Southern California Anti-submarine Warfare Range (SOAR).Bathymetry data were obtained from the GEBCO_2022 Grid (GEBCO Compilation Group, 2022).

Figure S1 :
Figure S1: Length-weight relationship of various stranded beaked whales

Figure S2 :
Figure S2: Tracks of a) seven Blainville's beaked whales (Mesoplodon densirostris) at AUTEC and b) 12Cuvier's beaked whales (Ziphius cavirostris) at SOAR that were used to estimate transition rate between discrete geographical areas used in the model.Red lines demarcate the boundaries that were used to allocate beaked whale tracks to different model areas.The parts of the tracks that overlap with gray shadings are considered to be on the hydrophone range and all other parts of the tracks (no shading) are off-range areas.In the right panel, these boundaries coincide with the hydrophone array at the Southern California Anti-submarine Warfare Range (SOAR).Bathymetry data were obtained from the GEBCO_2022 Grid (GEBCO Compilation Group, 2022).

Figure S3 :
Figure S3: Smooth terms of the GAM on the hydrophone data from AUTEC

Figure S4 :
Figure S4: Smooth terms of the GAM on the hydrophone data from SOAR

Figure S6 :Figure S7 :
Figure S6: Distribution of MFAS use aggregated by day (daily sum of standardized sonar area) for both western and eastern areas of the AUTEC and SOAR ranges

Figure S8 :
Figure S8: Lifetime reproductive output in terms of the mean number of calves born and weaned during the entire life of a single female.Bars represent mean values from all females living in stationary populations under a certain disturbance regime.

Figure S9 :
Figure S9: Distribution of female body condition as a function of reproductive status with and without disturbance from MFAS.Triangles indicate mean values and shapes indicate the distribution and total range of the data.Females with body condition below the starvation body condition threshold of 0.15 suffered from increased mortality.The lactating female group also represent females that are simultaneously 'waiting' or 'pregnant'.Data represent body condition values of all females in a stationary population from single simulation of 100 y, with output collected each day.For Zc a Rmaxvalue of 2 was used for the cessation of displacement scenario, as the population went extinct with the default value of Rmax = 1.35.

Table S2 :
Equations of the energy budget model for an individual whale

Table S3 :
Model parameters

Table S4 :
Results of the Generalized Additive Model (GAM) on the hydrophone data from AUTEC.

Table S5 :
Results of the Generalized Additive Model (GAM) on the hydrophone data from SOAR.

Table S6 :
Mean (sd)of several life history variables and body condition (BC) of modelled Md and Zc females.
* Disturbance leads to cessation of foraging