Estimating and decomposing total uncertainty for survey-based abundance estimates of Norwegian spring-spawning herring

Abstract

Løland, A., Aldrin, M., Ona, E., Hjellvik, V., and Holst, J. C. 2007. Estimating and decomposing total uncertainty for survey‐based abundance estimates of Norwegian spring‐spawning herring. – ICES Journal of Marine Science, 64: –.

A method is developed to quantify the uncertainty of abundance and age-specific abundance estimates of Norwegian spring-spawning herring, accounting for the combined effect of several error sources. We present an abundance estimate adjusted for vessel avoidance, acoustic shadowing, and depth-dependent target strength, and highlight the potentially dramatic combined effect of these adjustments. Uncertainty in this estimate is affected by the precision of the adjustment formulae as well as by the sampling errors in acoustic and trawl measurements. This total uncertainty is decomposed into contributions from each source. Here, the method is applied to data from surveys on the overwintering stock in the Vestfjord system and off Vesteraalen in northern Norway in November/December of the years 2001–2004.

acoustic abundance estimation, total estimation uncertainty, uncertainty decomposition

Introduction

The Norwegian spring-spawning herring (Clupea harengus) stock has been, at times, the largest fish stock in the Northeast Atlantic [6.8 (8.9, 10.3) million tonnes in 1994 (1995, 1996); Dommasnes et al., 2004]. The adults spawn along the Norwegian coast, feed in the Norwegian Sea and adjacent areas, and overwinter in various parts of the Northeast Atlantic.

On a regular basis, the Norwegian Institute of Marine Research (IMR) surveys the stock to estimate age-specific abundance, a task vital to the assessment and management of the stock. A winter survey is carried out in the Vestfjord system and off Vesteraalen in November/December each year. The spatial density of fish is measured by echosounder, and age/length data are gathered by pelagic trawling. These data are combined into estimates of abundance-at-age.

The standard total abundance estimate (Thoresen et al., 1998) is based on the assumption that the acoustic signal is proportional to fish density when the backscattering is corrected for fish size. Further, abundance-at-age is computed using the observed age proportions from trawl samples. However, because of factors such as the acoustic shadowing effect (Zhao and Ona, 2003), depth-dependent target strength (TS) (Ona, 2003), and vessel avoidance (Vabø et al., 2002), this method of abundance estimation may be seriously biased (Aglen, 1994). IMR has for some years used an improved estimator to correct for the acoustic shadowing effect caused by fish schools being dense and inducing attenuated signals deeper down. Further improvements, accounting for depth-dependent TS and vessel avoidance, are presented here. This improved estimator is a preliminary version, and we demonstrate the effects of bias correction.

Traditionally, the uncertainty in acoustic estimates of abundance has either not been or only partly quantified (Rivoirard et al., 2000; Rose et al., 2000; Tjelmeland, 2002), despite the increasing recognition of its importance in fish stock management (Aglen, 1994; Frederick and Peterman, 1995; Potter et al., 2003). The focus of this paper is to present a method for quantifying the uncertainty of the total and age-specific acoustic estimates of abundance, accounting for the combined effect of systematic and random sampling variability.

We chose five important sources of uncertainty: trawl sampling, acoustic sampling, and the three acoustic signal corrections mentioned earlier. The uncertainty from echosounder calibration and species identification is not included, because they are assumed to be negligible for this particular survey. For other species and surveys, the identification of species may be an important source of bias and uncertainty (Tesler, 1989; Aglen, 1994; Demer, 2004; O’Driscoll, 2004). The echosounder calibrations were conducted inside Ofotfjord (in the Vestfjord system) before or during each survey, according to standardized procedures (Foote et al., 1987). During the surveys, herring are easily identified from the raw signal. The potential mixture of predators, mainly saithe (Pollachius virens) and cod (Gadus morhua), is estimated to be < 1% (roughly estimated from acoustic survey and catch data) of herring abundance, whereas the herring abundance estimate was some six million tonnes. Similar quantities of blue whiting (Micromesistius poutassou) are found under the herring layers. The water masses in which the herring overwinter are quite similar for the deep fjords and outside Vesteraalen, with temperatures between 7°C and 8°C and salinity between 34.5 and 35.1 (Dommasnes et al., 1994).

