Representative range of acoustic point source measurements in the Chukchi Sea

Detecting and predicting biological responses to episodic perturbations (e.g., oil spills, marine heatwaves) and longer-term trends (e.g., global warming, ocean acidification) in marine environments requires a characterization and understanding of the ecosystem’s “natural” or baseline variability. However, characterizing biological variability is not easy because marine ecosystems are comprised of numerous physical and biological processes operating over multiple spatial and temporal scales (Stommel, 1963; Haury et al., 1978; Levin, 1992; Schneider, 1994).

Characterizing and monitoring biological variability in marine ecosystems requires high resolution, long-term (i.e., high scope) datasets that can be challenging to obtain due to limited resources or constrained accessibility. Such challenges are amplified in high latitude marine environments where the presence of sea ice during most of the year limits vessel-based sampling of the water column (e.g., Mueter et al., 2017; Spear et al., 2019). In these areas, data acquisition is typically limited in duration and/or resolution, fragmenting our understanding of important bio-physical processes over the annual cycle.

The use of echosounders attached to bottom-mounted platforms or moorings (i.e., stationary acoustics) provides a non-invasive technology to characterize and monitor temporal variability in abundance and behavior of pelagic organisms (e.g., Urmy et al., 2012; Horne and Jacques, 2018; Gonzalez et al., 2021). Stationary acoustics can: 1) characterize baseline conditions through continuous sampling of fish and zooplankton year-round; 2) quantify the amount of change relative to baseline variability to help determine thresholds of environmental effects; and 3) monitor changes, episodic or trending, in the abundance or behaviors of pelagic organisms. The inclusion of active acoustics in instrumented platforms (i.e., marine observing systems; e.g. Godø et al., 2014) also allows the simultaneous and continuous collection of biological and physical data, which can be used to understand underlying mechanisms of observed biological patterns, a prerequisite to predicting biological responses to environmental change (Linder et al., 2017).

Despite covering a wide spectrum of temporal scales (e.g., seconds to years), acoustic measurements from stationary echosounders are point source measurements and are limited in their spatial coverage. Because fish and zooplankton distributions are heterogeneous (i.e., patchy), biological similarity is expected to decay with geographical distance (Soininen et al., 2007). Beyond a certain distance, known as the representative range, meaningful inferences cannot be derived as uncertainty and interpolation errors are expected to increase (Martin et al., 2005; Anttila et al., 2008; Milewska and Hogg, 2010). When monitoring a large spatial domain, multiple point measurements within the representative range provide redundant information and a reallocation of sampling resources would lower sampling costs or could be used to expand the sampling domain. A better understanding of a site’s representative range ensures appropriate characterization and monitoring of the environment and, at the same time, optimizes the cost-effectiveness of monitoring through the deployment of one or multiple sensor packages.

Located over the Northeast Chukchi shelf on the southern flank of Hanna Shoal, the Chukchi Ecosystem Observatory (CEO; Figure 1) is a set of instrumented moorings that has been collecting high-resolution, continuous biological, biogeochemical, and physical measurements since 2014 (Danielson et al., 2017; Hauri et al., 2018; Lalande et al., 2020). The location of the CEO enables tracking of temporal variability and potential trends (i.e., ocean monitoring; Danielson et al., 2017) in a biological hotspot (Grebmeier et al., 2015). In this study we use concurrent mobile and stationary acoustic data to characterize spatial variability in fish and zooplankton densities and to quantify the representative range of CEO’s point source, active acoustic sampling. Comparison of stationary acoustic data with data from surface, mobile surveys will further the understanding of spatiotemporal biological variability in Arctic ecosystems and inform the design of effective monitoring networks to monitor and detect potential impacts of environmental change on pelagic communities.

