Impacts of infauna, worm tubes, and shell hash on sediment acoustic variability and deviation from the viscous grain shearing model

Underwater acoustics in littoral and bottom-interacting ocean environments is strongly controlled by seabed properties, and high-resolution geoacoustic data are needed to accurately predict acoustic propagation and scattering in such environments.¹ Furthermore, seabed geoacoustic variability introduced by geological, hydrodynamic, and biological processes can occur over a wide range of spatial scales ranging from fine-scale variability on the order of meters or less to large-scale variability over kilometers.² Understanding seabed geoacoustic variability at multiple spatial scales and over a wide frequency band is critical for the accurate prediction of underwater acoustic propagation and unbiased estimation of seabed properties in heterogeneous environments using remote sensing methods.^3,4

Marine sediments are complex physical and biological systems. Bottom types can range from mud to sand to rock, with various combinations of geologic and geomorphologic features (e.g., smooth, rippled, or jagged). Sediment grain size distributions can be very narrow, as commonly seen with well sorted sand, or be inclusive of diverse particle sizes, as in heterogeneous muddy sediments which can include clay, silt, and sand in varying proportions. Additionally, organic matter can have a strong influence on sediment cohesive properties, particularly for fine-grained sediments,⁵ and infaunal organisms can contribute to horizontal and vertical variability near the seabed surface, known as the sediment-water interface.⁶ 

Infauna are marine organisms that live within sediments, generally within the top 10 cm of the surface.⁷ Depending on the community of organisms present, infauna can re-work the sediment in different ways, altering sediment physical properties, contributing to surface and volume scattering, and affecting sediment cohesion.^2,8–10 Bioturbation (the mixing of sediments) and other biogenic processes can influence sediment physical properties, such as porosity, permeability, and the presence of gas due to organic decomposition by microbial processes.² Infauna can affect grain size distributions through the production and distribution of fecal pellets, which can be deposited on the sediment surface or in subsurface sediments and may be transported through bioturbation.¹¹ Suspension feeders, including bivalves and polychaetes, can capture fine suspended particles, aggregate them, and deposit them as fecal pellets on the sediment surface. Bioirrigation of burrows transports oxygen to anoxic subsurface sediments, creating sharp redox gradients and modifying sediment chemistry.¹² Increased cohesion by exopolymeric substances (EPS) secreted by microbes and microalgae can contribute to enhanced shear strength near the sediment surface;¹³ however, infaunal organisms can substantially reduce the stabilizing effects of EPS through bioturbation.¹⁴ These combined effects can manifest as temporal and spatial variability in compressional and shear wave propagation parameters in the seabed surficial sediment,^15–19 which can impact acoustic scattering from and propagation into the seabed.^20–28 Further investigation is needed to determine the effects of infauna on sound propagation in ocean bottom sediments, which could better inform sonar operation in various shallow-water applications, as well as the use of acoustics for remote sensing of benthic ecosystems.^29–31 

A deficiency that exists in our present understanding of sediment acoustics is that predictive models do not explicitly account for the physical impacts of volume heterogeneities (e.g., infauna bodies or burrows, worm tubes, or shell hash) on geoacoustic properties.^32–34 Part of this stems from a paucity of broadband geoacoustic data that can be directly compared with data on physical sediment properties from cores and adequate characterization of the distribution of infauna and other potential heterogeneities within the upper portion of the seabed. A goal of this work is to provide a new data set that facilitates such a comparison between acoustic variability and data-model deviation with sediment biological properties.

The main hypothesis of this paper is that seabed biological activity (1) increases variability in acoustic properties over small spatial scales and (2) increases the deviations from sediment acoustics model predictions. This paper reports on a field experiment that was conducted to test the prediction that (a) infaunal presence and activities, (b) worm tubes, and (c) shell hash have the potential to affect sediment acoustic properties. In Sec. II, the field experiment, in which in situ sediment acoustics measurements and diver cores were collected in natural marine sediments, is described. The viscous grain-shearing (VGS) sediment acoustics model^33,35–40 is summarized in Sec. III. Core analysis results describing vertical profiles of sediment bulk properties (density and porosity) and grain size distribution are presented in Sec. IV, in addition to measurements characterizing the vertical distributions of infauna, worm tubes, and shell hash. The measured sediment acoustic properties and data-model comparisons are presented in Sec. V, and further analysis of the data-model deviation in the context of the infauna, worm tube, and shell hash distributions is explored in Sec. VI. Finally, the results are discussed in Sec. VII, and conclusions are summarized in Sec. VIII.

A field experiment was conducted during May 8–11, 2017, in Petit Bois Pass, a channel that lies approximately ten miles off the northern Gulf of Mexico coast and in between the barrier islands of Petit Bois Island, Mississippi to the west and Dauphin Island, Alabama to the east (Fig. 1). A seabed lander with acoustic probes (Fig. 2) was deployed from the vessel (the R/V E. O. Wilson) to the seabed surface to collect in situ measurements of sediment sound speed, sound attenuation, and shear speed within the top 20 cm of the sediment, spanning the depth range within the sediment in which infaunal activities are expected, generally in approximately the top 10 cm.⁷ Diver cores were collected from each measurement location to sample the sediment's physical properties (bulk density, porosity, and grain size distribution), the local infauna distribution, and the biogenic features of shell hash and worm tubes. Prior to destructive sampling, the cores were vertically scanned with an acoustic core logger to provide additional sound speed and attenuation measurements. Vertical gradients were examined within two sites, as infauna are most abundant near the sediment-water interface, and the experiment targeted sites with different communities, specifically higher and lower abundances of tube-building worms, which have been shown to increase attenuation.¹⁹ 

The in situ measurements and cores were collected at two sites that differed in infaunal biomass, number of worm tubes, and amount of shell hash. The water depth at both sites was approximately 6 m. Site 1 (⁠$30°$ 14.754′ N, $88°$ 22.590′ W) was near the middle of the pass, and Site 2 (⁠$30°$ 13.782′ N, $88°$ 24.300′ W) was 3.3 km to the west near Petit Bois Island. Bottom water sound speeds calculated from a conductivity-temperature-depth (CTD) probe were 1518 m s⁻¹ and 1521 m s⁻¹ at Sites 1 and 2, respectively. Three separate deployments of the in situ acoustic measurement system were conducted at each site within approximately 3 m of each other to sample intra-site variability in sediment acoustic properties.