The uncertainty is initially quantified for each source separately. Acoustic and age/length sampling uncertainties are estimated by non-parametric bootstrapping. Acoustic signal corrections are described by parametric functions, assumed to have an error-free shape, but uncertain parameters (Foote, 1999; Ona, 2003; Appendix). The parameter uncertainty is represented by simulated parameter values from the estimated distributions. Finally, the total estimate of uncertainty is found by combining the non-parametric and parametric simulations into the total and age-specific abundance components. In addition, we deduce the proportion of uncertainty contributed by each source.

Methods

Surveys

The data consist of acoustics and age/length data from the Norwegian spring-spawning herring November/December surveys from 2001 to 2004. An overview of the surveys is given in Table 1. In 2001 and 2002, Vestfjord was the only system surveyed, because it was assumed that almost all adults overwintered there. However, in 2002, commercial fish vessels observed herring off Vesteraalen too (Figure 1), and the analysis performed in the Results section indicates that some herring may have been wintering off Vesteraalen already in 2002. Therefore, in 2003, equal sampling effort was put into surveying the concentrations outside Vesteraalen. In 2004, most adult herring wintered outside the Vestfjord system, so the survey concentrated its effort accordingly.

Figure 1.

Map of Norway and the survey areas in the Vestfjord system and off Vesteraalen.

Table 1.

Survey overview, 2001–2004.

Area	Date	Number of EDSUs	Number of trawl stations	Number of measurements with age and length	Number of measurements with length and non-readable age	Number of measurements with length only
Vestfjord system	21–29 November 2001	1 144	9	818	81	0
Vestfjord system	30 November—3 December 2002	350	10	907	94	0
Vestfjord system	30 November—4 December 2003	346	13	1 170	130	0
Vesteraalen	8–18 December 2003	1 000	10	781	219	0
Vestfjord system	2–3 December 2004	262	9	793	107	900
Vesteraalen	5–18 December 2004	1 555	7	677	23	601

EDSUs: number of elementary distance sampling units (1 nautical mile).

The meteorological conditions are quite different for the survey areas investigated. For the years when herring overwintered inside the Vestfjord system, the wind speed and wave height during surveying seldom exceeded 10 m s^–1 and 1–2 m, respectively, because most of the fjord system is quite sheltered. In the wintering area outside Vesteraalen, the survey was conducted under variable, but often rough sea conditions, with average wind speeds between 20 and 35 m s^–1 and corresponding wave heights of 3–7 m. In rough weather, bubble attenuation may give a significant bias in acoustic measurements (Dalen and Løvik, 1981). However, the vessel used here, the RV “Johan Hjort”, has a drop keel on which all the echosounder transducers are mounted (Ona and Traynor, 1990), so it has an excellent bad-weather acoustic performance. The drop keel is usually fixed 1 m below the vessel keel in good weather, and 3 m down when surveying in rough weather. The calibration data are valid for both positions, and transducer draft settings refer all echograms to sea surface.

During all surveys, a Simrad EK500 operating at 18, 38, and 120 kHz was employed, with all frequencies calibrated, usually in Ofotfjord before the surveys, but only the data from the 38 kHz system were used in the abundance estimation. The principal acoustic signal is the area-backscattering coefficient, s_A (m² Nm⁻²). All symbols are defined in Table 2. The result is stored in a database in depth channels 10 m thick, averaged over 0.1 nautical mile (185.2 m) along the survey transects, except for the Vestfjord system in 2002 and Vesteraalen in 2004, where data were available on a scale of 1 nautical mile only. For each measurement, the position x and time t is also recorded. Figure 2 shows the 2004 survey off Vesteraalen. As the exact location of the herring is not known before the survey, a large portion of the survey area in Figure 2 has no herring registrations, so there are many zero observations.

Figure 2.

Overview of the survey off Vesteraalen in 2004. The line is the survey line. The backscattering from herring is here integrated over the water column and averaged over each 1 nautical mile. The circles have diameter proportional to the acoustic signal energy log (s_A). Lines with no circle correspond to a zero measurement. The triangles represent trawl hauls.

Table 2.

Symbols and units.