The seasonally ice-covered NE Chukchi Sea receives a nearly continual input of heat, nutrients, organic carbon, and organisms from Pacific-origin water flowing northward in response to an oceanic pressure head that results from an elevation difference between the Pacific and Arctic Oceans (Stigebrandt, 1984) (Figure 1). This input from the Bering Sea, combined with shallow depths, enhances biological productivity in the Chukchi Sea (Grebmeier et al., 2015). A large phytoplankton bloom that occurs seasonally in late spring and summer (Questel et al., 2013) supports the largest soft bottom benthic faunal biomass in the world’s ocean (Grebmeier et al., 2006; Grebmeier et al., 2015), and corresponding populations of zooplankton (Ershova et al., 2015), seabirds (Kuletz et al., 2015), and marine mammals (Hannay et al., 2013). The Chukchi continental shelf is also home to Arctic cod (Boreogadus saida), a fish species that plays a key role in the transfer of energy from lower to higher trophic levels in high latitudes (Lowry and Frost, 1981; Whitehouse et al., 2014).

To characterize spatial variability in fish and zooplankton densities in the study area and to quantify the representative range of CEO point measurements, we analyzed active acoustic data from a stationary mooring (i.e., temporally indexed data) and coincident mobile surveys (i.e., spatially indexed data) conducted in proximity to the mooring (Figure 1).

Spatial variability in fish and zooplankton density was characterized using acoustic data from four transects (ML3, ML4, ML5, and ML6) measured during two Arctic Marine Biodiversity Observing Network (AMBON) surveys that included coverage of Hanna Shoal (Figure 1). These surveys were conducted during the periods of August 8 to September 5 in 2015 and August 5–24 in 2017. Spatially indexed acoustic data were collected using a Teledyne 307 kHz Workhorse Mariner acoustic Doppler current profiler (ADCP) hull-mounted to the survey vessel at a depth of approximately 3 m below the sea surface. Data were collected every 0.5 seconds (2 Hz) at a vertical resolution of 1 m. ADCPs can be used to measure acoustic backscatter (i.e., ensemble reflected energy) throughout the water column as a proxy for pelagic fish and zooplankton densities ( e.g., Cochrane et al., 1994; Ressler, 2002; Zedel et al., 2003) and have been used to describe diel migrations and distributions of biological scattering layers (e.g., Luo et al., 2000; Cisewski et al., 2010). ADCPs cannot be calibrated using standard reference targets (Demer et al., 2015) but comparisons of measurements by an ADCP to a calibrated scientific echosounder showed that data from the two instruments were highly correlated (Brierley et al., 1998). We assume that any bias in ADCP acoustic measurements is constant throughout the survey and that our measurements can be used to quantify relative change in acoustic backscatter.

All ADCP data were acquired using VMDas software (v. 1.46; RD Instruments, Inc.) and post-processed in WinADCP (v. 1.14; RD Instruments, Inc.) as part of standard quality control procedures. Data files were then exported and further processed in Matlab R2020a. Volume-backscattering strength (Sv measured in dB re 1 m^–1, hereafter dB) was calculated across a maximum of 25 vertical bins from the ADCP’s recorded raw echo intensity data using the sonar equation as described in Deines (1999) and updated by Mullison (2017). The backscatter sonar equation accounts for two-way signal transmission loss (i.e., time-varying gain) and sound absorption, and uses an instrument- and beam-specific returned signal strength indicator (RSSI) scaling factor to convert counts to decibels. This makes ADCP backscatter measurements more robust compared to raw echo intensity measurements, which are more readily extracted from an ADCP. Mean Sv was calculated for each vertical bin along each of the four beams of the ADCP. The mean of the four beams (Sv_mean) was used to calculate mean backscatter in the water column for each along-track 4-ping (i.e., 2 second) bin (corresponding to 8 to 10 m at vessel speeds of 8–10 knots). A double pass 1x3 median filter was applied to the Sv data to reduce contamination/interference from the ship’s echosounder.

Temporally indexed acoustic backscatter data was collected at the CEO (71° 35.976’ N, 161° 31.621’ W) using an upward-looking, Acoustic Zooplankton Fish Profiler (AZFP; ASL Environmental Sciences), deployed at a depth of 35 m (bottom depth of 45 m; Figure 1). The instrument operated at 38 (12°), 125 (8°), 200 (8°), and 455 (7°) kHz (nominal beamwidths, measured between half power points in parenthesis) since September 9, 2014. The AZFP collected data every 15 seconds (0.067 Hz) at a vertical resolution of 4 cm.