The in situ acoustic measurement system employed two sets of probes to measure compressional (60–100 kHz) and shear (0.2–0.9 kHz) waves in the sediment. Compressional waves were generated and received by custom-built, spike-like probes with pointed tips that housed cylindrical lead zirconate titanate (PZT) transducers [Fig. 2(a)], and circular-disk, bimorph piezoelectric benders housed in custom-built, spear-like probes [Fig. 2(b)] transmitted and received shear waves.^41–43 The lander employed to deploy the probes [Fig. 2(c)] was lowered from the ship's A-frame at each measurement location. Additional details of the measurements and data analysis used to extract sound speed and attenuation from the compressional wave measurements and shear speed from the shear wave measurements are given in Ref. 44.

Several cores were collected by scuba divers at each in situ measurement location so that sediment physical properties, shell hash, worm tubes, and infauna distribution could be characterized. Three 15.2-cm-diameter cores and one 7.6-cm-diameter core were collected from in between the in situ acoustic probes at each deployment location. The cores were capped on both ends underwater before being brought to the surface by the divers. Two of the 15.2-cm-diameter cores from each deployment location were cut into sections aboard the research vessel and sieved (500-μm mesh) for infauna, worm tubes, and shell fragments, which were subsequently preserved in 95% ethanol and stored for later analysis and classification. The remaining cores from each measurement location were stored vertically in an aerated seawater tank for transport back to the laboratory for further analysis and processing.

Diver cores returned to the onshore laboratory underwent non-invasive acoustic measurements using a custom-built, vertically oriented, broadband acoustic core logger called CARL⁴⁵ before any destructive measurements were performed on the cores (Fig. 3). CARL employed two modes of operation: (i) a high-frequency (100–300 kHz) transmission mode for depth-resolved measurements of travel time and amplitude along the length of the core to provide vertical profiles of sound speed and attenuation and (ii) a low-frequency (8–30 kHz) resonance mode for determining the effective sound speed of the entire core. Details of the core logger operation and data analysis are given in Ref. 44.

After the CARL measurements were completed, all 7.6-cm-diameter cores (three per site, one per deployment location) were sectioned and analyzed for vertical profiles of sediment bulk density, porosity, and grain size distribution. The sieved material from the 15.2-cm-diameter cores (nine per site) was quantified by mass as well as abundance or number of animals to account for differences in size and density of infauna, hard worm tubes, and shell hash. These quantities are reported both in terms of depth-integrated biomass per area for overall site comparisons and depth-separated biomass per core section volume for comparison with the measured acoustic gradients. The average worm tube diameter and shell fragment size are also reported.

In addition to considering total infauna biomass, infauna were sorted into functional groups based on characteristics such as body type (hard or soft), sediment activity (mixing or structuring), and type of worm tube (soft, hard, or no worm tube). The four functional groups identified in this work are Hard-bodied, No tube, Mixing (HNM), consisting of animals like brittle stars or mollusks that scatter sound and increase attenuation in the same way as other physical inhomogeneities, e.g., shell fragments of similar size; Soft-bodied, No tube, Mixing (SNM), which includes various burrowing worms that modify bulk sediment properties by dilation or compaction; Soft-bodied, Hard-tube, Structuring (SHS) worms that build tubes consisting of sediment, shell hash, and mucus that may increase attenuation due to the scattering of acoustic energy;¹⁹ and Soft-bodied, Soft-tube, Structuring (SSS) worms that build tubes of fine-grained sediment and mucus that may increase shear strength but are unlikely to substantially increase scattering. Additional details about the diver core analyses are provided in supplementary material.⁴⁴ 

It is useful to compare the acoustic data with the predictions of an established sediment acoustics model that relates wave speed and attenuation dispersion to sediment physical properties that can be obtained from the cores. The VGS model treats unconsolidated marine sediment as a dissipative two-phase continuum capable of supporting elastic wave propagation.^36,37 Internal losses are attributed to shearing and stress relaxation at grain-to-grain contacts and viscous dissipation in the thin layer of pore fluid between the grains. The VGS dispersion relations have been shown to agree with various data sets in sediments ranging from coarse-grained sand to fine-grained mud consisting predominantly of silt and clay.^33,40 The VGS model was chosen for this comparison because of its generality and because it could be efficiently implemented with very few fitting parameters. The model is briefly described here; however, for detailed derivation and description of the model, the reader is referred to the original literature.^33,35–40

The VGS expressions for sound speed c_p, attenuation α_p, shear speed c_s, and shear attenuation α_s are given by

and

where $i=−1$⁠.

Familiar physical parameters in the above expressions are the sediment bulk density ρ, the angular acoustic frequency $ω=2πf$⁠, and the sediment sound speed c₀ in the absence of inter-granular friction, i.e., the fluid limit, given by the Mallock-Wood equation,^46,47

Here, B is the effective bulk modulus of the sediment,

where B_f is the pore fluid bulk modulus, B_s is the mineral grain bulk modulus, and β is the fractional porosity. The bulk density is related to the porosity and the pore fluid and grain densities ρ_f and ρ_s by the mixture rule,

The functions $gp,s(ω)$ represent the effects of viscous dissipation,

where $τp,s$ are the viscoelastic time constants for the compressional and shear waves, respectively. The compressional and shear viscoelastic time constants are related by the expression³⁸ 

The remaining VGS parameters related to the stresses and strains associated with inter-granular motion are the compressional and shear moduli, γ_p and γ_s, and the strain-hardening index n. The arbitrary time normalization constant T was introduced by Buckingham in his original derivation of the grain-shearing theory to keep quantities raised to the fractional power of n dimensionless; however, he showed that the expressions for sound and shear speed and attenuation are independent of the value of T, which was taken to be T = 1 s for convenience.³⁶ 

The VGS model employs a sub-model that considers the effects of overburden pressure on the grain-to-grain contacts, which imparts dependence of compressional and shear moduli on depth in the sediment d, grain size u, porosity β, and the elastic properties of the mineral grains.^33,40 In the case where porosity, grain size, and elastic properties of the grains are homogeneous with depth, the expressions for γ_p and γ_s are given by

and

where the scaling coefficients $γp0$ and $γs0$ and the normalization parameters β₀, d₀, u₀, and H₀ are constants, and

where E_s and σ_s are the grain Young's modulus and Poisson's ratio, respectively. In the case that the sediment has M_T layers, the compressional and shear moduli in layer M are given by⁴⁰ 

