Environmental conditions on the Pacific halibut fishing grounds obtained from a decade of coastwide oceanographic monitoring, and the potential application of these data in stock analyses

Oceanographic instrument deployment

Seabird SBE19plus profiling units were deployed just prior to hauling the longline fishing gear at each FISS station from 2009 to 2018. The FISS was conducted using chartered longline fishing vessels that deployed standardised longline gear and bait during the summer months (primarily June, July and August) along a grid of fishing stations located on the continental shelf. Profiles were collected during periodic range expansions into areas and depths not ordinarily sampled, but only data from those stations scheduled to be profiled in all years (i.e. core stations) were retained for this analysis. The profiling instruments were equipped with sensors that collected measurements of pressure (decibars were used as a shallow-water proxy for depth in metres), conductivity (from which salinity in practical salinity units, PSU, were calculated), temperature (°C), dissolved oxygen (DO, mL L⁻¹), pH and fluorescence (mg m⁻³; indicating chlorophyll-a concentration). A pump ensured consistent sampling throughout the profile. An anchor weight was attached by 15 m of the fishing vessel’s buoy line to the bottom of the stainless-steel cage that surrounded the instruments. Attached to the top of the cage was a float assembly and buoy line that remained secured to the fishing vessel throughout deployment (Fig. 2).

Fig. 2.

Schematic of the water-column profiler set-up used for oceanographic sampling during the IPHC fishery-independent setline survey. The unit is equipped with floats attached to the top of a protective frame and a weight attached by buoy line to the bottom of the frame to ensure a vertical descent and to minimise the risk that it will affect the bottom.

Once on station and prior to hauling the longline gear, the profiler assembly was deployed. The float and weight rigging was designed to accomplish a vertical descent of 1–2 m s⁻¹. Data were recorded at a rate of four readings per second, and were uploaded periodically to onboard computers. Both in-season and annual maintenance was performed on the instruments as recommended by the manufacturer. In 2012, the anchors of the profiling units were equipped with an additional depth sensor (Star-Oddi DST logic model) and analysed in conjunction with the profiler-depth and vessel-depth sounder instruments to estimate proximity of the profiler to the sea floor at its deepest reading. The analysis produced an estimate that the instruments came within 5–15 m from the bottom, on average, depending on total depth of the station and likely other factors such as horizontal-current strength (Sadorus et al. 2016).