Acoustic backscatter from the AZFP was processed using Echoview software (v. 9.0). Background noise was subtracted and a minimum signal-to-noise ratio filter of 6 dB was applied. Echoes within 3 m from the face of the transducer were excluded from analyses to avoid integration of echoes in the acoustic nearfield (Foote, 2014). Seawater surface and sea ice edges were delimited using Echoview’s bottom detection and linear offset operator algorithm followed by visual inspection and manual correction. A surface exclusion line was set 0.5 m below the corrected surface and echoes above the line were excluded to ensure that backscatter from surface turbulence or sea ice was not included in analyses.

We used the 125 kHz AZFP mooring backscatter data corresponding to August and September of 2015 and 2017 to match the timing of mobile surveys in the area. The 125 kHz data were selected due to the presence of noise in the 200 kHz dataset. A threshold of –110 dB was applied and backscatter for the entire water column was integrated in 1-minute bins (i.e., 4 pings), matching the resolution of the spatial data.

Spatial and temporal data were collected during periods of midnight sun (i.e., 24 hours of daylight), and no significant biases are expected due to differences in availability of organisms to the acoustic beam throughout the day.

To characterize spatial heterogeneity and spatial dependence in acoustic backscatter we quantified and compared spatial variability in across-shore (ML3 and ML4) and along-shore (ML5 and ML6) transects within years, and from 2015 and 2017 AMBON surveys using autocorrelation functions (Legendre, 1993) and wavelets (Torrence and Compo, 1998).

Autocorrelation is a property of ecological variables that occurs when measurements that are close together in space or time are more similar than measurements taken farther apart. Spatial autocorrelation manifests as patches or gradients in biological distributions that can result from physical forcing of environmental variables or community processes (Legendre, 1993). Lagged autocorrelation functions quantify spatial dependency in acoustic backscatter by partitioning covariance over distance classes (Legendre and Fortin, 1989). Lagged Pearson’s product-moment correlation coefficients (ρ) define the correlation between all measurements at a given lag (h) and can be modeled using an exponential model:

where $ρˆ(0)$ is the autocorrelation at lag 0, d is the number of lags, and α is the range at which the autocorrelation decays by a value of e. We calculated the isotropic (i.e., direction-independent) and anisotropic (i.e., direction-dependent) correlograms of acoustic backscatter from transects conducted in 2015 and 2017.

Wavelet analysis can be used to explore spatial structure along one-dimensional transects (e.g., Torgersen et al., 2004; Grados et al., 2012; McGowan et al., 2019). A wavelet transform decomposes a data series as a function of scale and position through the convolution of a waveform—the wavelet—that is stretched or compressed (i.e., scaled) and slid through the series (i.e., translation). In this way, the wavelet decomposition identifies dominant spatial scales that account for variance throughout the series (Torrence and Compo, 1998). A continuous Morlet mother wavelet function (Torrence and Compo, 1998) was applied to each transect. Wavelet power was calculated using the R package WaveletComp (v. 1.1, Roesch & Schimidbauer, 2018). Statistical significance in localized wavelet power was evaluated through comparison to a white noise (constant value, equal to the time series variance) null hypothesis at a 95% confidence level (Torrence and Compo, 1998) using 100 simulations. Edge effects were minimized by adding zeroes at the beginning and end of each data series to increase the total length of the series to the next power of two (Torrence and Compo, 1998). Horizontal integration of wavelet power at each scale over the entire series—the global wavelet spectrum—allows the measurement of variance contributed by each scale across the entire transect. The global wavelet spectrum was calculated using the R package WaveletComp (v. 1.1, Roesch & Schimidbauer, 2018). Significance of this time-averaged variance was tested against white noise at a 95% confidence level (Torrence and Compo, 1998).