Symbol	Unit	Description
N	None	Number of fish
Na	None	Number of fish at age a
pa	None	Proportion of fish at age a
sA	m2 Nm−2	Area-backscattering coefficient
sA(d)	m2 Nm−2	Area-backscattering coefficient at depth d
s̃A(d)	m2 Nm−2	Area-backscattering coefficient at depth d, adjusted for vessel avoidance
$\stackrel{\approx }{s}$A (d)	m2 Nm−2	Area-backscattering coefficient at depth d, adjusted for vessel avoidance and the shadowing effect
A	Nm−2	Area
ρi	Nm−2	Average density of fish per area unit in area A at location i
ρi(d)	Nm−2	Average density of fish per area unit in area A at location i and depth d
l	cm	Fish length
σ(l)	m2	Back-scattering cross section of fish of length l
σ(d, l)	m2	Back-scattering cross section of fish of length l, per depth d, derived by Ona (2003)
σtotal	None	Standard deviation of the abundance estimate
x̂	–	Estimate of quantity x
x*	–	Bootstrap replicate of quantity x

Trawl hauls were taken at irregular intervals. Usually, 100 fish were chosen randomly from each haul to measure their length (and weight). When possible, age was determined from fish scales. However, for ∼10% of the fish, mainly older ones, age was not determined. For some trawl hauls, the length of a further 100 herring was measured, but the age of these fish was not read. Table 1 shows the number of trawl hauls during each survey, and the locations of such sampling are displayed in Figure 2.

Overall, few fish were older than 6 years. Moreover, the trawl hauls tended to be quite heterogeneous, some hauls consisting of young and others mainly of older herring, probably because fish of similar size/age school together.

Standard abundance estimators

Let ρ be the average fish density per unit area in an area of total size A. The total number of fish is then N = ρ · A. Further, let p_a be the true proportion of fish at age a, and the corresponding number of fish N_a = N · p_a.

The acoustic estimates of N and N_a, defined by Thoresen et al. (1998) and Rivoirard et al. (2000), are based on the n acoustic signal measurements s_Ai, i = 1, … , n, from various locations (and times) along the transects, integrated over the water column. The standard TS relationship for herring (Foote, 1987), TS = 20 logl−71.9 db, corresponds to the area backscattering strength from one fish of length l, given by σ(l) = 8.1135 · 10⁻⁷l². The fish density in location i can be estimated from

1

where l² is the sample mean of the squared fish lengths from all relevant trawl hauls. This leads to the standard estimate of total abundance:

2

where ρ = (1/n)Σ_iρ̂_i^standard is the sample mean of the estimated densities.

The probability of missing age measurements varies systematically with age. Therefore, the sample mean fish age is a biased estimator of the true proportion at age p_a. Instead, it is estimated by the less biased

where p̂_l is the sample proportion at length group l and p̂_a|l is the sample proportion at age a, given length group l. Each length group is 0.5 cm. The abundance-at-age is then given by

3

This calculation is performed per subarea A. Summing over subareas gives estimates for the total area. In the Vestfjord system, the subareas are three separate smaller fjords (outer Vestfjord, Ofotfjord, and Tysfjord). In Vesteraalen (Figure 2), the subareas are rectangles of 1° latitude by 30′ longitude.

Improved abundance estimators

Here, we describe the improved estimators of abundance, correcting for vessel avoidance and shadowing. They incorporate an improved TS relationship to convert the (corrected) signal to an estimate of fish abundance. The three correction factors are based on parametric functions, with parameters estimated from experimental (non-survey) data. The values of the parameters involved and their sources are given in Table 3.

Table 3.

Parameter value and sources.

Parameter	Estimate	Standard deviation	Source
λ	2.41	0.33	Foote (1999)
σ0	2.98	0.19	Ona (2003)
γ	−0.23	0.035	Ona (2003)
a(55)	2.626159	–	–
a(65)	2.463773	–	–
a(75)	2.450203	–	–
a(85)	2.215783	–	–
a(95)	1.704268	–	–
a(112)	1.194746	–	–
a(137)	0.5915432	–	–
a(162)	0.5962801	–	–
a(187)	0.4893917	–	–
a(212)	0.4174365	–	–
a(237)	0.3775844	–	–
a(262)	0.3595151	–	–
a(287)	–0.3884178	–	–
a(312)	–0.4704613	–	–
a(337)	–0.1545417	–	–
a(362)	–0.2495875	–	–
a(387)	–0.009496362	–	–
b	0.9215633	–	–
c	1.314305	–	–

Vessel avoidance

