Methodology:

Hurdle model approach:

To model the distribution of Antarctic krill (hereafter krill) we applied a two-part hurdle model or zero-altered model (Zuur et al. 2009), which allows for heavily skewed data distributions. The first model predicts the probability of presence of krill. To do this, the krill density data were transformed into a binary zero/non-zero form (n = 2709) and modelled against the environmental covariates using a binomial generalised additive model (GAM; Wood (2017)) with a logit link function. The second model investigates the relationship between non-zero data (i.e. presence-only data, n = 1833) and environmental covariates. This was carried out using a GAM with a Gaussian distribution and the default identity link function. Based on exploratory density plots, the presence-only data were log transformed to follow a normal distribution. Finally, outputs from both models were multiplied together. This allowed us to identify where krill were both likely to be present and occur at high densities.

Model fitting and selection:

All GAMs were fitted using the R package 'mgcv' (Wood 2019), with a Restricted Maximum Likelihood (REML) optimisation method to estimate splines, and penalised thin plate regression splines on all smooth terms (Eilers and Marx 1996). To reduce model overfitting the basis dimension (k) was limited to between 3 and 6 with the optimum number guided by edf values and associated p values reported in mgcv's gam.check, and by visualising the partial effects plots for each covariate with the raw data.

For both the binomial and Gaussian GAMs, model selection followed a forward stepwise selection approach with five-fold cross validation. Specifically, each environmental covariate was modelled against the response variable independently and repeated five times, each time withholding a different random subset of data (fold) for evaluation. The model coefficients from each run were used to predict the outcome for the withheld fold and performance metrics of the prediction - Root Mean Squared Error (RMSE) and R2 - were extracted. The best performing covariate (i.e. lowest RMSE and highest R2 averaged over five folds) was retained within the model. This selection process was repeated allowing for all possible combinations of environmental covariates at their different temporal scales (sample and decadal). At each iteration, the retained set of covariates were assessed for collinearity using Pearson correlation coefficients and Variance Inflation Factors (VIF), and for concurvity using the worst-case measure of overall concurvity for each smooth. If issues were identified (Pearson's r > 0.7, VIF >3, concurvity >0.8) the next highest-ranking covariate was selected. Forward selection continued until model performance metrics plateaued and/or issues of collinearity and concurvity could not be overcome. Predictions from the final Gaussian GAM were back transformed to obtain outputs on the original density scale. Once the final set of covariates was selected, predictions for the probability of occurrence and estimated krill density were projected onto year-specific grids at the scale of Subarea 48.2. The mean and +/-1 standard deviation of predictions and their product (interpreted as the krill density weighted by the probability of occurrence) across all years were generated.

Data collection:

All analyses were carried out in R version 4.3.3.

Data quality:

The quality of model outputs are, in part, dependent on the data used for model input, i.e. the krill density data and environmental covariates.

Krill density data:

As part of the ongoing Norwegian scientific contribution to monitoring distribution, abundance and population characteristics of Antarctic krill, an acoustic survey takes place annually in waters surrounding the South Orkney Islands (Longitudinal stratum boundaries at 43.5 degrees W and 48 degrees W, and latitudinal boundaries at 59.67 degrees S and 62 degrees S). In the present study we use data from the first ten-year time period of this survey (2011 to 2020). All survey transects occurred between January and February each year, but in two of the years - 2013 and 2015 - sea ice prevented full completion of the survey (Krafft et al. 2018). The surveys were conducted by using the Norwegian commercial fishing vessels 'Saga Sea' and 'Juvel' (Aker Biomarine ASA, Oslo, Norway and Rimfrost AS, Fosnavag, Norway) as research platforms except in 2019 when RV 'Kronprins Haakon' was used.

The acoustic data used in this study were collected using calibrated hull mounted Simrad echo sounders. We used the swarm-based approach for acoustic target identification of krill which is recommended by CCAMLR when the 38, 120 and 200 kHz frequency combination is not available. Details on the processing procedures used to determine krill density are reported and evaluated by Skaret et al. (2023). The retained nautical area scattering coefficient (NASC) allocated to krill per nautical mile was converted to biomass density (g m-2, hereafter referred to as density) using full Stochastic Distorted Wave Born Approximation (SDWBA) model runs to estimate backscattering cross-sectional areas (delta) for each krill length group of 1 mm increment present in the sample.

Each acoustic sample value was matched to the covariate raster data according to the latitude, longitude and year of collection. Finally, the combined krill density-covariate dataset was aggregated to the same spatial resolution as the environmental covariates (0.04×0.04 decimal degrees) by calculating the mean values within each grid cell for each year. This was done to avoid pseudo-replication given multiple acoustic samples within the same grid cell of environmental data, and to reduce any effect of spatial autocorrelation in model residuals.

Environmental covariate data:

Twelve environmental covariates were identified as candidate explanatory variables for the model. These included three static variables, i.e. unchanging with survey year: water depth (bathymetry), bathymetric slope and distance from shelf break defined as the 500m isobath (where values on-shelf were positive, and those off-shelf were negative). The nine remaining variables were dynamic across survey years: distance from sea ice edge (defined as 15 percentage ice concentration), seven sea surface variables (temperature, mixed layer thickness, sea surface height above geoid, salinity, chlorophyll a, primary productivity, and geostrophic current velocity) and bottom temperature. Raster grids of all covariates were obtained from a combination of empirical observations and model re-analyses (see Section 9. Related datasets). For each of the dynamic covariates, two different temporal scales were extracted. These were: i) a sample scale which averaged conditions during the sampling months (January to February) independently for each year; and ii) a decadal scale climatology which was the average of summer conditions (January to March) between 2011 and 2020.