Six different methods described in Horne and Jacques (2018) were used to quantify the representative range of point measurements from the CEO AZFP. The use of multiple methods enables assessment of consistency among calculated ranges and facilitates generation of guidelines to design the spatial distribution of monitoring sensors. All methods calculate the representative range of the mean or variance of a quantity. Both mean and variance are considered appropriate monitoring metrics to track changes in biological variables (Underwood, 1991; Osenberg et al., 1994). The six methods can be categorized into four approaches: 1) distance between sensors based on spatial correlation; 2) sample size calculations assuming random sampling to detect a minimum threshold of change; 3) maximization of spatial variance; and 4) scales at which spatial and temporal variability are equivalent. The first approach calculates the optimum distance between sensors based on the relationship of spatial measurements. Using the decay of spatial correlation with distance (e.g., Anttila et al., 2008), the representative range corresponds to the distance at which measurements become independent. The second approach quantifies the number of sensors needed to detect change using a paired t-test, sample size calculation assuming a random sampling framework (e.g., Rycroft, 1949). Methods that use this random sampling approach include the number of replicates using a derivative of minimum sample size calculations for a paired t-test (Gray et al., 1992), a paired t-test, repeated measures ANOVA (Sullivan, 2006), and a sample size calculation for a paired t-test including statistical power (Zar, 2010). The third and fourth approaches quantify representative ranges of temporal variance rather than the mean. The third approach models the theoretical power spectrum as a function of the spatial autocorrelation (modeled in the first approach). The spatial period at which 95% of the maximum observed variance in fish density is set as the representative period of variance (e.g., Gilman et al., 1962). The final approach compares empirically derived spatial and temporal power-spectra to identify equivalent scales of spatial and temporal variability by identifying periods at which identical magnitudes of spatial and temporal variability are observed (e.g., Wiens, 1976).

The autocorrelation structure of each of the four transects surveyed in 2015 and 2017 is presented in Figure 2. The spatial distribution of vertically integrated backscatter corresponding to pelagic organism density, described by the distance at which samples become independent (ρ = 0) and the % of variance associated with spatial autocorrelation (ρ(1)), varied widely among transects. Differences between calculated ranges within across and along-shore transect pairs were typically larger than differences observed between transects of different orientation. Across-shore measurements became independent (ρ(2/√l)) at distances of the order of 15,000–53,000 m whereas along-shore measurements became independent at approximately 11,000–74,000 m. In 2015, more than half (58–89%) of the variability observed between sequential observations was explained by autocorrelation (ρ(1)) whereas autocorrelation explained a smaller proportion of this variability (30–38%) in 2017.

Scales of spatial variability differed among transects and between years (Figure 3; Figures S1 and S2). In general, spatial variability was concentrated at fewer scales in 2015 than in 2017. Peaks in averaged wavelet power were observed at scales ranging from 1,024 m to 73,500 m in 2015 and from 64 m to 73,500 m in 2017. In 2015, the largest peaks in averaged wavelet power (i.e., highest spatial variance) were observed at large scales (greater than 30 km) in ML3 and ML5 and at intermediate (10–30 km) scales in transects ML4 and ML6. In 2017, largest peaks in averaged wavelet power were observed at smaller scales (1–4 km) compared to 2015 for ML3 and ML5 and at large scales (greater than 40 km) for transects ML4 and ML6. These observations indicate variability in the spatial distribution patterns of organisms throughout the area and between years that could be associated with changes in the spatial structure of the environment (e.g., water mass distributions).

Spatial autocorrelation decayed exponentially as a function of distance (Figure 4). Exponential decay models fit to the lagged coefficient of determination (squared correlation coefficient, R²) followed y = 0.685e^{−0.0006456x} in 2015 and y = 0.229e^{−0.0007704x} in 2017. The coefficient of determination at lag 1 (ρ(1)) indicated that more than half of the variability between sequential observations was explained by autocorrelation in 2015 (ρ(1) = 0.68) whereas a smaller proportion of this variability could be attributed to autocorrelation in 2017 (ρ(1) = 0.23). The mean number of observations per transect (26,488 in 2015 and 25,179 in 2017) was used to calculate a 95% confidence interval for the lagged coefficients of determination. Based on calculated thresholds of significance (0.012 in 2015 and 0.013 in 2017) the representative range was calculated to be 49,817 m in 2015 (Figure 4a) and 30,101 m in 2017 (Figure 4b). The representative area, defined as a circle of radius equal to the representative range, was 7,797 km² in 2015 and 2,846 km² in 2017.