The presence of a vessel may influence the behaviour of fish. Vabø et al. (2002) investigated this effect on overwintering Norwegian spring-spawning herring. They interpreted systematic changes in acoustically estimated density as vessel avoidance, probably caused by a combination of changes in herring orientation and horizontal movement. The avoidance reaction was particularly strong for the shallowest fish, and it gradually decreased down to ∼100 m. On the basis of the data of Vabø et al. (2002), we corrected the raw signal, s_A(d), measured at a transducer passing at depth d, to the undisturbed signal, s̃_A(d):

4

where c(d) and b are parameters. The correction depends both on depth and signal strength, i.e. fish density.

Vessel avoidance was estimated by fitting a linear model to the data described by Vabø et al. (2002), noting that these data are from experiments in Tysfjord in the Vestfjord system. Data from all three groups (night I, day II, and night III) were used. The model is

5

where d is depth, a(d) the intercept with the y-axis for depth d, s_A,i^ref(d) the s_A at depth d for experiment i averaged over a 70-s interval ending 88 s before the time of passage, s_A,i^pass(d) is the corresponding s_A averaged over a 7-s period at the time of passage, and ε_d,i are independent, identically distributed errors. The data and the fitted regression lines are shown in Figure 3. To find the predicted value ŝ_A^ref(d) corresponding to an observed value s_A^pass (d) at a given depth, Equation (5) is backtransformed as

6

The term c is equal to the sample mean of all exp(ε_i(d)), included to correct for bias since, generally, if y_i = log(x_i), then exp(y) ≠ x (Duan, 1983). Table 3 lists the parameter estimates from Equations (5) and (6).

Figure 3.

Data from Vabø et al. (2002) and fitted regression lines from Equation (5).

The data in Vabø et al. (2002) were stratified on depth intervals of 10 or 25 m. We needed depth intervals of 10 m. For depths below 100 m, we used a(112) for the 100–110 m and 110–120 m intervals, a(137) for the 120–130 m, 130–140 m, and 140–150 m intervals, and so on. In addition, we needed parameter values for depths above the highest interval (50–60 m) and below the lowest (375–400 m), so for those depths we extrapolated (used the estimated parameter values of the highest and lowest intervals, respectively).

Shadowing effect

Acoustic shadowing or extinction, which may appear in dense fish schools, induces attenuated signals deeper in the water column. The shadowing increases with density and vertical extension of the layer. The signal-correction formula for acoustic shadowing (first-order scattering), developed by Zhao and Ona (2003), is

7

where s̃_A(d) is the signal corrected for vessel avoidance, s̃_A(1:(d−1)) the cumulative apparent area-backscattering coefficient of the fish above depth layer d, and s̃_A(1:d) includes depth layer d. The constant K = 2/1852², and λ is a parameter, estimated with mean and standard deviation by Foote (1999). Table 3 lists the parameter values.

Depth-dependent TS: from signal to stock abundance

Ona (2003) derived a new depth-dependent TS relationship for herring:

Using the relationship between TS and back-scattering cross section, TS = 10 log₁₀(σ/4π), the area backscattering strength from "one" fish of length l at depth d is given by

8

where σ₀ and γ are parameters (and the TS is normalized to a common fish length of 32 cm). Table 3 again lists the estimated values. Ona (2003) does not present the estimated correlation between σ₀ and γ, so we used bootstrap samples obtained from the author.

Based on the correction formulae above, the fish density in location i and at depth d was estimated by

9

The improved density estimate in location i is obtained by summing over all depths by ρ̂_i = ∑_dρ̂_i(d). The abundance estimates are then given by Equations (2) and (3), as before.

Based on three independent (experimental) datasets, each set of parameters, c(d) and b (vessel avoidance), λ (shadowing), and γ and σ₀ (TS), was estimated by regression, along with parameter uncertainty. Parameter uncertainty is independent between, but dependent within, the three corrections.

Here, we corrected for vessel avoidance before correcting for shadowing, because the parameters in Equation (4) were estimated from raw data without a shadowing correction.

Measures of total uncertainty

We describe the uncertainty of the estimates of N and N_a by their estimated standard deviations σ_total and σ_a. The coefficient of variation for N is then CV_total = σ_total/N̂. There is no unique summary of uncertainty over all ages, so we defined the measure as

10