and

Here, $ρs0$ and $ρf0$ are additional normalizing constants, h_m is the thickness of the mth layer, and all other parameters are allowed to vary from layer to layer.

It is important to note that the VGS model does not explicitly include the effects of shell hash, worm tubes, infauna, or other inhomogeneities in the medium, affecting the accuracy of the model predictions when such inclusions are present. This has been noted previously, particularly in regard to comparison of sound attenuation predictions with previous sediment acoustics measurements.³³ Furthermore, the model does not account for grain-size heterogeneity but assumes a single grain size u. In practice, the mean grain size is often used for this model input parameter, as is done here, for sediments with a distribution of grain sizes.

Comparison of depth profiles of the sediment bulk density, porosity, and grain size distribution from each core collection location within both sites demonstrates the level of intra-site variability in sediment properties (Fig. 4). For each of these parameters, the measured profiles within a single site tended to be consistent with each other, and there were only a handful of outliers present.

In Site 1, the bulk density was highest near the sediment-water interface and decreased with depth into the sediment. In contrast, the bulk density profile at Site 2 was fairly uniform with depth in the top 20 cm. At Site 1, the porosity had a positive gradient with depth, ranging from approximately 50% near the seabed surface to 70% at 20-cm depth. In contrast, the porosity profile at Site 2 was nearly constant, with a mean value near 60% throughout top 20 cm, with only a slight negative gradient.

Both sites were predominantly sand and silt. For each site, the fraction of sand was the greatest near the seabed surface, about 80% sand for Site 1 and 60% sand for Site 2. As sediment depth increased, there was lower sand content (approximately 40%–50% at 20 cm) and higher silt content at both sites; however, the gradients were steeper at Site 1. Gravel and clay fractions were low at both sites. Shell hash was included in the grain size distribution as part of the gravel fraction, but it accounted for less than 1%–2% of the sediment content by mass. Based on the Folk's sand-silt-clay ratio ternary diagram,⁴⁸ the sediment at Site 1 was silty sand (50% < sand fraction < 90%; silt-clay ratio > 2:1) transitioning into a sandy silt (10% < sand fraction < 50%; silt-clay ratio > 2:1) as the sediment depth increases, and Site 2 was comprised of sandy silt sediment that was more uniform with depth.

Mean grain size is given both in units of $ϕ$ (see Ref. 44) to aid in comparison with sediment physical and acoustic properties compiled in the literature² and in micrometers because the VGS model uses SI units. Based on the Udden-Wentworth scale,^49–51 the mean grain size at Site 1 ranged from very fine sand near the surface (⁠$3<ϕ<4$⁠, or between 125 and 62.5 μm) to coarse silt near 20-cm depth (⁠$4<ϕ<6$⁠, or between 62.5 and 15.6 μm). Site 2 was coarse silt, with mean grain size in the 4–6-$ϕ$ range. The larger mean grain size at Site 1 (near the middle of Petit Bois Pass) was potentially influenced by the hydrodynamic environment. Site 2 was near the western edge of the pass and potentially less affected by current, and thereby better able to retain finer grains near the seabed surface. For the grain size distributions here, the inclusive standard deviation was $1<σI<2.5$⁠, indicating poorly to very poorly sorted grain size distributions based on the descriptive terminology of Folk and Ward,⁵² in contrast to sediment acoustic model assumptions of a well sorted distribution characterized by the mean grain size.

Because the distributions of the depth-integrated mass per surface area of infauna, worm tubes, and shell hash from the cores at each site tend to have a great deal of spread and contain several outliers, the results are reported in terms of median and interquartile range instead of mean and standard deviation (Table I). For comparison, the depth-integrated mass per surface area was estimated from cores that were processed for sediment properties. Infauna, hard worm tubes, and shell hash represented approximately 0.02%, 0.03%, and 0.3% of the mass per surface area of the overall seabed material, respectively.

The total infauna biomass per site was similar, with the Site 1 median value being approximately 1.5 times greater than Site 2. Most of the animal biomass at both sites fell into the hard-bodied, no tube, mixing (HNM) functional group, which includes larger hard-bodied or hard-shelled animals such as brittle stars, crustaceans, and mollusks that burrow in the sediment. The body sizes of many of the HNM animals were on the order of 1 cm or greater, thus they are potential acoustic scatterers at the upper end of the acoustic frequencies used in this experiment.

Another functional group of potential acoustic significance at these sites was the soft-bodied, hard-tube, structuring (SHS) group, which included hard-tube-building worms such as Owenia and Diopatra. These animals structure the sediment by building tubes made from shell fragments, sediment, and mucus. Although the biomass of SHS was lower than that of the hard-bodied group, the mass of tubes was high, especially at Site 1, which had an approximately six times greater median hard worm tube mass than Site 2, primarily from the polychaete, Owenia fusiformis (Table I). The diameter of the Owenia tubes ranged between 2 and 5 mm (2.6 ± 1.0 mm at Site 1 and 2.8 ± 1.6 mm at Site 2, mean ± standard deviation), but the tube can extend from above the surface down to 5–10 cm depth; therefore, they are also good candidate acoustic scatterers.

The soft-bodied, soft tube, structuring (SSS) and soft-bodied, no tube, mixing (SNM) functional groups were also present, although in small biomass quantities that were similar between the sites. Thus, it is assumed that these animals had less potential influence on the sediment and, hence, the acoustics at these sites.

Site 1 had roughly twice as much mass per surface area of shell fragments present compared to Site 2 (Table I). Site 1 also had larger fragments of shell hash, 3.5 ± 1.6 mm at Site 1 compared to 2.3 ± 1.3 mm at Site 2 (mean ± standard deviation).

Vertical profiles of infauna biomass density $ninf$ (Fig. 5) displayed negative vertical gradients at both sites, where most of the infauna resided in the upper 7.5 cm of the sediment, 72% at Site 1 and 69% at Site 2. Thus, the infauna biomass density in the upper portion of the sediment was comparable between the two sites. Greater differences were seen between the two sites in terms of hard worm tube mass density $ntube$ vertical profiles. The mean value of $ntube$ in the upper 7.5 cm of the sediment was approximately eight times higher at Site 1 than at Site 2. The shell hash mass density $nsh$ profiles displayed less steep vertical gradients than infauna or worm tube mass density, thus they were relatively constant over the upper 20 cm (Fig. 5). For the most part, there was approximately twice as much shell hash at Site 1 than Site 2 over all depths, with the exception of one high shell hash mass density value observed close to 20 cm from one core at Site 2.