Representative ranges calculated using this method were highly dependent on the length of transects used in the analysis (Figure 5). To examine the relationship between representative range and transect length we calculated the representative range for each of the four transects surveyed in 2015 using transect lengths from 1.6 km (200 observations) to approximately 250 km in 0.8 km (100 observations) increments. Representative ranges increased with transect length with values ranging from 16 m to approximately 70,000 m and constituted 11–32% of the transect length on average. In general, representative ranges were similar for all transects up to a transect length of 18 km after which ranges increased at different rates. Representative ranges became asymptotic at different transect lengths for each transect (Figure 5). Minimal changes in representative range (i.e., rate change smaller than 0.5%) were observed for transect lengths greater than 100 km for ML3, 144 km for ML4, 149 km for ML5, and 165 km for ML6.

The number of required sensors to detect a minimum threshold of change was calculated for the 2015 and 2017 acoustic grids using the Gray et al. (1992) equation and derived paired t-test and power test equations. To meet independence requirements of sample size calculations, data were averaged into 49,817-m bins for the 2015 and 30,101-m bins for the 2017 data. These bin sizes corresponded to representative ranges calculated using the spatial autocorrelation method in Section 3.2.1. This coarser data resolution resulted in only 4–6 observations per transect. However, no significant differences in sample size calculations were observed when using at least 30 observations (approximately 375 m bins) per transect (results not shown).

Representative ranges varied among sample size calculation methods and between years (Table 1). The representative range using the Gray et al. (1992) equation was 12,271 m in 2015 and 19,325 m in 2017. Paired t-test calculations resulted in representative ranges of 305 m in 2015 and 756 m in 2017. Representative ranges calculated using the power test equation at a significance level of α = 0.05 and a probability of type II error β = 0.90 were 2,677 m in 2015 and 4,215 m in 2017. Because the choice of α and β are arbitrary, representative ranges for different combinations of α and β are presented in Figure 6. Setting α = 0.05 and the effect size E = 1 dB constant, the representative range increased 30% when decreasing beta from 0.90 to 0.70, from 2,677 m to 3,492 m in 2015 (Figure 6a), and from 4,215 m to 5,500 m in 2017 (Figure 6b). If beta is held constant at a conservative value of β = 0.90 and α is increased from 0.05 to 0.1, the representative range increases from 2,677 m to 2,965 m in 2015 (Figure 6a) and from 4,215 m to 4,670 m in 2017 (Figure 6b).

The representative range of backscatter variance corresponding to the scale at which 95% of the maximum variance (21.9 dB² in 2015 and 2.4 dB² in 2017) was quantified using the spatially indexed acoustic data. Differences of one order of magnitude were observed in representative ranges between years (Figure 7). The relative variance exceeded the 95% maximum threshold at a period of 2,524 m in 2015 (Figure 7a) and 259 m in 2017 (Figure 7b) and negligible increases in variance are expected beyond these distances.

Spatial and temporal spectral power increased with period (Figure 8). Units of space and time were set relative to the 8-m and 1-minute resolutions of the spatial and temporal datasets. In both years the spatial spectra (y₂₀₁₅ = –0.58 + 1.25x and y₂₀₁₇ = 0.30 + 0.99x) increased more rapidly than the temporal spectra (y₂₀₁₅ = 0.16 + 0.96x and y₂₀₁₇ = –0.197 + 0.92x). In 2015, the modeled maximum temporal variance at the Nyquist frequency, equivalent to a period of 21 days, was 4.46 log₁₀(dB²). The corresponding spatial variability was observed at a period of 86,152 m. Equal spatial and temporal variance was observed at a range of 1.5–3 units (250–8,000 m; Figure 8a). In 2017, the modeled maximum temporal variance, also corresponding to a temporal period of 21 days, was 3.92 log₁₀(dB²) and the corresponding spatial variability was observed at a period of 35,780 m (Figure 8b).

The spatial representative range of point source acoustic measurements at the CEO was dependent on the metric property (mean or variance), analytic method, and year (Table 1). The representative range of mean backscattering strength varied over one to two orders of magnitude among methods with values ranging from 305 to 49,817 m in 2015 and from 756 to 30,101 m in 2017. Differences of similar magnitude were observed for the representative range of backscattering strength variance with values ranging from 2,524 to 86,152 m and from 258 m and 35,780 m for 2015 and 2017 datasets, respectively. Large variations in representative range were also observed between years but these differences were smaller (two to tenfold differences between years) than those observed among methods.