If the estimates of abundance-at-age were independent (which they are not), σ_age = σ_total. We decomposed the total variance σ_total² into the contributions from five sources: three correction factors with uncertain parameter values; acoustic sampling uncertainty, attributable to covering the area only partially; age/length sampling uncertainty, because the lengths and ages are measured only for a small number of fish and at a few locations. If the uncertainty of each source was independent and contributed additively to the total uncertainty, the total uncertainty would be the sum of the five contributions σ_k² (e.g. Var(X+Y) = Var(X) + Var(Y) if X and Y are two independent random variables). The relative contribution from source k would be

11

Although the various sources are independent, the individual uncertainties do not add up to the total uncertainty, because they are combined non-linearly (e.g. Var(X/Y) ≠ Var(X)+Var(Y), even if X and Y are independent). Therefore, the contribution to the total uncertainty from each source is not uniquely defined. Hence, we defined the contribution from the kth source as

12

where σ_total² (−k) is computed without it. The relative contribution from source k is then given by Equation (11).

Estimating total uncertainty by simulation

The five sources of uncertainty are independent, because they are based on different and independent data. Age/length and acoustic sampling can in principle be dependent if there is a spatial relationship between them. In the herring case studied here, there are large differences between the age/length distributions from the various trawl hauls, but no apparent spatial trends, so we assumed that the age/length distribution was homogenous over the entire study area and independent of the acoustic sampling. The length and age measurements are obviously dependent, so they were treated simultaneously.

We quantified the uncertainties by the following resampling method. Let θ denote all parameters in the three correction formulae [Equations (4), (7), and (8)]:

• Simulate a set of parameters θ* by sampling from their multivariate normal distributions (bootstrap distribution for the TS parameters).

• Draw a two-step bootstrap sample of age/length data, taking into account the lack of homogeneity between trawl hauls (see details below). Calculate the average of the squared bootstrapped length l²* and the estimated proportion-at-age p̂_a*.

• For each acoustic signal measurement i, compute the corresponding fish density ρ̂_i* by Equation (9) and sum over all depths, conditional on the simulated parameters and the bootstrapped length data.

• Draw a bootstrap sample ρ̂_i**, i = 1, … , n, from the densities ρ̂_i*, i = 1, … , n, calculated in step (iii), taking into account correlations in space (see details below). Calculate the average ρ**.

• Calculate the corresponding estimates N̂** of total abundance and N̂_a** of abundance-at-age [Equations (2) and (3)].

• Repeat the procedure B = 1000 times, and calculate standard deviations from the bootstrap distributions of total abundance and abundance-at-age.

In step (ii), the two-step bootstrap of the trawl samples is carried out as follows. First, M trawl hauls are drawn with replacement from the M trawl hauls. Assume that trawl haul m, m = 1, 2, … , M, consists of n_m^{age read} samples for which age is read, n_m^{non-readable age} samples for which age is impossible to read, and n_m^{length only} samples for which age is not read, but length is measured, which adds up to the total number of fish n_m in trawl haul m. Separate bootstrap samples (age/length) were drawn with replacement from each of the three data types. As all fish were measured for length, the bootstrap sample for length-measured fish only was drawn from all n_m observations.

At the end of step (iii), we have a simulated one-dimensional acoustic transect series ρ̂*_i, i = 1, … , n. In step (iv), we bootstrapped samples ρ̂_i** from ρ_i* by the block bootstrap for autocorrelated series (Carlstein, 1986; Künsch, 1989). The block bootstrap requires that the transect series is stationary. Therefore, we first divided the transect series into non-overlapping approximate stationary strata. Each stratum had at least 20 measurements (20 nautical miles). For the Vesteraalen area, the stratum limits were chosen automatically, as the limits of the subareas or rectangles in Figure 2. The stationarity assumption (constant mean and variance) within strata holds approximately, because they were relatively small. For the Vestfjord system, where the pattern of transect lines is more complex, we selected the stratum limits manually based on visual inspection. Within each stratum, we performed a block bootstrap. This consisted of dividing the transect series into blocks of consecutive observations, subsequently sampling a suitable number of blocks (with replacement). The correlation structure was preserved within each block, but not between blocks. If the block lengths were large enough, the break between blocks was less important. We used a block length of 5 nautical miles as the standard error of the sample mean stabilized for block lengths of five of more. With a block length of 5, there were at least four blocks within each stratum, because we required at least 20 measurements (20 nautical miles) within each stratum. We used overlapping blocks, and to avoid undesired edge effects, the data were wrapped around a circle to form additional blocks—a procedure known as circular bootstrapping (Lahiri, 2003).