Overall, Site 1 had more hard worm tubes, SHS functional group infauna, and shell hash than Site 2, but two sites had similar overall amounts of infauna. Both sites had similar depth-dependence of infauna (negative gradient) and shell hash (no gradient). Finally, Site 1 had a higher mass of hard worm tubes near the surface than Site 2.

Acoustic results are presented in this section from the in situ field measurements and the laboratory CARL measurements. By combining the in situ and CARL data sets, the acoustic behavior of the sediment can be examined over many decades of frequency, between 8 and 300 kHz for sound speed ratio, 60–300 kHz for attenuation, and 0.2–0.9 kHz for shear speed.

Rather than simply presenting the measurements by themselves, they are compared with the VGS model to place the acoustic measurements in context with the measured profiles of sediment bulk properties and grain size distribution. Therefore, the selection of model inputs is first discussed. We then present broadband comparison of the acoustic data with the VGS model, and we examine vertical profiles of sound speed ratio, attenuation, and shear speed.

As implemented here, the VGS model equations [Eqs. (1)–(14)] are expressed in terms of a set of eight measurable physical parameters (β, ρ_f, ρ_s, B_f, B_s, H, d, and u), ten fixed constants (⁠$γp0, γs0$⁠, β₀, d₀, u₀, H₀, $ρs0$⁠, ρ₀, n, and T), and the two viscoelastic time constants $τp,s$⁠. All of the measurable parameters can be obtained from either the literature or analysis of the sediment cores, pore fluid, or bottom water properties. Tabulated values for pore fluid density ρ_f, grain bulk modulus B_s, grain density ρ_s, and H were taken from the literature.^2,40 Pore fluid bulk modulus was estimated using the pore fluid density and mean measured bottom water sound speed, $Bf=ρfcf2$⁠. For depth-dependent data-model comparisons, vertical profiles of porosity β and mean grain size u averaged over the various core measurements shown in Fig. 4 were directly input into the VGS model. Although the grain size fraction varied with depth in the sediment, the predominate grain size classes were sand and silt; therefore, single values of ρ_s, B_s, and H for quartz were used. The fixed constants were given in previously published literature⁴⁰ and are listed in Table II, and the viscoelastic time constants were determined by a fit to the broadband sound speed ratio data.

The compressional viscoelastic time constant τ_p was obtained by fitting Eq. (1) to the broadband (8–300 kHz) combined in situ and CARL sound speed ratio data sets. For this fit, the model was implemented with a homogeneous sediment with compressional and shear moduli given by Eqs. (10) and (11). Furthermore, the depth-dependence of the in situ and high-frequency CARL pitch-catch data were not taken into account for this fit because the low-frequency CARL resonator measurements return effective sound speeds for the entire core and are inherently depth-independent in the 20-cm measurement interval. Therefore, to find the best-fit τ_p, the in situ and high-frequency CARL pitch-catch measurements were depth-averaged at each frequency before being input into the fitting routine. Best fits between the modeled and measured sound speed ratio were obtained using an unconstrained nonlinear curve-fitting algorithm to minimize the objective function

where $M(τ′,fj)$ is the modeled sound speed at the jth frequency for the value $τ′$ of the viscoelastic time constant and $c( fj)$ is the measured sound speed. After the best-fit value of τ_p was found for each site, computed high-frequency values of sound speed and shear speed well above the compressional threshold frequency (⁠$ft,p=1/2πτp$⁠) were used to calculate τ_s from Eq. (9).

The best-fit values of τ_p and τ_s for each site are listed in Table III, and these parameters were then combined with vertical profiles of porosity and mean grain size from the cores as inputs to the VGS model. Average vertical profiles were computed from the porosity and mean grain size data shown in Fig. 4 and were used as direct inputs to the VGS model to produce depth-dependent predictions of c_p, α_p, and c_s for comparison with the acoustic data. To calculate the depth-dependent compressional and shear moduli, each core section was treated as a sediment layer, and the porosity, grain size, and thickness of each section were input directly into Eqs. (13) and (14). The depth-dependent VGS models were then computed at the same frequencies as the in situ and high-frequency CARL measurements for direct comparison.

The low-frequency asymptote of the VGS sound speed prediction is the Mallock-Wood sound speed, given by Eq. (5), depth-averaged values of which are within 1% of each other for the two sites (Fig. 6). The VGS model predicts a transition to higher sound speed ratio between low- and high-frequency, and based on the best-fit values of τ_p, the transition between low- and high-frequency regimes (at the threshold frequency $ft,p=1/2πτp$⁠) occurs approximately at 27 ± 3 kHz for Site 1 and 26 ± 2.7 kHz for Site 2. The VGS compressional and shear moduli [Eqs. (10) and (11)] both increase with increasing depth, increasing mean grain size (u), and decreasing porosity, leading to a complex effective dependence of sound speed on depth in the presence of both porosity and mean grain size gradients. At Site 1, the mean grain size decreased and the porosity increased with increasing depth, leading to a predicted sound speed ratio profile with a negative gradient (i.e., sound speed decreases with depth). Site 2 had more uniform porosity and mean grain size profiles, thus the predicted sound speed increases with depth (a positive sound speed gradient).

The low-frequency sound speed data from CARL's resonance mode (Fig. 6) represent the effective sound speed of the entire core averaged over the entire 20-cm measurement depth interval. The measured low-frequency sound speed ratios tended to be less than unity, indicating that at low frequencies the sediment was acoustically similar to mud at both sites, and the low-frequency CARL measurements approached the depth-averaged VGS low-frequency limits for each site. Similar to the VGS prediction, there was an increase in sound speed between the low-frequency CARL measurements and the higher frequency in situ measurements, and sound speed variability was greatest in the high-frequency CARL measurements. The high-frequency sound speed values at Site 1 were approximately 2.5% higher than those at Site 2. The deeper high-frequency CARL measurements from Site 1 agree better with the deeper VGS-predicted sound speed ratios, but they exceed the model predictions at shallower depths. Scatter in the high-frequency sound speed data at Site 2 encompasses the range of depth-dependent VGS predictions; however, there is little correlation between the shallower measurements and predictions.