Our approach comparing two data years and multiple methods to calculate spatial representative ranges at the CEO illustrates the occurrence of: 1) variations in calculated ranges between years (1–10-fold differences), 2) relatively large representative ranges (up to approximately 90 km), and 3) large variations in calculated ranges among methods (1–2 orders of magnitude differences). Each of these observations impacts the choice of representative range calculation method and the interpretation and application of results derived from point source observing systems such as the CEO in the Chukchi Sea.

Differences in representative range values between years are expected in response to changes in the timing and spatial distribution of environmental structure (e.g., sea ice melt fronts) and processes (e.g., water inflow from the Bering Sea) that determine water mass characteristics and associated species assemblages in the region (Spear et al., 2019). Interannual variations in sea ice extent and water inflow from the Bering Sea (Chen and Zhao, 2017; Zhang et al., 2020) result in interannual variations in water mass spatial distributions (Yang and Bai, 2020). Zooplankton assemblages in the Chukchi Sea can vary from predominance of small zooplankton species of Pacific origin (e.g., copepods Neocalanus spp., Eucalanus bungii) to large lipid-rich Arctic zooplankton species (e.g., Calanus hyperboreus) depending on interannual changes in advection pathways into the Arctic (Pinchuk and Eisner, 2017). Variations in water mass distributions and associated species compositions result in changes in acoustic backscatter spatial patterns that affect calculated representative range values. In the study area, differences were observed in current direction and water column characteristics (e.g., surface and near-bottom temperature and salinity, primary production; Danielson, 2021a, 2021b) between the two survey years, which could be contributing to differences in the spatial distribution of pelagic organisms and the resulting calculated representative ranges.

Relatively large representative ranges of acoustic measurements at the CEO, compared to similar studies, could be a result of low spatial heterogeneity in physical and/or biological characteristics, low abundance of organisms, and/or potential limitations of our sampling approach. Representative ranges calculated by Horne and Jacques (2018) for a stationary echosounder deployed in Admiralty Inlet, WA, were much smaller (maximum 1.4 km) than those observed in our study (maximum of approximately 90 km) but were calculated over a much smaller area (8 km²), in an environment with different physical characteristics and dynamics (temperate latitude, high tidal current speeds), and at a much higher spatial resolution (transects placed every 0.25–0.5 km). Large representative ranges observed in our study could be expected in environments with tight bio-physical associations and large characteristic environmental length scales; for instance, when spatial variability of physical characteristics that shape biological distributions is low, or in environments characterized by weak bio-physical associations when biological distributions are homogeneous and independent of the spatial structure of the physical environment. In the Chukchi Sea, zooplankton and fish distributions are tightly coupled with environmental features (e.g., water masses and frontal zones) that vary over scales consistent with our calculated representative ranges. Water masses in the study region often retain distinguishing characteristics over scales ranging from approximately 20 to 100 km (Day et al., 2013; Gong and Pickart, 2015; Danielson et al., 2020) and depend in part on the inflow of water masses from lower latitudes and the presence of sea ice and melt water in the summer (Spear et al., 2019). The length scale of lateral boundaries between oceanic water masses is set by the internal Rossby radius of deformation, which in the Chukchi Sea is often on the order of 1–5 km (Nurser and Bacon, 2014). Frontal zones in the Chukchi are found between water masses associated with sea ice melt, winter water, and summer shelf water (Lu et al., 2020). Zooplankton abundances, species compositions, and spatial distributions have all been shown to be influenced primarily by water mass properties (e.g., Hopcroft et al., 2010; Ershova et al., 2015; Pinchuk and Eisner 2017; Spear et al., 2019). Fish distributions in the Chukchi Sea are also influenced by physical environmental factors including bottom depth, water temperature, and salinity (Norcross et al., 2010; De Robertis et al., 2017; Iken et al., 2019; Forster et al., 2020). Iken et al. (2019) observed that distributions of demersal fish assemblages contained little spatial structure with variations in species composition and abundances observed mainly between coastal and offshore areas from samples collected during the 2015 AMBON survey. Day et al. (2013) reported variability in zooplankton and fish biomass across three regions located approximately 40 km apart in the NE Chukchi Sea. Large representative ranges could also be a result of weak spatial patterns due to low abundances and diversity of pelagic organisms in the benthic-dominated Chukchi Sea (Grebmeier et al., 2006).