Results

Effect of bias corrections

Abundance N was estimated using Equation (2) and the three corrections for vessel avoidance, shadowing, and depth-dependent TS (Table 4). For comparison, we include estimates without any correction, with only shadow correction and with both shadowing and depth-dependent TS corrections. As expected, the shadow and vessel avoidance corrections adjusted the estimates upwards, and the depth-dependent TS downwards. The cumulative effect was positive for some surveys (the Vestfjord system 2001, Vesteraalen 2003 and 2004) and negative for others (the Vestfjord system 2002, 2003, and 2004). The adjustment upwards was particularly apparent off Vesteraalen, but this was mainly due to the vertical distribution of fish during the survey. The effect of vessel avoidance correction was particularly strong in the first 100 m (Figure 4).

Figure 4.

Cumulative effect of corrections, relative to the original, uncorrected signal, Vesteraalen 2004. (Upper panel) The original signal is normalized to 100% on average over the water column. The numbers in parentheses are total cumulative effects. (Lower three panels) The original signal is normalized to 100% at each depth.

Table 4.

Estimated number of herring (10⁹), from Equation (2), 2001–2004.

Area and year	No correction	Shadow correction	Shadow correction and depth-dependent TS	Vessel avoidance, shadow correction and depth-dependent TS
Vestfjord system 2001	14.2	17.2	7.5	14.3
Vestfjord system 2002	5.9	7.1	3.3	3.6
Vestfjord system 2003	5.0	5.7	2.6	3.4
Vesteraalen 2003	28.0	30.0	11.7	54.7
Total 2003	33.0	35.7	14.3	58.1
Vestfjord system 2004	3.3	3.7	1.7	2.5
Vesteraalen 2004	21.3	21.9	8.7	46.4
Total 2004	24.6	25.6	10.4	48.9

Total uncertainty for estimates of abundance and abundance-at-age

The estimated coefficients of variation (CVs) for estimates of total abundance from each survey are summarized in Table 5. The CV estimated by the block bootstrap method is quite stable at 16.5–20%, except for Vesteraalen 2004, which had a CV of 33%. The high CV in 2004 was mostly due to the patchy distribution of herring that year, resulting in proximate very high and zero measurements. Of the 21 subareas with non-zero measurements (Figure 2), two contained about two-thirds of the total number of fish, and most of the total variance. Moreover, there were more fish in the upper layers in 2004, and hence greater uncertainty because of the vessel-avoidance correction.

Table 5.

Estimated CVs for total abundance estimates of spring-spawning herring.

Area and year	CV (%) based on block bootstrap	CV based on variogram
Vestfjord system 2001	16.5	13.0
Vestfjord system 2002	18.0	15.5
Vestfjord system 2003	18.5	18.5
Vesteraalen 2003	20.0	17.0
Total 2003	19.0	16.5
Vestfjord system 2004	19.5	19.0
Vesteraalen 2004	33.0	19.0
Total 2004	31.5	18.5

For comparison with the proposed approach, we also calculated the acoustic sampling uncertainty, following the recommendations of Rivoirard et al. (2000), using variogram modelling (see Appendix for details). The CVs calculated by the variogram approach for the acoustic sampling uncertainty were usually in accordance with the results from the block bootstrap approach, although usually slightly lower. For Vesteraalen 2004, however, the difference was large (33% vs. 19%). Closer inspection shows that the variogram method clearly underestimated the uncertainty in the two previously mentioned subareas. The variogram model does not fit the data for these subareas well, because they contain both extremely high and zero values.

Components of total uncertainty

Table 6 shows the contribution of each uncertainty source to estimates of total abundance. Generally, shadowing contributed very little to the total uncertainty (0–3%). The contribution from depth-dependent TS was higher (1.5–12.5%), but was low when the fish resided mainly in the upper depth layers. The parameter uncertainty from the vessel avoidance correction was high (16.5–41%). The contribution from age/length sampling was quite small, but not negligible, 4% in Vesteraalen 2004, stemming from the heterogeneous length distributions that year. Finally, acoustic sampling uncertainty was the largest contributor (46–78%).

Table 6.

Relative contribution (%) of different factors to the uncertainty of the total abundance estimate and estimates by age (in parentheses).