Site 1 sound speed ratio decreased with sediment depth (Fig. 7), which was more evident in the CARL measurements because the in situ measurements had much coarser depth-resolution and did not extend closer to the interface than 5 cm. As noted previously, this decrease in sound speed ratio with depth was in part due to the increase in porosity and decrease in mean grain size with depth at this site. In situ sound speed ratio measurements at Site 1 tend to be higher than the VGS model predictions by approximately 2%–5%, and this difference is greatest at the lowest frequencies where the in situ measurements partially overlapped the transition region predicted by the VGS model. Site 1 high-frequency CARL measurements in the 4–20 cm depth interval are in good agreement with the predicted sound speed ratio profile; however, the measurements deviate from the predicted sound speed ratio in the 1–4 cm interval closest to the seabed surface, where the measured sound speed ratio gradient steepens and the variability in the measurements increases. In contrast, both the in situ and CARL sound speed ratio measurements are in better agreement with the VGS model at Site 2; however, variability in the measurements is greatest at this site in the 1–4 cm interval CARL measurements, as is the discrepancy between the data and the model. Finally, both in situ and high-frequency CARL data sets include measurements at a frequency of 100 kHz, and it should be noted that the data points from these different measurement systems overlap at this frequency.

Below the threshold frequency $ft,p$⁠, the VGS model predicts that the attenuation scales with frequency squared, whereas above the threshold frequency the attenuation is expected to scale approximately linearly with frequency. The lowest frequency in situ attenuation measurements (60 kHz) are just above the expected threshold so no data are available for comparison with the predicted sub-threshold frequency dependence (Fig. 6). The attenuation data between 60 and 140 kHz roughly agree with the VGS prediction, although a good deal of variability is present in the data, which is not captured by the depth-dependent model predictions. Above 180 kHz, the attenuation data appear to depart from the approximately linear frequency dependence predicted by the VGS model. Indeed, comparing the high-frequency CARL data above 180 kHz to a power law fit, $αp∼fk$⁠, we find that k = 1.31 for Site 1 and k = 1.23 for Site 2, revealing a departure from linear frequency dependence.

Based on the porosity and grain size vertical profile measurements from the cores, the VGS model predicts increasing attenuation with depth at both sites, in which the attenuation at Site 1 is expected to double with increasing depth in the 1–20 cm interval and increase by a factor of 2.3 in the same interval at Site 2 (Fig. 8). Note that the measured and predicted gradients do not change much across the 60–300 kHz band for both sites, whereas the measurement-model deviation increases with frequency. The in situ attenuation measurements at both sites are over-predicted by the VGS model by factors between 1.4 and 2.7, with larger data-model deviation observed at the lowest frequencies. High-frequency CARL measurements at both sites are in good agreement between 100 and 140 kHz below 4 cm, although the model tended to under-predict the data above 4 cm. Like the sound speed profiles, the variability in the attenuation measurements is highest in the 1–4 cm interval. As noted, before, the high-frequency CARL attenuation measurements are under-predicted by the VGS model in the 180–300 kHz band, and this disagreement occurs over a greater portion of the measurement depth-interval as the frequency increases (Fig. 8). This is most apparent at 300 kHz, where there is higher attenuation measured than predicted for most of the measurement depth interval. There is greater deviation between the attenuation measurements and model at Site 1, particularly at higher frequencies, which could be influenced by the greater abundance of hard worm tubes at Site 1, which might have more of an effect as frequency increases.

Values of shear speed measured at both sites range between 10 and 40 m s⁻¹ (Figs. 6 and 9). The VGS model shear speed prediction overlaps with the measurements at the lower end of the 0.2–0.9 kHz band; however, the measured shear speed displays little variation with either depth or frequency, which is in disagreement with the model. Note that the shear threshold frequency, $ft,s=1/2πτs$⁠, is 1.7 kHz for Site 1 and 1.5 kHz for Site 2, placing the measurement band in the transition regime where greater dispersion would be expected. If the shear viscoelastic time constant is set to $τs=∞$⁠, however, the resulting VGS model prediction (dashed lines in Fig. 6) exceeds the data by more than a factor of two with no overlap between the model and data. It is worth noting that in other recent applications of the VGS model, the shear viscoelastic time constant is taken to either be infinite^39,40 or orders magnitude larger than the best fit values here;^53,54 however, there have been very little shear wave data in those applications to justify otherwise. The predicted shear speed gradient displays an approximately cubic-root power-law depth dependence (Fig. 9). The in situ shear speed data tend to cluster around the mean measured shear speed of 25.3 m s⁻¹, which is close to the predicted shear speed at 0.2 kHz at both sites; however, the data-model deviation increases as the frequency increases for the entire measurement depth interval. The measured shear speeds tend to skew roughly 50% lower than the best-fit model predicts on average, indicating the sediment is slightly softer than is predicted by the VGS model based on the porosity and mean grain size inputs from the cores.

In this section, relationships between the deviation from the acoustic data to the VGS model and infauna, worm tube, and shell hash mass density profiles are examined to test the hypothesis that the infaunal community alters sediment acoustic properties. Although correlations between the data-model deviation and benthic parameters (e.g., sound speed deviation vs infauna biomass density) are not necessarily an indication of cause-and-effect relationship, such correlations or the lack of correlation are useful in establishing what parameter combinations might be important and could lead to more directed future research. Some general data-model deviation metrics are introduced first, followed by results of the correlation analysis focused on the high-frequency CARL data.

To quantify the data-model deviation at each measurement frequency and depth, the following quantity $δk,j$ is defined:

where $Dk,j$ and $Mk,j$ are the measurement quantity and model prediction, respectively, at the kth frequency and jth depth. Note that $δk,j$ accounts for the sign of the deviation (i.e., measurements higher or lower than predicted), which can be readily observed in Figs. 7–9.

Another useful metric is the cumulative data-model deviation, which is the root mean square over all frequencies for a given depth:

where N_j is the number of data-model pairs over all frequencies at a given depth. Vertical profiles of the cumulative data-model deviation can be computed separately for each of the three measured acoustic parameters (c_p, α_p, and c_s) to give a succinct picture of the data-model deviation as a function of sediment depth. In the following analysis and discussion, the subscripts j and k will be suppressed for conciseness. The frequency-dependent data-model deviation will be denoted $δcp, δαp$⁠, and $δcs$⁠, and the cumulative absolute-value data-model deviation will be denoted $|Δcp|,|Δαp|$⁠, and $|Δcs|$ for sound speed, attenuation, and shear speed, respectively.