Limitations of our sampling approach might also contribute to large representative ranges observed in our study. The lower sensitivity of ADCPs, compared to scientific echosounders, could result in weak biological spatial patterns that result in large representative ranges. In addition, acoustic backscatter measured using a 300-kHz ADCP includes small particle aggregates, bubbles, and electrical or mechanical noise that would reduce signal-to-noise ratios, and potentially obscure biological patterns. We observed strong biological scattering layers throughout all transects, and biological patterns are not expected to be affected by the sensitivity of the ADCP or by the inclusion of non-biological backscatter sources in our study. Finally, acoustic backscatter measurements integrate all sources of energy reflected from particles and animals in the water column. These values correspond to assemblages of zooplankton and fish species of different age/size classes where each category potentially has a different spatial distribution. As a result, strong heterogeneous spatial patterns of some species could be masked by homogeneous distributions of abundant species. A higher discrimination of integrated backscatter measurements into species and/or groups of species using calibrated, multi-frequency scientific echosounders could enable a more accurate characterization of spatial patterns and provide representative ranges specific for taxa of interest.

Differences in calculated representative ranges among methods are attributed to differences in the methodologies, rationale, and associated assumptions of each approach. Variations of 1–2 orders of magnitude in our calculations are consistent with differences of up to 2 orders of magnitude (approximately 30 to 1,400 m) observed by Horne and Jacques (2018) in a much smaller study domain. A comparison of methods to calculate representative ranges of air quality monitoring stations in EU countries also resulted in variations of 1–3 orders of magnitude in calculated ranges that were attributed to differences in the basic principles of the methods, in the definitions of similarity criteria and thresholds, and in the underlying definitions of representative ranges (Kracht et al., 2017). The choice of method will depend on the purpose of the study and the available data. Horne and Jacques (2018) provide a decision tree for method selection depending on the study/monitoring objective including detection of change (sample size calculation approach), mapping of spatial distributions (autocorrelation function approach), characterization of spatial variance (maximization of spatial variance approach), and interpolation of temporal variability over space (equivalent scales of spatial and temporal variability approach). This decision tree can be applied to any size sampling domain of interest.

Monitoring programs using stationary platforms in the Arctic will play a key role in the detection of biological responses to climate change, oil and gas exploration and exploitation, increased marine traffic and other ongoing environmental changes in high latitudes. Biological monitoring in remote areas needs to be efficient and cost-effective to be sustainable in the long term. Stationary sensors are assumed to be more cost-effective than repeated, mobile surveys and enable year-round sampling in areas with limited access due to seasonal sea ice cover. Quantifying the representative range removes uncertainty in the monitoring network design and reduces monitoring costs by identifying the minimum number of sensor packages for complete coverage of a site and provides a quantitative measure of the spatial scope of ongoing monitoring efforts. Stationary acoustic measurements at the CEO are representative of circular areas up to 24,000 km² (7,000 nmi²), depending on the calculation method, which correspond to areas of up to 46% of that covered by our spatial surveys. The study/monitoring goals, including the extent of the spatial domain, will impact the number of monitoring packages and their spatial distribution.

The calculation of representative range and design of a monitoring network system are best completed as part of a baseline characterization of spatial and temporal biological patterns in the study domain. All methods used in this study require synchronous spatial and temporal characterizations of biological distributions prior to the calculation of representative range(s). Spatial datasets used in this study come from surveys conducted in the vicinity of the CEO but were not explicitly designed for this study. Differences in equipment and frequency between spatial and temporal measurements could have introduced error in the calculation of representative ranges when using the equivalent spatial and temporal scales method. However, we expect differences among methods and between years to be larger than any differences associated with sampling. In our study, representative ranges were calculated based on spatiotemporal patterns of biological characteristics, and we did not consider physical features (e.g., depth, water mass characteristics) that may influence biological distributions. Also, best practice dictates that temporal and spatial data should be collected using instruments with the same characteristics. For acoustic data, calibrated, scientific echosounders operating at the same frequencies should be used for both temporal and spatial sampling. In all cases, a prior delineation of the study domain and a definition of what constitutes “similar” or “representative” should be aligned to meet the purpose of the study.