Area and year	Shadow correction	Depth-dependent TS correction	Vessel avoidance correction	Age/length sampling	Acoustic sampling
Vestfjord system 2001	0.0 (0.0)	12.5 (6.5)	31.5 (17.0)	2.0 (52.5)	54.0 (24.0)
Vestfjord system 2002	3.0 (1.0)	8.0 (4.0)	28.0 (17.5)	1.5 (45.5)	59.5 (32.0)
Vestfjord system 2003	1.0 (0.5)	11.0 (8.0)	41.0 (30.5)	1.0 (28.5)	46.0 (32.5)
Vesteraalen 2003	0.0 (0.0)	5.0 (4.0)	16.5 (14.0)	0.5 (17.0)	78.0 (65.0)
Vestfjord system 2004	0.0 (0.0)	8.0 (5.5)	28.5 (22.0)	0.0 (27.5)	63.5 (45.0)
Vesteraalen 2004	0.0 (0.0)	1.5 (1.0)	21.0 (13.0)	4.0 (46.0)	73.5 (40.0)

Similarly, the age-specific uncertainty σ_age can be decomposed by source (Table 6). The general picture is the same as for N, but the age/length uncertainty is now more dominant (17–52.5%), so it is comparable with the acoustic sampling component. However, if we focus on one particular age class, the uncertainty may be greater or less. Comparing Vesteraalen 2003 and 2004, the uncertainty contributions from age/length were 17% and 46%, respectively. Whereas the trawl hauls in Vesteraalen 2004 were heterogeneous, those in Vesteraalen 2003 were particularly homogenous (not shown).

Discussion

The study has demonstrated the combined effect of correcting for vessel avoidance, shadowing, and depth-dependent TS on the estimate of adult Norwegian spring-spawning herring abundance in the overwintering areas. Applied to data from the November/December survey from 2001 to 2004, we have demonstrated how the bias corrections significantly change the estimates of abundance. Furthermore, we have demonstrated how the total uncertainty is affected by parameter uncertainty from vessel avoidance, shadowing, and depth-dependent TS corrections, and uncertainty from trawl hauls and acoustic sampling. Note that we did not include modelling error for the acoustic signal corrections.

The shadowing correction in formulation (7) is actually only valid for 0.1 nautical mile, because its parameters are estimated on 0.1 nautical mile data (Foote, 1999). Part of our data series is only available at a resolution of 1 nautical mile. However, we believe that the effect on the total survey result is negligible. If there were to be an effect of using 1 nautical mile instead of 0.1 nautical mile, it would increase the abundance estimate slightly.

A major advantage of our approach is its simplicity. The many zero measurements pose no problems, as long as the transect series is stationary within each stratum and our choice of block length is appropriate. Furthermore, there is no need for data transformation. Handling the two spatial dimensions (and the time dimension) by one dimension alone is of course a drawback. The stationarity assumption is more valid in some subareas than in others. The assumption of independence between subareas is more appropriate with longer strata, but longer strata lead to less stationarity.

Many of these error sources in fish abundance estimation are discussed by Aglen (1994), but not quantified, at least not for a case study. Rose et al. (2000) and O’Driscoll (2004) apply a Monte Carlo simulation method to estimate and diagnose the sources of uncertainty in acoustic-survey estimates of fish density. Demer (2004) utilizes Monte Carlo simulation to assess the combined sampling and measurement error for an acoustic estimate of krill abundance. Previous work on Norwegian spring-spawning herring includes that of Høst et al. (2002), who describe a model-based approach to quantifying the uncertainty of the abundance estimate. The CVs are all very high compared with the CV of 4% reported by Høst et al. (2002), who used data from the December 1996 survey. However, those authors only considered age/length and acoustic sampling uncertainty as sources of uncertainty, they had more acoustic data, and the data were probably more evenly distributed than in the surveys between 2001 and 2004. Gimona and Fernandes (2003) consider North Sea herring and study two sources of uncertainty: the spatial distributions of length and the nautical area scattering coefficients. Their approach is probably biased through normal-score data transformation. Our approach is more non-parametric in nature, and takes more sources of uncertainty into account.