There was lower cumulative data-model deviation (Fig. 10) associated with the high-frequency CARL sound speeds compared to the in situ sound speeds for the 5–20 cm depth interval, but absolute values of the cumulative sound speed error tend to stay below approximately 2% for both data sets in this interval. The sound speed deviation increases in the 1–4 cm interval, particularly for Site 1, where the absolute value of the sound speed deviation reaches as high as 5% in the top centimeter for the high-frequency CARL data.

Cumulative attenuation deviation (Fig. 10) is similar at both sites for the in situ and high-frequency CARL measurements in the 5–20 cm depth interval with absolute values under 50%, but the CARL attenuation measurements in the 1–4 cm interval show a steep increase in cumulative data-model deviation with absolute values reaching the 100%–200% range. The cumulative attenuation deviation skews higher at Site 1 than Site 2 in the upper 3 cm of sediment.

Finally, because shear speed data were only collected using the in situ system, measurements were only available in the 5–20 cm interval. The cumulative shear speed deviation is between 7 and 27% in this interval and suggests that deviation decreases with depth at Site 1 and increases with depth at Site 2. However, because shear speed data were not available in the 1–4 cm interval, it is not known whether any significant increase in the shear speed data-model deviation would be seen in this interval, as in the case of the sound speed and attenuation deviation.

A few comments are made on the acoustic data used for the correlation analysis. The high-frequency CARL data set provides greater depth resolution (1-cm steps in the 1–4 cm interval, 2-cm steps in the 6–20 cm interval) than the in situ measurements (5-cm steps in the 5–20 cm interval), providing a more detailed description of the depth-dependence, particularly in the uppermost layer of sediment. Only the high-frequency CARL data extend into the 1–4 cm depth interval and having data in this interval is crucial for comparison with the worm tube mass and infauna biomass, since this interval has the highest concentration of those quantities. The high-frequency CARL measurements extend to the frequency band where a significant deviation of the measured attenuation from the VGS model was observed (f > 180 kHz). Therefore, this data set is most useful in investigating whether the data-model deviation is correlated with the presence of the infauna, worm tubes, or shell hash.

Only CARL data from the 15.2-cm cores were used for this analysis since these cores provided the measurements of infauna biomass, worm tubes, and shell hash. The high-frequency CARL sound speed and attenuation data were sorted into bins that matched the 15.2-cm-diameter core sections so the two data sets could be compared. These bins were in the following intervals: 0–2.5 cm, 2.5–7.5 cm, 7.5–12.5 cm, 12.5–17.5 cm, and 17.5–22.5 cm.

The acoustic variables considered in this analysis were the depth-dependent cumulative and frequency-dependent sound speed and attenuation data-model deviation (⁠$|Δcp|, |Δαp|, δcp, δαp$⁠), and the biologically related variables considered were depth-dependent densities of infauna biomass, hard worm tube mass, and shell hash mass (⁠$ninf, ntube, nsh$⁠). For each acoustic-biological variable pair, the data from both sites were combined to increase the number of data points N, and the Pearson correlation coefficient r was computed. Additionally, p-values were computed for each variable pair to assess whether the reported correlation was statistically significant, using the p < 0.05 criterion.

There was statistically significant positive correlation between infauna biomass density and the cumulative data-model deviation for sound speed (r = 0.32, N = 76) and attenuation (r = 0.42, N = 78) (Fig. 11). The highest correlation was found between hard worm tube mass density and cumulative attenuation error (r = 0.62, N = 61), although this was driven primarily by two samples from Site 1 with very high worm tube density (Fig. 5). No significant correlation was found between hard worm tube mass density and cumulative sound speed deviation or between shell hash mass density and either sound speed or attenuation deviation.

Correlation coefficients were also computed as a function of frequency for the frequency-dependent data-model deviation $δcp$ and $δαp$⁠. Each correlation coefficient value was obtained from data pairs spanning all available sediment depths, similar to the scatter plots shown in Fig. 11, but for each measurement frequency (Fig. 12). Correlation of sound speed deviation with infauna biomass density had values between r = 0.3 and r = 0.4 for the 140–300 kHz band, whereas there was no significant correlation at any frequency between sound speed deviation and hard worm tube mass density or shell hash mass density. For attenuation data-model deviation, there was significant correlation with both infauna biomass and worm tube mass densities in the entire 100–300 kHz band (r > 0.5).

Although the two sites in this study were geographically very close and both contained primarily silt and sand, there were several key differences between the sites. Site 1 had a gradient from sandy to siltier sediment with increasing depth, whereas the sediment at Site 2 was more vertically mixed. Site 1 had nearly the same total infauna biomass as Site 2; however, the infaunal community at Site 1 exhibited a higher fraction of tube-building worms (group SHS) and approximately eight times higher hard worm tube density near the surface than Site 2, which is notable because hard worm tubes have been shown to increase attenuation at acoustic frequencies greater than 100 kHz.¹⁹ Also, noted, Site 1 had about twice as much shell hash as Site 2. Consistent with our hypotheses, Site 1, with greater biological activity, also had higher variability in acoustic properties and deviation from the acoustic model.

A deficiency of the data set is that there were only two sites, and they vary in more than two parameters. Therefore, it is difficult to isolate the effects of the worm tubes, the infauna themselves, specific functional groups, shell hash, or gradients in grain size and porosity between the two sites and which of these features result in higher variability and model deviation. We note that the grain size and porosity profiles are direct inputs to the VGS model; therefore, the depth-dependent heterogeneity in sediment type is captured to some extent in the depth-dependent profiles of acoustic parameters predicted by the model. Thus, by looking at where the data-model deviation is highest, we focus on biogenic effects not accounted for by the model.

Comparison of the measured sound speed and attenuation with the VGS model predictions indicated better overall data-model agreement at sediment depths greater than 5 cm below the sediment-water interface than in the near-surface layer. Increased variability in the sound speed and attenuation measurements also occurred in the upper 1–4 cm of the sediment at both sites. At both sites, most of the infaunal animals (∼70%) were found in the upper 7.5 cm of the sediment, consistent with global patterns.⁷ Of parameters not included in the model, neither shell hash mass density nor mean shell fragment size had significant vertical gradients, yet there were strong vertical gradients in the deviation and variability, suggesting that infauna or worm tubes had more influence than shell hash.

