StoX applied to cod and haddock data from the Barents Sea NOR-RUS ecosystem cruise in autumn

North East Arctic (NEA) cod ( Gadus morhua ) and NEA haddock ( Melanogrammus aeglefinus ) are the largest stocks of these species in the world and have great economic value. In 2018, the Norwegian catches alone had a first-hand landing value of over 8 billion NOK. Indices from three surveys are used for the catch advice of both these stocks: the “ Barents Sea NOR-RUS ecosystem cruise in autumn” (BESS) , the “ Barents Sea NOR-RUS demersal fish cruise in winter” (winter survey), and the Russian bottom trawl survey for demersal fish (Lepesevich and Shevelev 1997, discontinued after 2017). Indices from these surveys are used as tuning series for cod ages 3-12 and haddock ages 3-9 in the stock assessment models (SAM, www.stockassessment.org ) and the indices for pre-recruits are used as input to the recruitment forecast models (ICES 2019).

StoX is a software developed at the Institute of Marine Research for marine survey analysis and index calculation. StoX is available for free ( ftp://ftp.imr.no/StoX/Download/ ) and is relatively well documented (Johnsen et al 2019). StoX is currently used for bottom trawl index calculation from winter survey data (Mehl et al 2016) but not from BESS data. BESS indices for cod and haddock are calculated with the BIOFOX software developed and run by the Polar Branch of FSBIU VINRO, Murmansk Russia (Prozorkevich and Gjøsæter 2014). Here we use StoX to calculate BESS swept area indices for haddock and cod 2004-2017. StoX needs a predefined strata system to estimate swept area indices. Whereas the winter survey has had a predefined strata system and a stratified survey design since the onset in 1981, BESS, starting in 2004, has had no predefined strata system. BESS has employed a regular grid design, with notable deviations, especially in deeper areas and around Svalbard.

Here we develop a strata system and use the StoX software to estimate haddock and cod abundance indices from “Barents Sea NOR-RUS ecosystem cruise in autumn”.

2.1 - Input data and strata system

BESS started as a compilation and extension of several surveys with different aims, gears and survey designs. BESS has been running annually August – October 2004-2018 with 2-3 Norwegian vessels and 1-2 Russian vessels (Appendix Table 1). The number of survey days and the number of bottom trawl stations has declined over the years (Appendix Table 1). With the exceptions shown and described in Appendix Figures 1-15, the whole Barents Sea shelf has been covered with stations set out in a regular grid with inter-station distance of 35 nautical miles. From 2017 and onwards the stations in the station grid was set out using the Albers equal-area projection. We do not include data from 2018 further due to a nearly complete lack of coverage in the eastern Barents Sea (Appendix Figure 15). The 2018 BESS indices were not used in the assessment (ICES 2019). BESS monitors plankton, oceanographic conditions, pelagic and demersal fish, pollution, benthos, shellfish, sea birds and marine mammals using an array of sampling methods. Here we used only data from bottom trawls. The standard towing time for bottom trawls is 15 minutes at 3 knots, but in deeper areas along the shelf break, tows has been longer due to lower fish densities. The bottom trawl gear, rigging, handling of the catches etc are similar as for the winter survey and are detailed in Mehl et al (2016, see also Appendix Table 2). The exception is that from 2018 and onwards gonadal weights and stages have not been registered at BESS. The number of aged and length measured cod and haddock is shown in Appendix Table 1.

Strata are used in StoX when filling in missing ages for length groups and in bootstrapping (see below). The strata should therefore ideally be homogenous with respect to the relationships between age and length and age and weight, as well as age specific abundances. We defined a basic strata system from the geographical distribution patterns of cod and haddock, observed at BESS and aggregated over the period 2009-2018 (Figure 1).

Figure 1. The basic BESS strata system used here with modifications by year accounting for the variable coverage and station densities in some areas and years (Appendix Figures).

The 250 m depth contour was used as guideline for the borders between strata, as a preliminary analysis show that this depth separates quite well between low (deeper) and high-density (shallower) areas for cod and haddock. The 500 m depth contour was used as a limit for the strata system towards the Norwegian Sea in the West and the Polar Ocean to the North as the densities of cod and haddock at greater depths are very low. The southernmost strata were divided by the border separating the Russian-Norwegian Exclusive zones, since the survey effort varied between Russian and Norwegian vessels in this area some years.

The result is a basic strata system that served as a starting point when generating the yearly strata systems used for estimation (Figure 1), modified to account for irregular sampling, ice hindrance, and limited survey coverage due to other reasons (Appendix Figures 1-15). The size of the strata was defined so that each stratum ideally contained at least 20 trawl hauls (although difficult to satisfy some years). We also modified the basic strata system to account for area specific differences in station density, since if not, areas with denser station coverage grid will be given more weight when calculating totals for the whole area.

2.2 - StoX calculations and settings

The data was downloaded from: https://datasetexplorer.hi.no/apps/datasetexplorer/v2/navigation and the folder “cruise series\Barents Sea NOR-RUS ecosystem cruise in autumn”.

The StoX projects presented here is stored in folders:

Surveytimeseries\Barents Sea Northeast Arctic haddock bottom trawl index in autumn

Surveytimeseries\Barents Sea Northeast Arctic cod bottom trawl index in autumn

We used the same station filter in StoX to select stations as Mehl et al (2016), except that stations with the gear code “3293” (“Tromsø rigging”) was included in the StoX filter here (Tromsø rigging is used for bottom trawls at soft bottoms and is described in Mehl et al 2016 and in Appendix Table 2).

Furthermore, we used 1 cm length groups rather than 5 cm length groups as in Mehl et al (2016) when applying the length dependent sweep width correction. The equations for length dependent sweep width for cod and haddock can be found in Mehl et al (2016).

2.2.1 - Swept area indices by length

The swept area density ( ρ , individuals per square nautical mile, inds nmi ^-2 ) by stratum (k), station (s) and length group l (1 cm), is given by

, (eqn 1)

where f _k,s,l is the number of individuals standardized over a towing distance of 1 nmi by k, s and l, and sw _l is the adjusted swept width in nmi’s by length group calculated using

, (eqn 2)

where EW _l is the length dependent effective swept width

The abundance ( N , inds) by l and k is calculated using

(eqn 3)

where A is stratum area (nmi ² ), and ρ _k,l is the average swept area density by l and k, given by

(eqn 4)

where n is number of stations.

2.2.2 - Swept area indices by age

The sampling protocol for the BESS survey is to sample one individual from each 5 cm length group at each trawl station for aging and individual weights. A two-stage conversion process is used to convert the abundance of fish by length group to abundance of fish by age group.

Firstly, the abundance ( N _k,l ) by length group l (5 cm) and stratum k is distributed the length-measured individuals (j) to generate so-called “Super-individuals” (super-individuals represent fractions of a total, our use corresponds to a probability based design where is the inverse of the inclusion probability for a single fish sample), each representing an abundance estimated as:

, (eqn 5)

where

(eqn 6)

and m is the number of length-measured individuals

Secondly, in instances where a super-individual is not aged, the missing age is filled in by a random data imputation. The imputation of missing age is principally carried out at the station level, randomly selecting the value from aged super-individuals within the same length group. If no aged super-individual is available at the station level, the imputation is attempted at strata level, or lastly on survey level. In instances where no age information is available at any level for a specific length group, the abundance estimate is presented with unknown age (Johnsen et al 2019).

2.2.3 - Length and weight at age

Length and weight at age was calculated using the weighting factors defined in eqn 6 (the “super -individuals”).

2.2.4 - Uncertainty and consistency of abundance indices

Uncertainty was estimated as the coefficient of variation (ratio of standard deviation to the mean, CV). StoX calculates CV using bootstrap runs by stratum, treating each trawl station as the primary sampling unit. Here we used 500 bootstrap runs

Consistency was quantified here as the correlation between the log abundance of all cohorts at age a in year y with the log abundances of the same cohort at older ages and later years. It is most common to consider only the correlation between the abundance of the cohorts at age a in year y, with the abundance at age a+1 in y+1, but here we present the correlations for the cohorts at all ages.

3.1 - Estimates at age

3.1.1 - Abundance at age

The numbers at ages in millions are given for cod and haddock in Table 1 and 2 respectively.

Table 1. Cod numbers at age in millions, swept area estimates, mean values from 500 bootstraps.

*2014 not used in cod assessment due to poor coverage of the Northern Barents Sea, see App. Fig 11.

Table 2 Haddock numbers at age in millions, swept area estimates, mean values from 500 bootstraps.

The numbers at age (ages 1-15 for cod and 1-12 for haddock) for the cohorts 2000-2016 are plotted in Figure 2 (cod) and Figure 3 (haddock).

Figure 2. Log abundance of cod (individuals), swpet area estimates from the Barents Sea Ecosystem Survey. Each line represent a cohort, first cohort the 2000-cohort, last cohort 2016. The numbers correspond to age (1-15).

Figure 3. Log abundance haddock (individuals), swept area estimates from the Barents Sea Ecosystem Survey. Each line represent a cohort, first cohort the 2000-cohort, last cohort 2016. The numbers on the lines are the ages from 1 to 12. The dip at age 9 in 2012 might be caused by errors in age reading e.g. misreading haddock age 9 as haddock age 8 or 10. Age 9 in 2012 is from a weak cohort, and the uncertainty is high (see Table 7). Also note the upwards jump for many age classes in 2016, a year with lack of coverage in central/eastern Barents Sea.

3.1.2 - Mean length at age

Length at age of cod and haddock are given in Table 3 and 4, respectively.

Table 3. Mean length at age of cod (cm), data from bottom trawls Barents Sea Ecosystem survey.

Table 4. Mean length at age of haddock (cm), data from bottom trawls Barents Sea Ecosystem survey.

3.1.3 - Weight at age

Weight at age of cod and haddock are given in Table 5 and 6, respectively.

Table 5. Mean weight at age of cod (g), data from bottom trawls Barents Sea Ecosystem survey.

Table 6. Mean weight at age of haddock (g), data from bottom trawls Barents Sea Ecosystem survey.

3.2 - Quality evaluation

3.2.1 - Coefficients of variation

Table 7. Cod. Estimates of coefficient of variation (CV) for swept area abundance - see Table 1.

*not included in cod tuning

Table 7. Haddock Estimates of coefficient of variation (CV) for swept area abundance - see Table 1.

3.2.2 - Consistency

The consistency was quantified as the correlation between the log abundance of each cohort at age a with the log abundances of the same cohort at older ages, and is presented here for cod ages 1 to 15 years (Figure 4 and 5), and haddock aged 1 to 12 (Figure 6).

Figure 4. Cod internal consistency Barents Sea ecosystem survey 2004-2017, ages 1-15. The ages are shown in the grey squares on the diagonal. The numbers in the lower part of the matrix are the correlation coefficient between the log10 abundance of a cohort of age a with the log10 abundances of the same cohort at older ages. The upper part of the matrix plots the log10 abundance of a cohort of age a against the log10 abundances of the same cohort at older ages, plotted with the fitted regression line. Colours goes from deep red (high consistency, >0.8), to pale yellow (low consistency <0.1) and white, negative correlation.

In Figure 5, 2014 was removed from the cod data, since this year had poor coverage to the North (see Appendix Figure 11) Indices from BESS in 2018 is not included in the tuning series for the cod assessment (ICES 2019).

Figure 5. Cod internal consistency Barents Sea ecosystem survey 2004-2013 and 2015-2017, ages 1-15. The year 2014 was excluded due to lack of coverage in the northern Barents Sea (see appendix Figure 11) and is not used in the cod assessment. The ages are shown in the grey squares on the diagonal. The numbers in the lower part of the matrix are the correlation coefficient between the log10 abundance of a cohort of age a with the log10 abundances of the same cohort at older ages. The upper part of the matrix plots the log10 abundance of a cohort of age a against the log10 abundances of the same cohort at older ages, plotted with the fitted regression line. Colours goes from deep red (high consistency, >0.8), to pale yellow (low consistency <0.1) and white, negative correlations.

Figure 6. Haddock internal consistency Barents Sea ecosystem survey, 2004-2017, ages 1-15. The ages are shown in the grey squares on the diagonal. The numbers in the lower part of the matrix are the correlation coefficient between the log10 abundance of a cohort of age a with the log10 abundances of the same cohort at older ages. The upper part of the matrix plots the log10 abundance of a cohort of age a against the log10 abundances of the same cohort at older ages, plotted with the fitted regression line. Colours goes from deep red (high consistency, >0.8), to pale yellow (low consistency <0.1) and white, negative correlation.

The Barents Sea ecosystem survey (BESS) has many applications, one of the most important is the use of abundance indices from this survey in the annual cod and haddock assessments. A proper, comprehensive description and overview of the changes in survey design and sampling methods for BESS, similar to what has been provided for the winter survey by e.g. Jakobsen et al (1997) and Mehl et al (2016), has been lacking. We have here presented and compiled some of the relevant information on changes in design and sampling (see Appendix Table 1 and 2, Appendix Figures 1-15). Most notable, BESS has had incomplete spatial survey coverage due to ice hindrance in 2014, and due to technical problems and restricted access and trawl permits particularly in 2015, 2016 and 2018. In 2018 the coverage was too poor to be consider for index calculation for the cod and haddock assessments. We constructed a simplified strata system based on average densities of cod and haddock, that we modified by year to account for variation in sampling intensities. We used this strata system in StoX and Rstox to estimate abundance indices for cod and haddock and lengths and weights at age. We also quantified the uncertainty and consistencies of the survey indices.

Uncertainty measured as the coefficient of variation (CV) were much higher for haddock than for cod. This might be due to overall lower haddock abundance, a more restricted spatial distribution, and that haddock is patchier distributed. Higher CV for haddock than cod were also evident for the winter survey (Mehl et al 2016). Compared to the winter survey, that has a denser station grid, the BESS CVs for both species were higher, but this was clearly much more pronounced for haddock. Some years appeared to have higher uncertainties than others. Cod uncertainties for older ages were high in 2010, a few very high catches in the eastern Barents Sea is likely responsible for this. Also, CVs for all ages of haddock were >0.2 in 2016, a year with a large hole in coverage in the central/eastern Barents Sea. Despite high CVs for haddock, the internal consistency was quite good also for haddock and comparable to the consistencies for cod. The survey tracked the declines in cohorts over time, and the correlation between the log abundance of a cohort at age a in year y with the log abundances of the same cohort at age a +1 in year y +1, was ≥ 0.73 for both cod and haddock of ages used in tuning.

Abundance estimates from BESS used in the Arctic Fisheries Working Group (ICES 2019) for assessment tuning is provided by PINRO, Russia and estimated using the BIOFOX software (Prozorkevich and Gjøsæter 2014). The internal consistency is higher for the haddock abundances estimated using BIOFOX compared to StoX, whereas consistency is almost identical for cod (Johannesen et al 2019). BIOFOX included extra demersal stations set out at acoustic registrations some years (Russian vessels only) and BIOFOX uses spatial interpolation by depth and distance (Prozorkevich and Gjøsæter 2014). This method appears to partly adjust for holes in the survey coverage and the large distance between the stations relatively to the patchy distribution of haddock.

We recommend that the BESS survey design should be kept constant from year to year, to make inter-annual variation in distributions etc easier to detect. If changes are to be made to the design, sampling gear etc, we recommended that the changes are based on relevant data analyses and biological knowledge and that they are well documented. BESS is a cooperation between Russia and Norway. Any change in the survey design and sampling must be agreed by both parties. Sampling methods are harmonized between the two countries, still there are issues related to e.g. sampling gear (Appendix Table 2), frequency (and use) of stations at dense acoustic registrations and data exchange and conversion that should be improved. A denser station grid in the main distribution area of haddock would probably improve the indices for this species and should be considered. More sophisticated methods for index calculations than presented here, e.g. geostatistical methods/delta-gam methods taking e.g. depth or other habitat characteristics into account, could perhaps amend some of the issues related to area and year variation in sampling intensities, large interstation distance and holes in the coverage, and should be considered. However, incomplete survey coverage (reduced extent and holes) will always increase uncertainty, there are no straightforward way to account for it, and it should therefore be avoided.

ICES (2019). Arctic Fisheries Working Group (AFWG). ICES Scientific Reports. 1:30. 930 pp.

Jakobsen T, Korsbrekke K, Mehl S and Nakken O (1997). Norwegian combined acoustic and bottom trawl surveys for demersal fish in the Barents Sea during winter. ICES CM Y: 17.

Johannesen E, Johansen GO and Prozorkevich D (2019). WD_03. Cod and haddock abundance indices by age from the ecosystem survey: comparing current indices from BIOFOX and new indices from StoX. AFWG 2019.

Johnsen E ,  Totland A ,  Skålevik Å , et al.  (2019). StoX: An open source software for marine survey analyses . Methods in Ecology and Evolution 2019 ;  10 :  1523 –  1528 .  https://doi.org/10.1111/2041-210X.13250

Lepesevich YM and Shevelev M.S. (1997). Evolution of the Russian survey for demersal fish: from ideal to reality. ICES CM Y: 09.

Mehl S, Aglen A and Johnsen E (2016). Re-estimation of swept area indices with CVs for main demersal fish species in the Barents Sea winter survey 1994 – 2016 applying the Sea2Data StoX software. Fisken og Havet 10 2016.

Pavlenko A, Likhoshapko A, Prozorkevich D, Prokhorova T, Engås A, Aasen A and Eriksen E (2014). Report from the meeting between experts from PINRO and IMR to discuss sampling trawl methodology. Murmansk 28-29. January 2014. Unpublished report. 4 pp.

Prozorkevich D and Gjøsæter H (2014). WD_02 cod_BESS_assessment. AFWG 2014.

Appendix table 1. Year: survey year; NO : name of participating Norwegian vessels; RU: Russian vessels; Days: the total number of survey days in that year – nb taken partly from survey reports and partly from data and therefore not 100% accurate; Stations: the number of all valid bottom trawl hauls, includes also hauls taken outside the strata system, see Appendix figures 1-15, Measured: total number of aged/length measured/ weighed individuals

*Jan Mayen renamed Helmer Hanssen in 2012

**Commercial vessel Atlantic Star took some pelagic hauls as part of the survey

***Commercial vessel Christina E and Eros replaced GO Sars in 2011 and 2016

****Access restricted to part of the REZ for Vilnius, so survey was extended one month for this vessel

Table 2 a. Overview of rigging and trawling procedure for the Campelen 1800 trawl used by IMR and PINRO, from Pavlenko et al (2014).

Table 2 b. Overview of rigging and trawling procedure for the Campelen 1800 trawl used by IMR and PINRO, from the March meeting protocol 2019.

Appendix Figure 1. 2004. Stations grid set out according to 0-group survey design with Mercator projection. Approximately 35 nautical mile grid station distance. Denser grid was due to Greenland halibut juvenile investigation in Storfjordrenna, and west and east of Svalbard, and due to this the “Great Bank” stratum is split to account for the denser sampling in the western part of this stratum. The extent of the “ South East” stratum is reduced due to limited sampling coverage .

Appendix Figure 2. 2005. Stations set out in a modified UTM with approximately 35 nautical mile grid station distance. Denser grid was due to Greenland halibut juvenile investigation in Storfjordrenna, and west and east of Svalbard, and shrimp investigations in the Bear Island trough. Due to this, the “Great Bank” stratum is split to account for the denser sampling in the western part of this stratum, and the eastern part of this stratum is added to the “North East” stratum to obtain sufficient sample size. The “South West” stratum was split in two for the same reasons. Northern and eastern edges was reduced due to limited sampling coverage.

Appendix Figure 3. 2006. Stations set out in a modified UTM with approximately 35 nautical mile grid station distance. Denser grid was due to Greenland halibut juvenile investigation in Storfjordrenna, and west and east of Svalbard and shrimp investigations in the Bear Island trough. The “Great Bank”, “Bear Island Through”, “South East” (flatfish investigations) and “South West strata were modified to account for denser sampling. Northern and eastern edges was reduced due to limited sampling coverage.

Appendix Figure 4. 2007. Stations set out in a modified UTM with approximately 35 nautical mile grid station distance. Denser grid was due to Greenland halibut juvenile investigation in Storfjordrenna, and west and east of Svalbard, and shrimp investigations in the Bear Island trough. The “Great Bank” stratum was split to account for denser sampling. North-eastern and eastern edges was reduced due to limited sampling coverage.

Appendix Figure 5. 2008. Stations set out in a modified UTM with approximately 35 nautical mile grid station distance. Denser grid was due to Greenland halibut juvenile investigation in Storfjordrenna, and west and east of Svalbard and shrimp investigations in the Bear Island trough. The sampling effort was regarded as fairly regular compared to earlier years, so the only modification was reduction of the northern edges due to limited survey coverage. This year suffered a strong reduction in Norwegian survey time due to budget cuts.

Appendix Figure 6. 2009. Stations set out in a modified UTM with approximately 35 nautical mile grid station distance. Denser grid west of Svalbard due to steep depth gradients. Reduction of the northern edges due to limited survey coverage. Russian investigations in Kara Sea (North east) not included in StoX calculations.

Appendix Figure 7. 2010. Stations set out in a modified UTM with approximately 35 nautical mile grid station distance. Reduction of the northern edges due to limited survey coverage. Russian investigations in Kara Sea (North east) not included in StoX calculations.

Appendix Figure 8. 2011. Stations set out in a modified UTM with approximately 35 nautical mile grid station distance. Denser grid west of Svalbard due to steep depth gradients. Reduction of the northern edges due to limited survey coverage. Russian investigations in Kara Sea (North east) not included in StoX calculations.

Appendix Figure 9. 2012. Stations set out in a modified UTM with approximately 35 nautical mile grid station distance, except around Svalbard. Denser grid west of Svalbard due to steep depth gradients. Reduction of the north-eastern edges due to limited survey coverage. Russian investigations in Kara Sea (North east) not included in StoX calculations.

Appendix Figure 10. 2013. Stations set out in a modified UTM with approximately 35 nautical mile grid station distance, except around Svalbard. Denser grid west of Svalbard due to steep depth gradients. Reduction of the north-eastern edges due to limited survey coverage.

Appendix Figure 11. 2014. Stations set out in a modified UTM with approximately 35 nautical mile grid station distance, except around Svalbard. Denser grid west of Svalbard due to steep depth gradients. Reduction of the northern edges due to limited survey coverage. NB! Ice restrictions in the north – this year the survey indices were not used in the cod assessment.

Appendix Figure 12. 2015 Stations set out in a modified UTM with approximately 35 nautical mile grid station distance, except around Svalbard. Denser grid west of Svalbard due to steep depth gradients. No coverage in the loop hole due to unclear jurisdictions. Reduction of the northern edges due to limited survey coverage.

Appendix Figure 13. 2016. Stations set out in a modified UTM with approximately 35 nautical mile grid station distance, except around Svalbard. Denser grid west of Svalbard due to steep depth gradients. No coverage in loop hole due to unclear jurisdictions. Russian military exercise limited survey access to the south eastern BS. Reduction of the north-eastern edges due to limited survey coverage.

Appendix Figure 14. 2017. Stations set out with approximately 35 nautical mile grid station distance, according to Albers equal-area projection (also called: Albers Equal-Area Conic Projection) with the following parameter settings: Centre latitude: 75 N, Centre longitude: 30 E, 1st standard Latitude: 70 N, 2nd standard latitude: 80 N. Modified station grid in the western part to avoid deeper stations than 500 m, and around Svalbard due to restrictions in natural reserves. Some additional stations in “Spitsbergen Bank”, “Svalbard” and “Great Bank” strata to cover Greenland halibut. Slight reduction of the north-eastern and south-eastern edges due to limited survey coverage.

Appendix Figure 15. 2018. Stations set out with approximately 35 nautical mile grid station distance, according to Albers equal-area projection (also called: Albers Equal-Area Conic Projection) with the following parameter settings: Centre latitude: 75 N, Centre longitude: 30 E, 1st standard Latitude: 70 N, 2nd standard latitude: 80 N. Modified grid in the western part to avoid deeper stations than 500 m, and around Svalbard due to restrictions in natural reserves. Some additional stations in “Spitsbergen Bank”, “Svalbard” and “Great Bank” strata to cover Greenland halibut. Slight reduction of the north-eastern edges due to limited survey coverage. Lack of coverage in east was due technical problems with the Russian vessel. No attempt was done to calculate survey indices in 2018, due to this severe coverage problem.