After correcting for vessel avoidance, shadowing, and depth-dependent TS, we assume that our improved estimator is unbiased. However, there may be further biasing factors, such as the stock boundary definition and diurnal variation. Diurnal variation (Huse and Korneliussen, 2000; Orlowski, 2005) is due inter alia to vertical fish migration from lower layers during daylight to upper layers at night, such that s_A is lower at night. Diurnal variability has already been reported both for the TS and extinction coefficients (Ona, 2006; Utne and Ona, 2006). Vabø et al. (2002) reported greater vessel avoidance by night. Implicitly, the vessel-avoidance correction should reduce the diurnal effects. However, there is good reason to believe that the corresponding formula is not fully developed. For instance, the experiments of Vabø et al. (2002) were all conducted during 2 d in Tysfjord in the Vestfjord system in 1996. In a more recent experiment in Ofotfjord (2004), the observed avoidance was much smaller than reported by Vabø et al. (2002). Similar experiments off Vesteraalen, with possibly greater background noise from waves, will be carried out in the future, which may result in a lower correction for avoidance.

The most biasing factor is inadequate coverage of the distribution area of the fish. More specifically, there may be dramatic changes in the wintering area, as found (at least) in 2002, when almost the entire survey was conducted in the Vestfjord system, but some herring were registered by the commercial fleet outside Vesteraalen. In addition, migration during the survey might induce a bias that is difficult to estimate and model. However, during winter one would expect slow movement, with swimming at minimum energy costs. This has been confirmed by moored observatory systems (Godø et al., 2005).

The approach taken here may be improved in various ways. In addition to vessel-avoidance and diurnal effects, modelling the acoustic sampling uncertainty is an area for further research. We assume stationarity within and independence between strata. These assumptions are not always satisfied in practice, and the consequences thereof are unclear. The large number of zeros, combined with extreme observations and the concurrent sampling in time and space, are challenging to model. An answer might be a space-time model, in which the naive plug-in estimate is replaced by an integral over time and space.

Nonetheless, the framework is believed to be a promising tool for uncovering the largest contributors to the uncertainty in abundance estimates. This may in turn be used to reduce the uncertainty. The present results imply that of our five sources of uncertainty, the focus should be on vessel avoidance and acoustic sampling, if reducing the uncertainty in total abundance is most important. If abundance-at-age is equally important, the number of trawl hauls should be increased, at least if there are heterogeneous trawl hauls. The number of trawl hauls needed is an area of further research.

Acknowledgements

This research was funded by the Institute of Marine Research and the Research Council of Norway, research grants 143249/140 and 154079/420. We thank Rolf Korneliussen and Kaare Hansen for great help with the data, and Dag Tjøstheim, Alf Harbitz, David Hirst, and Ingrid Hobæk Haff for helpful discussions and suggestions.

Appendix

Quantifying acoustic sampling uncertainty by variogram modelling

For comparison, we calculate the acoustic sampling uncertainty using variogram modelling, as Rivoirard et al. (2000).

Let ρ(x) denote the true density in a general location with spatial coordinate x, and assume that ρ(x) follows an intrinsic geostatistical model. The intrinsic hypothesis of the random field is that the increments ρ(x+h)−ρ(x) are stationary: E[ρ(x+h)−ρ(x)] = 0, and the variogram is Var [ρ(x+h)−ρ(x)] = 2γ(h).

We compute the estimation variance by σ² = Var(ρ−ρ̄) = 2γ̄(I,V)−γ̄(V,V)−γ̄(I,I), where I describes the survey line, V the domain of interest, and γ̄(I,V) the mean of γ(y−x), when y∈I and x∈V. We use a variogram of the exponential class with a nugget effect:

A1

We estimate the parameters (and the nugget effect) by fitting the parametric model (A1) to an empirical variogram of observed densities, using weighted least squares with weights suggested by Cressie (1985). The empirical variogram was computed by the modulus estimator of Hawkin and Cressie [Cressie, 1993, Equation (2.4.12)]. We used the geoR package (Ribeiro and Diggle, 2001) for variogram fitting.

As the true densities are not observed, the procedure above is instead performed on the simulated densities ρ̂_i*, giving the variance estimate σ̂²*. The quantity ρ̄** in step (iv) of the proposed procedure (Methods) is simulated as ρ̄** = ρρ̄*+ε*, where ρ̄* is the average of the ρ̂_i*s, and ε* is drawn from a normal distribution with mean 0 and variance σ̂²*.

This procedure is performed per subarea. As it is more computer-intensive than the block bootstrap method, it is repeated only B = 500 times.