First, we examine volume scattering from shell hash as a potential cause of the data-model disagreement and the variability in the measurements. The shell hash mass density was relatively constant in the top 20 cm (Fig. 5), compared to stronger vertical infauna and worm tube gradients that were observed. The high-frequency attenuation data-model deviation occurs over greater portion of the measurement depth-interval instead of being limited to the top 1–4 cm, as the frequency increases (Fig. 8). This disagreement is most pronounced for the CARL measurements at the highest two or three frequencies (in the 180–300 kHz band). Multiple scattering effects are expected to become important for $kd≳0.5$⁠, where k is the wavenumber and d is the size of the scatterer, and can lead to negative sound speed dispersion and enhanced attenuation.⁵⁵ The shell fragment size in the measurement depth interval was $dsh=3.5 ± 1.6$ mm at Site 1 and $dsh=2.3 ± 1.3$ mm at Site 2 (mean ± standard deviation). Based on the average shell hash size and sound speed at each site, the corresponding frequencies above which scattering from shell hash might be expected to become important are approximately 350 kHz and 530 kHz at the two sites, respectively, which are higher than the upper end of the high-frequency CARL measurements. Although very little evidence of negative sound speed dispersion was observed (Fig. 6), it is clear that the measured attenuation dependence has a steeper slope than the approximately linear dependence predicted by the VGS model; therefore, it is possible that either (a) scattering from shell hash contributes to the higher than predicted attenuation or (b) the high-frequency attenuation dependence predicted by the VGS model is not correct for these sediments. The lack of correlation between the shell hash and the data-model deviation is likely influenced by the fact that there was less depth-dependence in the shell hash vertical distribution, whereas the acoustic variability and data-model deviation were greatest near the seabed surface.

Infauna biomass density, sound speed, and attenuation variability, and the data-model deviation were all highest in the sediment surface layer at both sites. Furthermore, Site 1 displayed greater deviation from model predictions in both sound speed and attenuation and greater variability in the measurements in the 1–4 cm depth interval than Site 2. The greater abundance of hard worm tubes near the surface at Site 1 offers a possible explanation for this difference between the two sites. Previous laboratory studies with simulated infauna in homogenized sediment indicated that the presence of a single worm tube comprised of shell material can impact sound speed and attenuation at high frequencies (f > 100 kHz).¹⁹ As much as 30 dB/m greater attenuation was observed with the structured tube in place in the laboratory experiment. In this previous study, the effect on sound speed was negligible; however, multiple worm tubes and burrows within the acoustic propagation path are expected to lead to a larger effect. Although such a controlled comparison (tubes/no tubes) could not be made in the present field experiment, the increased data-model attenuation deviation observed is likely, at least in part, due to the presence of multiple tubes and burrows. Further, the variability in the measurements in the top 1–4 cm interval is expected to result from patchiness in the distribution of worm tubes and burrows. The correlation analysis reported in Sec. VI suggests that the sound speed and attenuation data-model deviation was correlated to the total overall infauna biomass density and that there was an even stronger correlation between the attenuation data-model deviation and hard worm tube mass density, which provides a plausible explanation for the difference in variability and data-model deviation between the sites.

As part of the data-model deviation discussion, several possible sources of error in the measurements should be acknowledged, along with steps that were taken to mitigate that error. Specifically, we discuss the following here: (a) measurement of porosity from the cores, (b) the effect of coupling between the in situ probes and the sediment on the amplitude measurements, and (c) potential timing errors in the time delay measurements. Further detail on the measurement procedures and data analysis for all of the measured quantities are given in Ref. 44.

Regarding the porosity and bulk density measurements, evaluating water content from core samples is notoriously difficult because manipulation of the sediment invariably involves loss of intergranular water; however, efforts were taken to minimize water loss during collection and laboratory sampling of the cores.⁴⁴ Furthermore, the expressions used to derive porosity and bulk density from the water content measurements did not rely on precise knowledge of the volume of the sediment samples, but rather on the wet and dry weights of the samples, which were measured, and the densities of the pore water and mineral grains, which are well known,^2,44 thereby removing uncertainty contributed from the sample volume.

Another potential source of error in evaluating the data-model deviation could arise from inefficient coupling between the in situ probes and sediment, particularly in the top several centimeters of sediment where the sediment matrix is easily disturbed. Below 5 cm (the depth range of the in situ measurements), it was assumed that the loading of the sediment was the same on both in situ receivers, but it is possible that heterogeneity in the sediment physical properties or the distributions of shell hash, worm tubes, or infauna could cause the two receivers to sense different loading. However, we note that the sound speed and attenuation data from the upper few centimeters, where the data-model deviation was greatest, came from the high-frequency CARL measurements. For these measurements, constant coupling was maintained between the transducers and the core liner so this would not have been a factor in the data-model deviation observed in the top several centimeters (0–5 cm).

Time delay measurements depend on accurate determination of the cross correlation peak, which can be complicated when there is high attenuation or distortion of the signal by scattering from inhomogeneities. Time delays were computed using the cross correlation between the two receivers for the in situ probe measurements and between the source and receiver in the case of the high-frequency CARL measurements. To mitigate this problem, the data we manually checked and adjusted if obvious misalignment from selecting the wrong peaks was observed. In a very small number of ambiguous or unreliable cases, those data were omitted.

Very few previous sediment acoustics data sets include measurements of sound speed, attenuation, and shear wave speed in as large of a frequency band as presented here. These acoustic data coupled with the sediment physical parameter profiles and distributions of infauna, worm tubes, and shell hash provide a unique data set for testing sediment acoustics models such as VGS or others and for evaluating potential reasons that the models fail. The increased data-model deviation for at high frequencies is particularly interesting, as these are the frequencies at which we would expect inhomogeneities like infauna, tubes, and shell hash to have the largest effect.