Data conversion and preliminary data-quality assessment and control (QC) of the profile data were accomplished using Seabird Scientific’s software SBE Data Processing (ver. 7.26.7, see https://software.seabird.com/). Profiles were averaged to 1-m bins. No calibration water samples were collected because of FISS limitations, and variables were not laboratory calibrated. Data were visually inspected and QC-edited as needed post-season. Initial accuracy for the instruments used was as follows: pressure (+0.1% of full-scale range), temperature (+0.005°C), conductivity (+0.0005 S m⁻¹), dissolved oxygen (+2% of saturation) and chlorophyll-a (sensitivity of 0.02 μg L⁻¹). Temperature data rarely fell outside set parameters and the sensors proved stable. Most problems with salinity data occurred in the upper water column, and values were stable in the deeper water layer where the majority of Pacific halibut habitat occurs. Oxygen and chlorophyll-a values would likely have been improved with in-season calibration, but were stable overall, i.e. the sensitivity of the instruments did not appear to drift considerably between annual calibrations. The pH data are useful as a general indicator, but are the least reliable of these data and were not included in this analysis, because the pH sensors were known to have drift problems. Full profiler datasets for each year are available for download from the IPHC website (see https://www.iphc.int/datatest/data/water-column-profiler-data), and near-bottom measurements for temperature and DO are also available on the website using the FISS data-visualisation tool (see https://www.iphc.int/data/FISS-catch-per-unit-effort).

Oceanographic and catch data

Because Pacific halibut is demersal, the deepest readings at each station were the most indicative of conditions that the animals in that area were experiencing. Profiles were excluded from this analysis if it was suspected that the instrument failed to reach the bottom. The profiler data for pressure, temperature, DO and salinity were averaged and binned into 1-m depth increments, checked for accuracy, and the deepest binned values for each profiler cast were used in this analysis. Averaging of near-bottom data included only those stations that had been profiled successfully five or more times during the 10-year period. Profiler data from 2009 to 2014 were pre-processed by IPHC and, subsequently, underwent rigorous quality control (QC) processes by the NOAA Pacific Marine Environmental Laboratory (PMEL) in collaboration with the University of Washington Joint institute for the Study of the Atmosphere and Ocean (UW, JISAO) (Sadorus et al. 2016). More rigorous QC included the use of various tools to account for pressure reversals owing to weather and sea conditions, to align data relative to pressure, and to derive variables (salinity, density, σ_t, oxygen saturation). Data from 2016 to 2018 were processed at the IPHC only.

Chlorophyll-a measurements were averaged within each 1-m depth bin (there were typically 3–5 measurements per bin), so as to smooth any anomalous readings that may have been due to movement of the profiler in the current. Depth-integrated chlorophyll-a density at each station was estimated by using numerical integration with the trapezoidal rule over the station depth range, by using the 1-m averages as input data. Thus,

where x_i, i = 0, …, n, is the ith depth with mean chlorophyll-a value of f(x_i). In our case, Δx=b1-b0n=1, since n, the number of depth bins, is equal to the depth (m), which equals b₁ – b₀, the difference between the greatest and smallest (b₀ = 0) depth at a station. This calculation yielded chlorophyll (mg m⁻²) throughout the area. Station outliers were examined by inspecting the chlorophyll-a pattern across all depths for that station for reasonable values and by comparing to nearby stations. Data that did not fit into the pattern examined were removed. Prior to the calculation, slightly negative readings within the casts resulting from calibration differences of the sensor were changed to zero. Only those stations that had been fished five or more times during the 10-year period were included.

The spatio-temporal modelling included near-bottom temperature and DO data collected by the profilers and Pacific halibut catch data, which were recorded during the FISS and uploaded to the IPHC relational database post-season (catch data available at https://www.iphc.int/data/FISS-catch-per-unit-effort). Catch data were standardised to weight-per-unit-effort (WPUE) of Pacific halibut with fork length at least 81.3 cm (~32″, the commercial fishery size limit), i.e. pounds caught per 1800-foot (~1645-m) skate of fishing gear (International Pacific Halibut Commission 2020).

In total, 10,573 surface-to-bottom useable-quality profiles were collected for this analysis over the 10-year period from 2009 to 2018. In any given year, a small number of profiler deployments were missed because of poor weather or data were discarded post-season because of sensor malfunction. Following the editing process, 1268 unique stations remained for use in the analysis, each station having had five or more useable profiles over the study period. For the spatial-modelling example of hypoxia (defined here as DO < 1.4 mL L⁻¹) off the western coast of USA, all viable stations were included regardless of the number of times they were sampled over the 10-year period, such that all useable data available for 2017 and 2018 (years when additional stations were fished in IPHC Regulatory Area 2A) were utilised.

Statistical analysis

The oceanographic variables of temperature, salinity, DO and depth-integrated chlorophyll-a were compared across regions and over the 10-year time period by fitting linear models (function lm in the ‘stats’ package in R, ver. 3.6.1; see https://stat.ethz.ch/R-manual/R-devel/library/stats/html/00Index.html). Residual plots showed that the underlying model assumptions were met for the temperature and DO data and no prior data transformation was required. Chlorophyll-a data were highly skewed and normality and variance stability (homoscedasticity) were improved with a fourth-root transformation; although these transformations may lack a biological motivation, they were necessary for insuring the validity of the statistical inference. Model assumptions were violated for salinity data regardless of any transformation. Therefore, the untransformed-data curves were determined to provide the most useful description of the time-series. Of those variables examined, overall time trends and comparisons among biological regions were estimated by fitting a sequence of increasing complex linear models (Kunter et al. 2004), with F-tests and Akaike’s information criterion (AIC, Akaike 1973) used to select the best-fitting model. Analyses were performed in R (ver. 3.6.1, R Foundation for Statistical Computing, Vienna, Austria, see https://www.r-project.org/). We fitted models that assumed a simple linear relationship with time within each biological region over the 10-year data period, because our goal was to understand overall temporal trends (i.e. determining whether the means of these variables were increasing or decreasing, on average, over time), while recognising that variation over time is more complex than a simple linear process (see Supplementary Fig. S1–S4 for raw data plots).

Spatio-temporal modelling

Spatio-temporal modelling allows for the integration of environmental information as covariates within models for catch rate, while accounting for spatial and temporal correlation among the observations. Specifically, we used a delta-type model (Shelton et al. 2014) in which the probability of catching zero fish and the distribution of non-zero catches are modelled as connected spatio-temporal processes for recent data from IPHC Regulatory area 2A. This area was chosen to illustrate the potential usefulness of including environmental covariates in such models for Pacific halibut for the following two main reasons: its small size means we can fit spatio-temporal models quickly relative to other management areas, and notable recent hypoxic events mean we have a least one covariate that spans a range of values likely to have a clear effect on Pacific halibut density.

The dependent variable in our models is weight-per-unit-effort (WPUE), the catch weight of Pacific halibut with fork length 81.3 cm and over (the commercial size limit) divided by the number of 100-hook skates set. Whereas we recognise that WPUE is a function of both density and catchability, for the purpose of this analysis, we are assuming that WPUE provides an index of Pacific halibut density (see Discussion). Let w(s,t) be the IPHC FISS WPUE value at location s (a vector of coordinates) in year t, where s represents the spatial locations of the fished survey stations, taking values s₁, …, s_n (vectors of coordinates) and t = t₁, …, t_T. In our model, each s_i ε S², the set of points on the surface of a sphere. Data from the FISS contain observations of zero WPUE, owing to stations in low-density areas catching no Pacific halibut. Two new variables are defined, x(s,t) for presence or absence of Pacific halibut in the catch and y(s,t) for the WPUE value when Pacific halibut is present, as follows:

The NA indicates that y(s,t) is a random variable that can take only non-zero values, and is therefore undefined when w(s,t) = 0. The variable x(s,t) has a Bernoulli distribution, x(s,t) ~ Bern(p(s,t)), whereas a gamma distribution is used for the y(s,t), y(s,t) ~ gamma(a(s,t), b(s,t)), which has mean μ(s,t) = a(s,t) ÷ b(s,t). The gamma is parameterised in terms of the mean and variance, with only the former allowed to vary; the variance (σg2 = a(s,t) ÷ b²(s,t)) is assumed constant over space and time.

Next, let the ε(s,t) be a Gaussian field (GF) that is shared by both component random variables in the following way:

The parameter β_ε is a scaling parameter on the shared random effect. Environmental covariates are introduced into each model component through f_x() and f_y(), functions of a spatially and temporally indexed covariate data matrix (z) and covariate vectors β_x and β_y.

Temporal dependence is introduced through a simple autoregressive model of order 1 (AR(1)), as described in Cameletti et al. (2013), as follows,

where ρ denotes the temporal correlation parameter and |ρ| < 1. For a given year, t, the spatial random field, η(s,t), is assumed to be a Gaussian field with a mean zero and covariance matrix Σ. We assume a stationary Matérn model (Cressie 1993) for the spatial covariance model, which specifies how the dependence between observations at two locations decreases with an increasing distance.

Two covariates that appear linked to Pacific halibut density were included in the models for data from IPHC Regulatory area 2A, namely, bottom temperature and DO. Models were fitted with either linear or quadratic relationships for temperature, because data exploration throughout the species’ range indicated that Pacific halibut density can be negatively affected by both extremely low and very high temperatures. Thus, for temperature, the most complex model considered included βx1T+βx2T2 as part of f_x() in Eqn 1 and βy1T+βy2T2 as part of f_y() in Eqn 2, where T = T(s,t), the bottom temperature measurement at location s and year t. We also considered linear and quadratic models for the relationships with DO; even though we do not have much data for either the very lowest or highest values of DO, a linear function did not appear to adequately describe the relationship. Thus, β_x3DO + β_x4DO² were added as part of f_x() and β_y3DO + β_y4DO² as part of f_y(). Finally, we also fitted a model with a different intercept term for DO with values less than 0.9 mL L⁻¹, which past and current data imply is a cut-off level in DO below which Pacific halibut catches are almost always zero (Sadorus et al. 2014). This meant adding the term β_x5I(DO < 0.9) to f_x() and β_y5I(DO < 0.9) to f_y(), where I() is the indicator function taking value one when its argument is true and the value zero otherwise. A similar ‘breakpoint’ model was fitted by Essington et al. (2022) to sablefish (Anoplopoma fimbriadata) data off the western coast of USA, although they estimated the cut-off level within the model.

Models were fitted in a sequence (Table 1), starting with a ‘base’ model that included depth, a binary variable for latitude (1 if <40°N, 0 otherwise), and a temporal-trend effect; this is the model the IPHC fits annually for the purposes of producing a density index for use in the stock assessment. We model both the depth relationship and the trend in a flexible manner by a random-walk process as described in Webster et al. (2020). We note that depth, temperature and DO are correlated (Fig. 3), potentially making it more difficult to interpret the results of the modelling; both higher bottom temperatures and higher DO concentrations are found in shallower waters, although the lowest DO values also most commonly occur in depths from 0 to 200 m. This was managed through the sequential model fitting and selection process, which allowed us to compare the fit of models with and without each covariate. The latitude covariate is to account for extremely low-density habitat at the southern limits of the species’ range that was surveyed only once in 2017; its inclusion means that the low WPUE predictions south of 40°N persist throughout the time series, whereas predictions in this region would otherwise drift towards the long-term mean.

Fig. 3.

Plots of (a) depth against dissolved oxygen, (b) depth against bottom temperature and (c) bottom temperature against dissolved oxygen, as recorded by IPHC water-column profilers in the IPHC Regulatory Area 2A FISS from 2009 to 2018.

Some station data from the southern end of adjacent IPHC Regulatory Area 2B were included in the space–time models because spatial dependence extends beyond the regulatory-area boundaries and, therefore, such data provide additional information for IPHC Regulatory Area 2A models. In total, 1265 observations from 2009 to 2018 were used in the space–time modelling. Fig. 4 shows station locations, with circle size indicating the sampling frequency over the 10-year period.

Fig. 4.

Maps showing location of stations, with data used in the space–time modelling that were collected during the IPHC FISS. Circle size indicates the number of times each station was both successfully fished and produced valid water-column profiler data. Stations selected for model prediction are shown in orange.

Models were fitted in R using the R-INLA package (ver. 22.12.16, see https://www.r-inla.org/; Lindgren and Rue 2015), which uses a computationally efficient Bayesian approach to fitting spatial and spatio-temporal models. Further details are available in Webster et al. (2020), which includes a description of the process we used for creating a non-convex triangulated mesh for the habitat region of interest used by the INLA approximation in the models. Independent, non-informative prior distributions were used for the model parameters, including Gaussian priors for fixed effects and log-gamma priors for variance parameters.

The fits of the models were compared using the deviance information criterion (DIC, Spiegelhalter et al. 2002) and a k-fold cross-validation. For the latter, the data were divided into k = 5 subsets of equal size. The models were then refitted with one subset omitted, and the refitted model was used to obtain predictions of WPUE (posterior means) for each omitted observation. This process was repeated for each of the five subsets, providing a prediction for each observation in the original data set. An overall measure of model fit was then computed as the mean square error, as follows:

where WPUEobsi and WPUEpredi are the ith observed and predicted values of WPUE and N is the total number of observations used in the modelling.

To help compare the observed data with modelled covariate relationships, we selected 10 prediction stations spanning the latitude range of the core of the IPHC Regulatory Area 2A stock north of 42°N (shown as orange circles in Fig. 4). Means of posterior predictions for these stations for 2013 (chosen as it was an ‘average’ year for mean WPUE over the 2009–2018 period) over a range of covariate values were calculated and plotted for the best-fitting models.

Measured environmental variables

Mean station bottom temperatures over the study period (2009–2018) ranged from 0.4 to 10.1°C, with the highest temperatures consistently occurring off the western coast of USA, and off Vancouver Island, British Columbia, Canada in Biological Region 2, and around Kodiak Island, Alaska, USA, in Biological Region 3 (Fig. 5a). Coldest mean near-bottom temperatures by station were found in the Bering Sea, particularly at stations surrounding St Matthew Island and along the northern Bering Sea continental shelf edge in Biological Region 4. The highest regional mean temperature overall was in Biological Region 2 and the lowest in Region 4 (Fig. 5a; Tables 2 and 3). Region 4 displayed the greatest variability among biological region station means, and Region 4B displayed the smallest range of variability (Fig. S1). Mean temperatures over the study decade varied, but there was strong evidence that the overall trends varied among regions (Table 4a), with a greater average increase in Biological Region 4 than other regions (Fig. 6a). A recent study (Amaya et al. 2023) examining bottom temperatures in the Gulf of Alaska and Bering Sea supports the finding that the Bering Sea experienced higher average intensities of marine heatwave temperatures than in the Gulf of Alaska.

Fig. 5.

Maps showing mean oceanographic conditions for each station profiled for the Years 2009–2018 during the IPHC FISS: (a) near-bottom temperature (°C); (b) near-bottom dissolved oxygen (mL L⁻¹); (c) near-bottom salinity (PSU); (d) depth-integrated chlorophyll-a concentration (mg m⁻²). Maps generated using Ocean Data View software (ver. 4.7, R. Schlitzer, see https://odv.awi.de).

Fig. 6.

Annual means and trends of (a) near-bottom temperature (°C), (b) near-bottom dissolved oxygen (mL L⁻¹), and (c) depth-integrated chlorophyll-a concentration (mg m⁻²), from 2009 to 2018, as recorded by IPHC water-column profilers during the FISS. Note that (c) shows the back-transformed values.

Elevated temperatures corresponding to the marine heat wave of 2014–2016 (Bond et al. 2015; Gentemann et al. 2017) were detected in the data collected in the near-bottom waters of Biological Regions 3 and 4 (Fig. S1), although they were not apparent until 2015. Danielson et al. (2022) found a similar pattern of delayed warming in continental shelf bottom waters during this marine heat wave. Maximum temperatures in those regions decreased notably in the following years, whereas minimum temperatures had small increases (Region 4) or smaller decreases (Region 3), resulting in higher mean temperatures.

Mean station DO ranged from 0.7 to 7.6 mL L⁻¹ coast wide (Fig. 5b). The best-fitting linear model for DO included a year by region interaction (Table 4b), although the AIC difference between this model and the one without an interaction was relatively small. Thus, there was some evidence for regional differences in DO trends over the 10-year study period. Fig. 6b, showing model predictions for DO, shows slightly declining mean DO for Regions 3 and 4, whereas DO increases somewhat, on average, (with high uncertainty) for Region 4B. Region 2 had notably lower DO than did the other three regions (Table 2, Fig. 6b, S3), and additionally, hypoxia was found at both shallow and deep stations in Region 2, but only at deep stations in the other areas (Table 3).

Fig. 7 shows maps of DO off the western coast of USA, with the original data being smoothed using a simple exponential geostatistical model (Cressie 1993) so as to make the DO data easier to interpret visually. (Our intention is to visually highlight the hypoxic event in 2017 and therefore the exact choice of smoother is not important.) In 2017, a large cluster of stations had near-bottom DO at or below the Pacific halibut minimum threshold of 0.9 mL L⁻¹ (Sadorus et al. 2014) detected during the FISS. Pacific halibut was not captured at stations within the threshold outline, whereas in other years there were moderate densities within this region. In 2018, low concentrations of DO were also observed over a wide area further south of the 2017 occurrence. However, the affected zone typically had low Pacific halibut catch rates even in years when hypoxia was not present, and the hypoxic stations were interspersed with stations above the minimum DO threshold. The indication was that hypoxia did not have a notable distributional impact in 2018, as it did in 2017.

Fig. 7.

Estimated near-bottom dissolved oxygen concentration (DO, mL L⁻¹) off the western coast of USA in 2017 (left panel) and 2018 (right panel) with Pacific halibut weight-per-unit-catch effort values from the IPHC FISS overlaid with black symbols. Stations with DO of ≤0.9 mL L⁻¹ are outlined in white. Weight categories per skate are 25 lb (~11.3 kg), 50 lb (~22.7 kg), 100 lb (~45.4 kg) and 250 lb (~113.4 kg).

Coast-wide mean station near-bottom salinity ranged from 30.8 to 34.1. Biological Region 2 had the highest mean salinity and Region 3 had the lowest (Table 2, Fig. 5c, S3). In Region 2, stations with higher mean salinity were both offshore and inshore, with the southern stations somewhat uniform and high relative to stations further to the north that displayed increasing spatial variability. Region 3 stations that had higher mean salinity tended to be offshore and lower-salinity stations tended to be nearshore. In Region 4, lower-salinity stations were located near-shore in the Gulf of Alaska and around the Pribilof Islands and St Matthew Island in the Bering Sea. Higher-salinity stations were located further offshore at deeper depths and along the Bering Sea shelf edge. Region 4B stations tended to have salinity in the higher end of the range seen coast-wide, with minimal variability. Owing to its bimodal distribution, salinity data were not conducive to linear modelling, but overall salinity values were consistent over the 10-year time span, with no noteworthy temporal trending.

Mean station depth-integrated chlorophyll-a concentration ranged from 20.3 to 409.8 mg m⁻². The coast-wide average over the 10-year time span was 99.60 mg m⁻² (Table 2). There was high variability both spatially and temporally across the geographical range (Fig. 5d, S4). On average, chlorophyll-a concentrations increased over time across all regions, but with strong evidence that rates of increase varied spatially (Table 4c). The plot of back-transformed predicted means (Fig. 6c) showed that estimated rates of increases were much higher in Biological Region 4B than other regions, which had similar rates of increase.

Spatio-temporal modelling using environmental covariates

The fitted model with the lowest DIC included DO, but not temperature (Model G), with DO modelled as a quadratic function with separate intercept parameters when DO is less than 0.9 mL L⁻¹ (DIC = 9887.9, Table 5). The model without the additional DO intercept parameters (Model F) was only a slightly poorer fit, with DIC = 9888.3. Although a model with bottom temperature on its own provided a somewhat better fit than did the base model (DIC = 9957.8 for quadratic Model C v. 9961.6, Table 5), much greater improvements in the DIC were obtained by adding DO and removing temperature. Almost all temperature measurements were within the range typically encountered by Pacific halibut (98% of values were between 4 and 10°C; see Sadorus et al. 2014 for a review of temperature-range studies), and it seems likely that any apparent improvement in model fit when adding bottom temperature to the base model was a result of the strong correlation with DO (Fig. 3). Whereas all parameter estimates for bottom temperature had relatively large standard deviations and were greatly affected by the addition of DO to the model (even changing sign in Model E; Table 5), DO parameters have low standard deviations relative to the mean and were much less affected by whether temperature was included in the model or not. The k-fold cross-validation MSE values selected the same model as did the DIC, but ranked models including temperature more poorly than those that omitted this covariate.

Fig. 8 plots the observed WPUE data against DO, together with posterior predicted means of the 10 selected prediction stations (and 95% posterior credible intervals) for Model F. Although Model F had almost the same DIC as did the best-fitting model, it appears that a simple quadratic model over the whole range of DO values does not capture the drop to near-zero WPUE as DO decreases below 0.9 mL L⁻¹ (Fig. 8a). Likewise, when considering only the binary zero or non-zero model component (x(s,t)) in Fig. 8b, mean posterior predicted values of the probability of non-zero WPUE for the selected stations are much higher than the raw observed proportions (as calculated over all stations and years). Note that each point in Fig. 8b and 9b represents a proportion calculated from between 19 and 131 observations of zeros and non-zeros.

Fig. 8.

Comparison of model predictions (Model F) with raw data. (a) O32 WPUE plotted against DO, with points indicating raw observed data, and the line showing predicted values for the selected prediction stations (see text, Fig. 3), with 95% posterior credible intervals (shaded region). (b) Raw proportions of non-zero WPUE values from all stations in all years for binned DO values (points) and posterior model predictions of the probability of non-zero WPUE for the selected prediction stations in 2013 (solid line), with 95% posterior credible intervals (shaded region). Note that each raw proportion in (b) was computed from at least 19 and up to 131 binned observations. WPUE conversion: 1 lb = ~0.45 kg.

Fig. 9.

Comparison of model predictions (Model G) with raw data. (a) O32 WPUE plotted against DO, with points indicating raw observed data, and the line showing predicted values for the selected stations (see text, Fig. 3), with 95% posterior credible intervals (shaded region). (b) Raw proportions of non-zero WPUE values from all stations in all years for binned DO values (points) and posterior model predictions of the probability of non-zero WPUE for the selected stations in 2013 (solid line), with 95% posterior credible intervals (shaded region). Note that each raw proportion in (b) was computed from at least 19 and up to 131 binned observations. WPUE conversion: 1 lb = ~0.45 kg.

Model G, with separate intercept parameter for data with DO < 0.9 mL L⁻¹, does a little better at predicting the low probability of non-zero WPUE for very low DO (Fig. 9b); however, when considering WPUE overall (Fig. 9a), predictions are very close to the near-zero observed WPUE values associated with a very low DO. That this model had almost the same DIC as Model F is presumably due to the relatively low number of observations with DO < 0.9 mL L⁻¹, only 5% of the values used in this modelling. Overall, the quadratic shape of the prediction curves appears to characterise observed patterns in the data quite well. However, we note that we lack a physiological explanation for decreased Pacific halibut density at high values of DO. This highlights the need for more data at both high and low concentrations of DO, so as to provide more confidence in our understanding of the relationship at both ends of the range of DO values.