Although there has been considerable interest in understanding spatial and temporal biological patterns, little attention has been paid to the consequences of observed patterns for monitoring design (Rhodes and Jonzén, 2011). The concept of representative ranges was initiated and predominantly used in association with meteorological (e.g., Ciach and Krajewski, 2006; Milewska and Hogg, 2010) and air quality (e.g., Piersanti et al., 2015; Yatkin et al., 2020) monitoring. Attempts to evaluate and standardize methods to calculate representative ranges have been advocated in air quality measurements (Martín et al., 2015; Kracht et al., 2017). However, this renewed interest is not evident within the oceanographic community. Despite continued use of moorings and increased use of ocean observing systems, quantifying spatial representative ranges of stationary sensors is not commonly used to optimize sensor layout nor to quantify the spatial scope of monitoring measurements. Our results emphasize the importance of defining the study/monitoring spatial domain in conjunction with study objectives to ensure that representative samples are obtained for analyses. A calculation of representative ranges improves monitoring and characterization of biological communities and enables the design of cost-effective monitoring including the number, location, and spacing of sensors. Understanding how representative ranges of point measurements change depending on the timing of the survey, the location of stationary platforms, and the design of baseline spatial characterization (e.g., grid resolution and extent) is key to ensure the collection of representative samples and/or measurements. Research to understand how spatial representative ranges are affected by these factors could be achieved using simulated sampling schemes. Ensuring the collection of representative samples in a cost-efficient way becomes even more relevant in the face of environmental changes that are driving shifts in species distributions within high latitude environments.

ADCP data from the Arctic Marine Biodiversity Observing Network projects can be accessed through NOAA National Centers for Environmental Information (DOI: https://doi.org/10.25921/ew7d-qw23 and https://doi.org/10.25921/00zm-r593). Chukchi Ecosystem Observatory data can be found at https://aoos.org/project-page/ecosystems/chukchi-ecosystem-observatory/.

The supplemental files for this article can be found as follows:

Figure S1. Wavelet decomposition of ADCP-derived acoustic backscatter data from the 2015 AMBON cruise transects: (a) ML3, (b) ML4, (c) ML5, and (d) ML6.

Figure S2. Wavelet decomposition of ADCP-derived acoustic backscatter data from the 2017 AMBON cruise transects: (a) ML3, (b) ML4, (c) ML5, and (d) ML6.

The authors would like to acknowledge the captains and crews of the M/V Norseman II, R/V Sikuliaq, R/V Ocean Starr, and USCGC Healy for CEO mooring turnarounds, along with chief scientists C Ashjian, R Hopcroft, K Iken, R McCabe, and P Winsor.

SG is funded by NPRB Graduate Student Research Award, Oil Spill Research Institute Graduate Research Fellowship and Fulbright-ANII graduate fellowship program. This publication is partially funded by the Cooperative Institute for Climate, Ocean, & Ecosystem Studies (CIOCES) under NOAA Cooperative Agreement NA20OAR4320271. SLD receives funding for the CEO from AOOS grants G9046 and G11133 and NPRB projects #1426 and L36-00A (NPRB publication # 1901). LL is funded by The Bryden Centre project which is supported by the European Union’s INTERREG VA Programme, managed by the Special EU Programmes Body (SEUPB).

The authors declare no competing interests.

Main concept and key ideas: SG, JKH.

Collection and compilation of observational data: SG, SLD.

Analysis and interpretation of data: SG, JKH, SLD, LL, GL.

Drafting of article and contribution of intellectual content: SG, JKH, SLD, LL, GL.

Final approval of the version to be published: SG, JKH, SLD, LL, GL.

How to cite this article: González, S, Horne, JK, Danielson, SL, Lieber, L, López, G. 2022. Representative range of acoustic point source measurements in the Chukchi Sea. Elementa: Science of the Anthropocene 10(1). DOI: https://doi.org/10.1525/elementa.2021.00055

Domain Editor-in-Chief: Jody W. Deming, University of Washington, Seattle, WA, USA

Associate Editor: Julie E. Keister, Biological Oceanography, University of Washington, Seattle, WA, USA

Knowledge Domain: Ocean Science