The low-frequency CARL resonance mode sound speed measurements are consistent with the VGS model transition from high-frequency behavior to the Mallock-Wood sound speed low-frequency limit, which is $c0/cf≈0.985$ for both sites. The transition occurs at higher frequency (⁠$ft,p≈26.5$ kHz) than previously estimated for sand (⁠$ft,p≈1.3$ kHz)^37,40 or mud (⁠$ft,p<100$ Hz),^40,53,54 as shown in Fig. 13. Furthermore, the low-frequency shear speed measurements, with an average value of approximately 25 m s⁻¹, are consistent with a shear threshold frequency of $ft,s≈720$ Hz, suggesting a softer sediment than would be predicted by lower values of $ft,s$ that have been previously suggested for muddy sediment.^40,53,54 In the VGS model, $ft,p$ and $ft,s$ are threshold frequencies above which the effects of pore fluid viscosity become negligible for compressional and shear waves, respectively. This can also be interpreted as a threshold wavelength λ_t, related to the threshold frequencies via the wave speeds.³⁸ For the sediments here, $λt≈0.4$ cm, which represents a much smaller spatial scale compared to threshold wavelengths based on previous estimates for viscoelastic time constants for either sand or mud.

There are two plausible physical explanations for the best-fit values of the viscoelastic time constants being inconsistent with previous findings. One possible explanation is that organic matter and clay flocs likely fill much of the interstitial spaces between the larger silt and sand grains and the intergranular fluid behavior has changed. Quantification of the organic matter from the core samples, for instance in terms of total organic carbon or EPS, and comparison with reported values of organic carbon or EPS for other sediments would provide some insight. Although these data are not available for the present experiment, previous work has indicated that organic matter has an effect not only on the cohesive nature of sediment but potentially on the sediment acoustic properties as well.^42,56 An alternative possible explanation is that the fit is erroneously capturing some physical effect that is separate from the viscoelastic response of the sediment. Because the inhomogeneities (infauna, worm tubes, shell hash) are expected to have an impact on the sediment acoustic properties, it is reasonable to assume that they may impact those properties in ways that change the frequency-dependence of the sediment acoustic response and, in turn, impact the values of the viscoelastic time constants determined from the model-data fit. This does not necessarily indicate that the physics behind the time constants have changed, but instead it suggests that erroneous values may have been found from the fitting process that does not reflect the actual physics taking place between the grains.

Finally, the choice of sediment acoustics model is also a factor that should be considered regarding the data-model deviation. The VGS model does not include losses arising from poroelastic effects,³⁴ which would occur at higher frequencies if relative motion between the pore fluid and the sediment skeletal frame is important in this sediment. This additional loss mechanism leads to a different frequency dependence at high frequency (proportional to f²), and hence using a poroelastic model for mud, such as the one proposed by Chotiros,⁵⁷ could potentially improve the data-model agreement; however, how important poroelastic effects are in sediments in which the interstitial spaces are largely filled with clay or organic matter is an open question. Current implementations of the VGS and Biot models take only a single grain size as an input parameter, most often the mean grain size, which is most appropriate for well sorted sediments, but perhaps not as appropriate for poorly sorted sediments, as found in this work. The silt-suspension mud model^58,59 explicitly incorporates the full grain size distribution into a card-house matrix of clay platelets and flocs; however, the silt-suspension model likely does not apply to sediment containing as much sand as the one here, as large sand grains would not be expected to hold in suspension. None of the models mentioned here explicitly include effects of soft organic matter in the interstitial spaces or on the surfaces of larger silt or sand grains. A model that incorporates the full grain size distribution, effects of clay flocculation and organic matter in interstitial spaces, and larger scatterers (i.e., worm tubes or shell hash) would likely better account for the variability observed in the acoustic data, as long as the model inputs could be well defined; however, such a model does not currently exist to the best knowledge of the authors.

The measurements of sediment acoustic, physical, and biological parameters presented in this paper suggest that infauna and their activities (a) increase variability in sound speed and attenuation and (b) lead to deviation from the VGS model predictions. Measurements of sound speed and attenuation indicated a higher degree of acoustic variability near the sediment surface where most of the infauna were present than at greater depths within the sediment. The acoustic variability was greatest at the site with more hard worm tube mass. The sound speed and attenuation data were compared with the VGS model to assess where the model was in agreement with the measured wave speeds and attenuation and where it was deficient. Depth-dependent porosity and mean grain size were input directly into the VGS model to account for variation in the sediment acoustic parameters due to variation in sediment physical parameters. Overall, the VGS model compared well with the measured sound speeds across a wide range of frequencies; however, significant deviation between the measured high-frequency (100–300 kHz) sound speeds and the VGS model was observed in the depth interval closest to the sediment surface. Attenuation measurements were typically higher than the VGS predictions, again particularly close to the sediment surface where infauna and worm tubes were most abundant and also at the highest frequencies (> 180 kHz). We note that the VGS model is intended to provide a prediction of the minimum intrinsic attenuation in clean sediment without any other components, and should not be expected to predict the overall effective attenuation when infauna, worm tubes, shell hash, or other inhomogeneities are present in the medium,³³ and the measurements presented in this paper are consistent with that notion.

Another finding of this work was that the VGS compressional viscoelastic time constant determined from the data-model fit was inconsistent with previous estimates for either sand or mud, and hence the associated threshold frequency was much higher. There are two possible explanations for this inconsistency that are related to the infauna: (1) the properties of the inter-granular pore fluid have been altered by seabed biological activity or (2) the fit has erroneously captured an unknown physical effect that is separate from the viscoelastic response. Additionally, the shear speed measurements indicated that the sediment at both sites was softer than would be predicted by previous implementations of the VGS model that did not include sufficiently small shear viscoelastic time constant. The fact that the shear viscoelastic time constant is coupled to the compressional viscoelastic time constant though the threshold wavelength suggests that the aforementioned possibilities also apply to the shear wave. However, the true causes of these disparities cannot be determined from the present work, and their resolution is left for a topic of further research.

This work was supported by the Office of Naval Research Grant Nos. N00014-18-1–2227, N00014-21-1–2254, N00014-15-1–2602, N00014-17-1–2625, N00014-21-1–2245, and the ARL:UT Independent Research and Development Program. The authors would like to thank the participants in the field experiment: Will Ballentine, Grant Lockridge, Cy Clemo, Kara Gadeken, and Ryan Parker from Dauphin Island Sea Lab; and the crew of the R/V E. O. Wilson. The authors also thank Kimbell Bui and Jeremy King, formerly of ARL:UT, who assisted with analysis of physical properties from the cores; Dan Duncan of the University of Texas at Austin Institute for Geophysics for assistance with the laser diffraction particle size analyzer; and Dr. Allen Reed for lending us the core-logger transducers for this experiment. The contributions of the anonymous reviewers are also gratefully acknowledged.