# An explicit granular-mechanics approach to marine sediment acoustics

Abram H. Clark<EMAIL_ADDRESS>Physics Department, Naval Postgraduate
School, Monterey, CA 99343  Derek R. Olson<EMAIL_ADDRESS>Oceanography
Department, Naval Postgraduate School, Monterey, CA 99343  Andrew J. Swartz
Physics Department, Naval Postgraduate School, Monterey, CA 99343  W. Mason
Starnes  Physics Department, Naval Postgraduate School, Monterey, CA 99343

###### Abstract

Here we theoretically and computationally study the frequency dependence of
phase speed and attenuation for marine sediments from the perspective of
granular mechanics. We leverage recent theoretical insights from the granular
physics community as well as discrete-element method simulations, where the
granular material is treated as a packing of discrete objects that interact
via pairwise forces. These pairwise forces include both repulsive contact
forces as well as dissipative terms which may include losses from the fluid as
well as losses from inelasticity at grain-grain contacts. We show that the
structure of disordered granular packings leads to anomalous scaling laws for
frequency-dependent phase speed and attenuation that do not follow from a
continuum treatment. Our results demonstrate that granular packing structure,
which is not explicitly considered in existing models, may play a crucial role
in a complete theory of sediment acoustics. While this simple approach does
not explicitly treat sound propagation or inertial effects in the interstitial
fluid, it provides a starting point for future models that include these and
other more complex features.

##  I Introduction

The dispersion relation for acoustic waves is a fundamental observable for
materials. The functions $c(f)$ and $\alpha(f)$ connecting the phase speed $c$
and attenuation coefficient $\alpha$ to the frequency $f$ depend strongly on
the microscopic processes that control wave propagation. Marine sediments are
fundamentally fluid-saturated granular materials, and there are many complex
physical aspects to consider, even for an idealized scenario (e.g., neglecting
biological matter ballard2017acoustics ). For example, the nonlinear Hertzian
daraio2005strongly ; gomez2012shocks  and inelastic contact forces
force_schemes  between the grains, the structure and other properties of the
granular contact network majmudar2005contact ; Silbert_PRL_2005 , and the
details of the interactions with the interstitial fluid chotiros2017Acoustics
; Williams2009b  all may affect the acoustic response.

One of the important applications of dispersion models is their use in
acoustic remote sensing of the seafloor, which commonly uses measurements of
the reflection coefficient Holland2012 ; Dettmer2010 ; Quijano2012 ;
Holland2017 ; Olson2023  or waveguide dispersion Bonnel2013 ; tan_ambient_2022
, combined with an acoustic model, to estimate the sound speed and density as
a function of depth within the seafloor. If these measurements are made over a
wide frequency band, then dispersion will be important Holland2013 . Another
application of sediment dispersion models is the prediction of the bottom-
interacting acoustic field due to a source (i.e. propagationPorter1985  or
reverberationHolland1998 ; LePage2000 ; LePage2003 ; Tang2017 ; Olson2020 ;
Jackson2020 ; Olson2020a ), with the goal of performance modeling of remote
sensing of objects in the ocean Abraham2019 . For these applications, linking
the microscopic properties with the observables is a key step.

Real-world data for $c$ and $\alpha$ can be difficult to obtain for a wide
range of $f$ and there is still debate about the functional forms that best
describe these materials. One compilationzhou_low-frequency_2009  of data
shows distinct behavior at low $f$, with $c(f)$ constant and $\alpha(f)\propto
f^{2}$, and high $f$, with $c(f)$ increasing and $\alpha(f)\propto f^{1/2}$,
with some crossover region characterized by a frequency $f=f^{*}$. This
behavior is also typical of plane wave propagation in a viscous medium, and
will be termed “classical viscosity.” These data can be fit to the Biot-Stoll
model and its extensions biot_theory_2005-1 ; biot_theory_2005 ;
biot_generalized_2005 ; stoll_acoustic_1977 ; stoll_theoretical_1980 ;
stoll_reflection_1981 ; stoll_marine_1985 ; chotiros2017Acoustics , which
treat the sediment as a porous medium in which the frame can support
elasticity, and losses result from the relative motion between the fluid and
the frame due to the viscosity of seawater. Its adaptation to sediments
introduced losses within the poroelastic frame using a viscoelastic solid
modelstoll_marine_1985 . The majority of the Zhou, et al. zhou_low-
frequency_2009 , data that show this dependence are based on indirect
dispersion measurements, such as bottom reflection or waveguide propagation.
It is also assumed that the sediment is a homogeneous half space. At those low
frequencies, the acoustic waves interact with at least several meters of the
sediment, which often contains gradients and layering structure. As recent
work Godin_2021 ; godin_effects_2023  has shown, slight layer contrasts within
a weakly shear supporting sediment may have significant impact on modal
attenuation. Therefore, we focus on only direct measurements of sediment
dispersion, either in situ, laboratory, or using an idealized sediment
composed of glass beads.

Many of the direct measurements of dispersion in sandy sediments depart from
classical viscosity, most notably showing an approximate $\alpha(f)\propto f$
relationship. The field studies cited here are from the measurements compiled
by Hamilton Hamilton1980 , the large SAX99 williams2002comparison , and SAX04
Hefner2009  experiments, as well as field measurements by Turgut and
YamamotoTurgut1990  and Simpson, et al.Simpson2003 . Laboratory data using
sand are collected from Wingham Wingham1985 , and laboratory measurements
using glass spheres were made by Hefner and Williamshefner2006sound  (using
both water and oil as the saturating fluid). These measurements are collected
in Fig. 1, and show a deviation from the Biot-Stoll model, and agreement with
$\alpha\propto f$ for frequencies $1<f<400$ kHz. We note that several
experimentsLee2007 ; Sessarego2008  were excluded from our compilation of
experimental data in Fig. 1 due to the presence of negative velocity
dispersion, an indication that the wavelength was close to the mean grain size
of the medium (see Hare and HayHare2020  and references therein).

Figure 1: A compilation of direct measurements of the attenuation coefficient
of saturated sand or glass beads. Field measurements from
SAX99williams2002comparison , SAX04Hefner2009 , HamiltonHamilton1980 , Simpson
et al.Simpson2003 , and Turgut and YamamotoTurgut1990  are included here. The
laboratory measurements are from from WinghamWingham1985  and Hefner and
Williamshefner2006sound . The solid black line shows $\alpha\propto f$.

Trying to understand how such an attenuation form (specifically the
$\alpha\propto f$ trend observed in SAX99Hamilton1980 ) might arise was in
part the motivation for Buckingham’s grain-shearing (GS) Buckingham1997 ;
buckingham_wave_2000 , theory of sediment acoustics, and the later addition of
viscosity buckingham_pore-fluid_2007 . In his justification of this theory,
Buckinghambuckingham_wave_2000  noted that “grains in an unconsolidated
granular material do not form a macroscopic elastic skeletal frame,” but
instead that “forces within the unconsolidated medium arise from grain-to-
grain interactions” which “give rise to force chains.” The term “force chains”
refers to the fact that granular contact networks are spatially disordered,
and the forces between grains are not homogeneous majmudar2005contact .

Buckingham engaged with some of the pioneering work in the field of granular
mechanics jaeger1996granular  to justify the need for a granular-focused
approach. However, the key mechanisms in the GS theory were the frictional and
time-dependent characteristics of the grain-grain contacts and not the spatial
structure of the “force chains” explicitly. The time-dependent damping law
used for grain-grain interactions was the key to producing $\alpha\propto f$.

However, Buckingham’s original works were published during the infancy of a
revolution forterre2008flows ; van2009jamming ; liu2010jamming ;
behringer2018physics  in our understanding of the mechanical properties of
dense granular media, facilitated by improvements in experimental
amon2017preface  and simulation radjai2017modeling ; kamrin2024advances
methods and technologies. One key result is that disordered packings often
have emergent mechanical properties that are not a priori obvious or
consistent with a simple continuum picture.

For example, disordered packings of elastic disks (2D) or spheres (3D) exhibit
an excess of low-frequency vibrational modes Silbert_PRL_2005 ; Wyart_PRE_2005
. In uniform elastic solids, low-frequency vibrations are assumed to be
acoustic plane waves obeying Debye scaling, which implies a density of modes
$\mathcal{D}(\omega)\propto\omega^{2}$ in three dimensions, where $\omega=2\pi
f$ is the angular frequency. In contrast, the density of vibrational modes in
disordered packings is empirically found to be very sensitive to the
dimensionless pressure $\hat{P}$, defined as the ratio of system pressure to
the bulk modulus of an individual grain. At large $\hat{P}$ (unrealistic for
geological granular materials including marine sediments), the density of
modes is consistent with the Debye picture, but these particles are “soft” and
possess many more contacts than necessary. As $\hat{P}$ is decreased, the
system develops an excess of low-frequency vibrational modes compared with
Debye scaling Silbert_PRL_2005 .

How packing structure affects dispersion and attenuation in disordered
packings is an open question, the answer to which would have obvious
implications for marine sediment acoustics as well as other systems. We note a
recent paper saitoh2023sound , which studied standing modes in disordered 2D
packings. These authors demonstrated that, even at large $\hat{P}$, the
disordered packing structure affects the dispersion relation and attenuation
rate (in time) of attenuation of plane wave modes. However, many open
questions remain, including the behavior at small $\hat{P}$, which is more
appropriate to geological granular materials like marine sediments.

To this end, we here consider discrete-element method (DEM) simulations of
granular packings and show that the granular packing structure plays a crucial
role in determining the dispersion and attenuation curves. DEM simulations
allow us to postulate grain-grain interaction rules and then solve the
equations of motion for each grain individually. Thus, our approach explicitly
considers the granular packing structure without assuming a continuum wave
equation. This approach has limitations, which we discuss below in Sec. II;
however it is able to isolate the effect of the granular packing structure on
sound speed and attenuation, which has not been explicitly considered
previously.

We use linearized forces (linear springs and dashpots) between every grain,
which is mathematically similar to viscous drag. The repulsive springs
represent repulsive intergrain interactions, and the dashpots represent the
combined effect of lossy grain-grain contacts and the effect of drag from the
interstitial fluid (see Sec. II for further discussion). Using these
linearized forces we theoretically and computationally obtain the standard
result of $\alpha\propto f^{2}$ at low $f$ and $\alpha\propto f^{1/2}$ at high
$f$ in the continuum limit for a 1D periodic chain of spheres. However, in
simulations of packings in two (2D) and three (3D) dimensions, our results
show emergent scaling, where $\alpha\propto f^{\beta}$ at low $f$ with
$\beta<2$; the 3D simulations demonstrate that $\beta\approx 7/6$, which is
very near 1. We find that 2D and 3D packings revert back to the continuum
behavior at large $\hat{P}$, similar to the disappearance of excess low-
frequency modes at large $\hat{P}$ discussed above.

Our results suggest that $\alpha\propto f^{\beta}$ with $\beta\approx 1$
scaling may be a direct result of the packing structure of granular materials,
even with simple linear forces. The finding of $\alpha\propto f^{1}$ is
consistent with many of the direct measurements of dispersion in sandy
sedimentswilliams2002comparison ; hefner2006sound ; Wingham1985 ; Hamilton1980
. Our results may be sensitive to other granular effects, including Hertzian
grain interactions daraio2005strongly ; gomez2012shocks , which are more
realistic and involve a contact stiffness that increases with relative
displacement. Other nonlinearities include contact breaking owens_sound_2011 ;
Schreck_PRL_2011 , frictional grain interactions silbert2010jamming . Non-
spherical grain geometry nguyen2014effect ; nguyen2015effects  may also play a
role. These effects may be considered in future studies. The problem of
integrating the inertia and compressibilityWilliams2009b  of the interstitial
fluid is key to explaining the behavior of in situ dispersion measurements and
is another area for future work. However, the key result is that the granular
packing structure plays a crucial role in determining the scaling laws for
$c(f)$ and $\alpha(f)$ and therefore should be more explicitly included in
theories of sediment acoustics.

The remainder of this paper is organized as follows. In Sec. II, we consider
the continuum limit of a 1D granular model. In Sec. III, we describe the
numerical simulations we use to compare to this continuum theory. Results of
the numerical simulations are shown in Sec. IV, where we demonstrate how
anomalous scaling laws emerge in 2D and 3D due to the granular packing
structure. In Sec. V we provide discussion and concluding remarks.

##  II 1D linear continuum theory

We first consider the continuum limit of a 1D granular model, based on the
dynamics of a single “force chain” of grains, as shown in Fig. 2. We use this
theory as a baseline to compare numerical simulations in 1D, 2D, and 3D. Each
grain has identical diameter $d$ and mass $m$. We assume grains $i$ and $i+1$
interact via dissipative, repulsive contact forces. These forces are
linearized by assuming that the repulsive forces use linear springs and the
dissipative forces use simple dashpots. Physical grains have nonlinear
repulsive forces between them (e.g. the Hertz contact law coste1999validity ),
but a linear approximation is reasonable for small displacements such as
acoustic strains.

One limitation of this approach is that the grains must be considered as
perfectly rigid, and compression of grains is modeled by the finite overlap
between them. This approximation is only valid in the limit where time scales
associated with dynamics inside the grains (in this case, internal elastic
waves) are very small compared with other dynamical time scales in the system
(e.g., shear rates or, in this case, propagation of waves in the bulk
material). Marine sediments are on the edge of validity of such
approximations, since the sound speed in, e.g., silica grains is only about
four times that of water. Future work should more closely examine if and how
treating the grains as objects with internal elasticity (e.g., an FEM
approach) affects the results.

The linear dashpot term assumes that the dissipation is directly proportional
to the relative velocity between the grains at any given time. Such a
functional form is meant to represent some combination of effects from the
interstitial fluid as well as inelastic effects that are inherent to grain-
grain contacts.

Figure 2: A one dimensional chain of identical spherical grains.

With these assumptions, the displacement $u_{i}$ of grain $i$ from its
unperturbed position $x$ on the chain is then governed by the following
differential equation, which includes a repulsive spring force and a dashpot
for each neighbor $i-1$ at position $x-d$ and $i+1$ at position $x+d$:

$\displaystyle m\ddot{u}_{i}=\kappa$
$\displaystyle(u_{i-1}-u_{i})-\kappa(u_{i}-u_{i+1})$
$\displaystyle+\gamma(\dot{u}_{i-1}-\dot{u}_{i})-\gamma(\dot{u}_{i}-\dot{u}_{i+1})\,,$
(1)

Here, $\kappa$ is a spring constant, $\gamma$ is the dashpot coefficient, and
dots denote time derivatives. Taking the long-wavelength limit (i.e., assuming
$d$ is small enough to transform the finite differences into derivatives)
yields

$mu_{tt}=\kappa d^{2}u_{xx}+\gamma d^{2}u_{txx},$ (2)

where subscripts $t$ and $x$ now denote derivatives with respect to time and
space, respectively. A decaying plane wave is guessed, given by

$u(x,t)=Ae^{-\alpha x}e^{i(kx-\omega t)},$ (3)

where $\omega$ is the angular frequency, $k$ is the acoustic wavenumber, and
$\alpha$ is the attenuation parameter with units of radians per meter. We can
then solve the system to arrive at the two following equations for non-
dimensional wavenumber and attenuation:

$\displaystyle\hat{\alpha}$ $\displaystyle=$
$\displaystyle\hat{\omega}^{2}\frac{\hat{\gamma}}{\sqrt{2}}\left\\{\left[1+(\hat{\omega}\hat{\gamma})^{2}\right]\left[1+\sqrt{1+(\hat{\omega}\hat{\gamma})^{2}}\right]\right\\}^{-1/2}$
(4) $\displaystyle\hat{k}$ $\displaystyle=$
$\displaystyle\hat{\omega}\frac{1}{\sqrt{2}}\frac{\left\\{\left[1+(\hat{\omega}\hat{\gamma})^{2}\right]\left[1+\sqrt{1+(\hat{\omega}\hat{\gamma})^{2}}\right]\right\\}^{1/2}}{1+(\hat{\omega}\hat{\gamma})^{2}}\,.$
(5)

The dimensionless quantities in Eqs. (4) and (5) are defined by the following
equations:

$\displaystyle\hat{k}$ $\displaystyle=kd$ (6) $\displaystyle\hat{\alpha}$
$\displaystyle=\alpha d$ (7) $\displaystyle\hat{\omega}$
$\displaystyle=\omega\sqrt{m/\kappa}$ (8) $\displaystyle\hat{\gamma}$
$\displaystyle=\gamma/\sqrt{\kappa m}\,.$ (9)

The wavenumber and attenuation coefficient each have dimensions of inverse
meters, and are easily non-dimensionalized by the grain diameter. Angular
frequency is non-dimensionalized by the resonance angular frequency for an
individual grain contact, and the dashpot term is non-dimensionalized using
the ratio of the energy lost and energy stored per cycle, which is related to
the Q factor of the contact Garrett2020a .

We note that Eq. (4) gives $\alpha\propto\omega^{2}$ for low frequencies and
$\alpha\propto\omega^{1/2}$ at high frequencies. For wave speed, the result is
constant at low frequencies and then increases according to
$c\propto\omega^{1/2}$ at high frequencies. Equation (2) has a similar
mathematical form to classic viscosity shitikova2022models , so these results
are expected.

##  III Methods

We now discuss the methods used to perform DEM numerical simulations in 1D,
2D, and 3D of granular packings subjected to harmonic driving at a boundary.

###  III.1 Numerical Integration Scheme

For each grain $i$, we use a modified velocity-Verlet binder2004molecular
scheme,

$\displaystyle\mathbf{r}(t+dt)$
$\displaystyle=\mathbf{r}(t)+\mathbf{v}(t)dt+\frac{1}{2}\mathbf{a}(t)dt^{2}$
(10) $\displaystyle\mathbf{v}(t+dt)$
$\displaystyle=\mathbf{v}(t)+\frac{1}{2}[\mathbf{a}(t)+\mathbf{a}(t-dt)]dt$
(11)

to integrate Newton’s equations of motion. At each time step, the acceleration
for each grain $i$ is calculated using

$m_{i}\mathbf{a}_{i}=\sum_{j}\mathbf{F}_{c}^{ij}+\sum_{j}\mathbf{F}_{d}^{ij}$
(12)

Here, $m_{i}$ is the mass of grain $i$; $\mathbf{r}_{i}$, $\mathbf{v}_{i}$,
and $\mathbf{a}_{i}$ are the vector position, velocity, and acceleration of
grain $i$; $\mathbf{F}_{c}^{ij}$ are repulsive contact forces between grains
$i$ and $j$ that depend on their positions; $\mathbf{F}_{d}^{ij}$ are
dissipative forces that depend on the velocities of grains $i$ and $j$; and
the sum over $j$ denotes a sum over all grains $j$ that are in contact with
grain $i$.

###  III.2 Intergrain forces

We use linearized repulsive forces that involve overlaps
$\xi_{ij}=(R_{i}+R_{j}-|\mathbf{r}_{i}-\mathbf{r}_{j}|)\Theta(R_{i}+R_{j}-|\mathbf{r}_{i}-\mathbf{r}_{j}|)$.
Here, $R_{i}$ and $R_{j}$ are the radii of particles $i$ and $j$ and $\Theta$
is the Heaviside function, which sets $\xi$ to zero when particles are not in
contact. The repulsive force is modeled in terms of this overlap,

$F_{c}^{ij}=-\kappa\xi_{ij}\hat{\mathbf{n}}_{ij}\,,$ (13)

where the normal vector $\hat{\mathbf{n}}_{ij}$ has magnitude of 1 and points
along the line from the center of particle $i$ to the center of particle $j$.
We note that the amplitudes involved in acoustic waves are sufficiently small
that individual grain-grain contacts are well approximated by linear repulsive
forces.

We also use linearized dissipative forces using a damping coefficient
$\gamma$. For the bulk of the results we show, the dissipative forces can be
simply written as

$F_{d}^{ij}=-\gamma(\mathbf{v}_{i}-\mathbf{v}_{j}).$ (14)

We note that this form damps all relative motion between grains, both normal
and tangential. Such a form is intended to be the simplest form that mimics a
combination of lossy grain-grain contacts combined with fluid-mediated damping
due to relative motion in saturated granular packings.

The functional form in Eq. (14) assumes that normal and tangential motion are
damped with the same magnitude (note that this only applies in 2D and 3D,
since there is no tangential motion in 1D). However, even if a linear
approximation is reasonable, there are likely differences in the magnitude of
damping for normal and tangential relative motion between grains. In this
case, we can write in 2D, using
$\delta\mathbf{v}_{ij}\equiv(\mathbf{v}_{i}-\mathbf{v}_{j})$,

$\displaystyle F_{d}^{ij}=$
$\displaystyle-\gamma_{n}\hat{\mathbf{n}}_{ij}(\delta\mathbf{v}_{ij}\cdot\hat{\mathbf{n}}_{ij})-\gamma_{t}\hat{\mathbf{s}}_{ij}(\delta\mathbf{v}_{ij}\cdot\hat{\mathbf{s}}_{ij}).$
(15)

Here, $\mathbf{s}_{ij}$ is a tangential unit vector such that
$\mathbf{s}_{ij}\cdot\mathbf{n}_{ij}=0$. If $\gamma_{n}=\gamma_{t}$, this form
reduces to the form given in Eq. (14). While we primarily use Eq. (14), we
also show simulations in 2D using Eq. (15) with $\gamma_{t}=0$ to support our
claims in Sec. IV.2. In 3D, two tangential unit vectors would be required.

We also note that use of $\gamma_{n}$ and $\gamma_{t}$ is common in DEM
simulations of dry granular media (i.e., with no interstitial fluid). This is
because grain-grain collisions involve restitution losses, i.e., the
conversion of some of kinetic energy of the grains into heat. The magnitude of
these losses is often parameterized by a restitution coefficient $e_{n}$,
which can be connected to the magnitude of $\gamma_{n}$ for perfectly normal
collisions. In particular, for linear interactions force_schemes ,

$e_{n}=\exp\left(-\frac{\gamma_{n}}{2m_{\rm eff}}t_{n}\right),$ (16)

with

$t_{n}=\pi\left[\frac{\kappa}{m_{\rm eff}}-\left(\frac{\gamma_{n}}{2m_{\rm
eff}}\right)^{2}\right]^{-1/2}$ (17)

where $m_{\rm eff}$ is the reduced mass. These normal losses have not been
considered by the acoustics community to date, although the work of
BuckinghamBuckingham1997  focused on modeling of a time-dependent
$\gamma_{t}$-like term. Some evidence for grain-based losses is also provided
by Hefner and Williamshefner2006sound  who found that increasing the kinematic
viscosity of the pore fluid by a factor of 100 only increased the attenuation
coefficient at the same frequency by a factor of ten. This situation is
possible if the effective $\gamma$ from grain contacts is ten times that of
viscous drag from water, and one tenth that of viscous drag from silicone oil.
Given estimates of $e_{n}$ from glass beads of approximately 0.98 for small
velocitiesMontaine2011  the resulting $\gamma_{n}$ results in the right order
of magnitude to explain the results of Hefner and Wiliamshefner2006sound .

Thus, a linearized dashpot, which applies a force proportional to the
magnitude of relative velocity between particles, is a reasonable model to use
both for viscous-like forces as well as inherently lossy grain-grain contacts.
Our $\gamma$ parameters are meant to capture some combination of these
effects. Using $\gamma_{n}$ to control $e_{n}$ models losses during
collisions, not necessarily due to relative motion during harmonic oscillation
of a packing, as we do here. In the latter case, the contacts are enduring.
Future work is needed to better understand the relative contributions of
losses due to viscous effects and lossy contacts.

###  III.3 Generating Packings at Different Pressures

In 1D, initial configurations are trivial to generate, as they constitute a
periodic chain. However, in 2D and 3D, grains must be prepared in a disordered
packing with all grains in force balance. Here we describe how disordered
packings are generated using DEM simulations, and we discuss both 2D and 3D
situations at the same time, calling out differences as appropriate. Unless
otherwise stated, the discussion below applies to both 2D and 3D packing
generation protocols.

Disorder is enforced using two diameters of grains with ratio 1.4, with equal
number of particles of each size ohern2003jamming . Using two grain sizes
ensures that 2D packings do not end up with a periodic, crystaline structure.
We consider $N$ grains, where $N$ varies between 5,000 and 20,000, half with
diameter $d$ and half with a larger diameter $1.4d$. Despite having different
diameters, all grain masses are set to the same value $m$. Grains are
initially placed randomly throughout a large domain with fixed walls in the
$x$-direction at $x=0$ and $x=L_{x}$ and periodic boundaries in the $y$ and
$z$ directions with length $L_{y}$ and $L_{z}$ (in 3D only). The periodic
boundaries simulate an infinite medium in the transverse directions, while a
nonperiodic (wall-based) boundary is required to observe the propagating and
decaying wave as it moves down the channel.

During the computationally intensive initial compression, we also include a
dissipative force that is linearly dependent on the absolute velocity of the
grain, as if the grains were in a highly viscous background fluid.
Additionally, we set the velocity and acceleration to zero for any particle
with no contacts, which significantly speeds up the compression phase. These
additional damping mechanisms are turned off after the initial packings are
generated.

Compression of the system is done in the following way. Grains are initially
placed randomly, without touching, with very large $L_{y}\approx L_{x}$. The
value of $L_{x}$ is held fixed throughout compression, typically between
$200<L_{x}/d<1000$. Additionally, in 3D, $L_{z}$ is held fixed throughout the
simulation, typically between $5\leq L_{z}/d\leq 10$. At each compression
step, we compress the $y$-direction using the substitution $L_{y}\rightarrow
L_{y}(1-\Delta)$ as well as substituting the $y$-position of each particle
using $y_{i}\rightarrow y_{i}(1-\Delta)$, where $\Delta$ is chosen to achieve
the desired $\hat{P}$ (typically, $\Delta\approx\hat{P}$. After such a
compression step, we allow the particles to rearrange and the energy to
dissipate by evolving the simulation in time. When the total kinetic energy
becomes sufficiently small, we compress again and again allow the system to
relax. During compression, we track the mean kinetic energy per particle
$\bar{K}=\frac{1}{N}\sum_{i}\frac{1}{2}mv_{i}^{2}$ and potential energy per
particle $\bar{U}=\frac{1}{2N}\sum_{i}\sum_{j}\frac{1}{2}\kappa\xi_{ij}^{2}$,
where the extra factor of 2 in the denominator of $\bar{U}$ accounts for
double counting of contacts.

The system pressure is defined $P=\frac{1}{d^{D-1}}\sqrt{2\bar{U}\kappa}$,
where $D$ is the number of spatial dimensions. We use the dimensionless
pressure $\hat{P}$ to characterize the packings, with $\hat{P}=P/\kappa$ in 2D
and $\hat{P}=Pd/\kappa$ in 3D, i.e.,

$\hat{P}=\sqrt{\frac{2\bar{U}}{\kappa d^{2}}}.$ (18)

Physically, $\hat{P}$ can also be thought of as an approximation of the
characteristic size of the overlap between two particles compared to a
particle size, $\bar{\xi}/d$. Realistic values of $\hat{P}$ for marine
sediments may be extremely small, since confining pressures are generated by
the weight of grains above them, e.g., $\rho gh\approx 25$ kPa for $\rho=2500$
kg/m3, $g=9.81$ m/s2, and $h=1$ m. The stiffness of the grains is tens of GPa,
so $\hat{P}\approx 10^{-6}$ may be a physically reasonable value. Although
packings at such low $\hat{P}$ are difficult to generate computationally, it
is often true that there is a low-pressure limit at a value of $\hat{P}$ that
is much higher and therefore much easier to generate. In this low-pressure
limit, the dynamics are insensitive to the value of $\hat{P}$ so long as it is
sufficiently low.

Thus, our approach here is to prepare packings at a range of values of
$\hat{P}$ spanning $10^{-1}$ down to $10^{-4}$ until a low-$\hat{P}$ limit is
reached. Despite the fact that the large-$\hat{P}$ packings are unrealistic,
our approach has two major benefits. First, it allows us to show that the low-
pressure limit occurs for our simulations at $\hat{P}\approx 0.005$, and
further decreasing $\hat{P}$ has no effect on the results. Second, our
approach allows us to show that the acoustic properties of granular packings
agree with the 1D continuum theory at large $\hat{P}$, whereas anomalous
scaling appears as the low-pressure limit is approached. Packings in 2D at
high and low $\hat{P}$ are shown in Fig. 3. The high $\hat{P}$ packing is far
from the unjamming transition Silbert_PRL_2005 ; van2009jamming , where the
packing would break apart and have no intergrain forces, which can be seen
from its large grain-grain overlaps and relatively large number of contacts
per grain. It is thus expected to behave more like a continuum elastic solid
than the low $\hat{P}$ packing.

(a)

(b)

Figure 3: Two packings, one at $\hat{P}\approx 0.1$ (a) and another at
$\hat{P}\approx 0.001$ (b).

##  IV Results

###  IV.1 1D Simulations

We first perform simulations in 1D. These simulations both confirm our
theoretical analysis above as well as demonstrate that the effects observed in
higher dimensions are not an artifact of DEM simulations but a result of
disordered granular mechanics in higher spatial dimension.

The initial positions of grains are prepared, in 1D, by assigning the grains’
positions at equally spaced positions similar to Fig. 2. In 1D,
$\mathbf{r}_{i}=x_{i}\hat{\mathbf{x}}$, where $\hat{\mathbf{x}}$ is a unit
vector along the chain. We study 10,000 equal sized grains placed in equal
spacing $d=1$ along a line. Each grain is subject to forces from two
neighbors, according to Eq. (1). The first grain in the chain is driven in one
of two ways. The first method, constant driving, begins driving at $t=0$
according to $x_{1}(t)=x_{1}(0)+A\sin(\omega t)$. We choose $A/d=10^{-4}$ and
have verified that our results are not sensitive to this choice. Second, we
drive the grain according to a Gaussian-enveloped pulse with a primary angular
frequency of $\omega$ and a bandwidth of $\omega/5$. These two types of
driving yield similar results.

We measure $\alpha$ for a given $\omega$ by measuring the peak height of
oscillation of each grain as a function of its distance down the chain and
fitting to a decaying exponential. For the constant-driving case, we measure
the wavenumber $k$ by observing the phase $\phi$ of each grain and fitting to
$\phi=kx$. This allows us to infer $c=\omega/k$. For the Gaussian-enveloped
pulse, we measure $c$ using the time corresponding to the largest displacement
of each grain. The results from these two methods are indistinguishable.

Results of these measurements are shown in Fig. 4. We systematically vary
$\kappa$, $\omega$, and $\gamma$ in order to access a wide range of values of
$\hat{\omega}\hat{\gamma}$. The measured data from simulations agree well with
the continuum theory from Sec. II over a large range of $\hat{\gamma}$.

Figure 4: Results for 1D simulations. Dimensionless wave speed, $\hat{c}$, and
ratio of attenuation to frequency, $\hat{\alpha}/\hat{\omega}$, both plotted
as a function of dimensionless frequency times dimensionless damping,
$\hat{\omega}\hat{\gamma}$. Solid black lines show the theoretical result from
Eqs. (4) and (5).

The fact that the 1D simulations do not significantly deviate from the 1D
theory is not surprising, since the theory was developed to model this exact
situation. These results are presented as an example to provide confidence
that the observed emergent scaling behavior in 2D and 3D is not an artifact of
the simulation methods but instead comes from the granular packing structure.

###  IV.2 2D Simulations

In 2D and 3D, we perform similar simulations. However, the initial
configurations of grains are not along a line but are generated by the packing
procedure described above in Sec. III.3. Before describing the results of the
phase speed and attenuation measurements, we first discuss some features of
the packings themselves.

A common way to characterize granular packings involves characterizing the
number of contacts per particle $Z$ using the concept of isostaticity
van2009jamming . When a system is isostatic, the number of degrees of freedom
is equal to the number of constraints. For $N$ frictionless particles in 2D,
there are $2N$ degrees of freedom. Thus, isostaticity for frictionless disks
in 2D requires $2N$ total grain-grain contacts, or $Z=4$ contacts per
particle, since each contact is shared between two grains. This relation is
valid for large systems where constraints related to the boundaries are not
important. The isostatic limit is denoted $Z_{c}$. Figure 5(a) shows that
$\Delta Z\equiv Z-Z_{c}\propto\hat{P}^{1/2}$, meaning that $Z\rightarrow
Z_{c}$, with $Z_{c}=4$ in this case, as $\hat{P}\rightarrow 0$; this agrees
with many previous studies ohern2003jamming ; van2009jamming ;
thompson2019critical . We find similar results for 3D packings, but with
$Z_{c}=6$ due to the extra degree of freedom for each particle.

(a)

(b)

Figure 5: (a) Excess contacts $Z-Z_{c}$ plotted as a function of dimensionless
pressure $\hat{P}$, showing $Z-Z_{c}\propto\hat{P}^{1/2}$ (solid line). (b)
The probability distribution $p(\hat{F})$ of dimensionless contact forces
$\hat{F}$. Changing colors (blue to red) represent increasing $\hat{P}$, with
the same color and symbol convension between panels (a) and (b).

For sufficiently low $\hat{P}$ and $\Delta Z$, the statistical properties of
the granular force network become nearly independent of pressure. For example,
Fig. 5(b) shows a probability distribution $p(\hat{F})$ of dimensionless
contact force $\hat{F}\equiv F_{c}/Pd$ (note that in 2D pressure has units of
force per length). These force distributions show exponential tails with some
dependence on pressure. However, the distributions become independent of
$\hat{P}$ for $\hat{P}\leq 0.005$. For the high-pressure packings, which have
large overlaps between particles, there are extra contacts in the force
network leading to a narrower distribution of $\hat{F}$. These packings with
large $\hat{P}$ are certainly not physically similar to marine sediments, due
to large particle deformations, and are also unphysical in the sense that real
grains would deform and fracture at such high strains.

Using these packings, we vibrate the wall at $x=0$ and measure the resulting
displacements of all grains as a function of time. A plot of the $x$
displacement of a single grain is shown in Fig. 6(a). For every grain, we fit
to a form $\Delta x/A=\mathcal{A}\sin(\omega t-\Delta\phi)$ and extract the
normalized oscillation amplitude $\mathcal{A}$ and the phase offset
$\Delta\phi$, as shown in the figure. These are plotted as a function of the
$x$ position along the channel in panels (b) and (c) of Fig. 6. The slope of
$\Delta\phi$ versus $x(t=0)$ gives us the wavenumber $k$, and the magnitude of
the slope of $\log\mathcal{A}$ versus $x(t=0)$ gives us $\alpha$. We perform
such measurements for multiple packings per pressure, widely varying $\omega$,
$\gamma$, $\kappa$, and $m$. We then can obtain the dimensionless quantities
$\hat{\alpha}$ and $\hat{k}$ as a function of the input parameters
$\hat{\omega}$ and $\hat{\gamma}$ to test the agreement with the continuum
theory.

(a)

(b) (c)

Figure 6: (a) Normalized displacement of an individual grain, $\Delta x/A$
(solid black curve) is plotted as a function of time, which can be fit to an
sinusoid of the form $\mathcal{A}\sin(\omega t-\Delta\phi)$ (dashed red
curve). (b,c) The constants $\mathcal{A}$ (b) and $\Delta\phi$ (c) can then be
plotted as a function of the initial $x$ position of each grain and then used
to extract the wavenumber $k$ and the attenuation coefficient $\alpha$.

Figure 7 shows the key result of this paper, which are measurements of sound
speed and attenuation for the packings characterized in Fig. 5, using only
linear springs and dashpots at each contact. In accordance with the 1D
continuum theory in Sec. II, we plot $\hat{c}$ and $\hat{\alpha}/\hat{\omega}$
as a function of dimensionless frequency $\hat{\omega}\hat{\gamma}$ for
packings with pressures ranging from $\hat{P}=0.001$ to 0.1. Each data point
represents an average measurement over multiple realizations of packings at
the same $\hat{P}$ as well as multiple combinations of parameter values to
form the dimensionless quantity $\hat{\omega}\hat{\gamma}$. We do this to
check that the particular combinations of dimensionless quantities being
plotted are still valid (despite the disagreement with the continuum theory
from which they follow). For this reason, we show error bars using the
standard deviation of the measurements being averaged (not the standard error,
which would be the same quantity divided by the square root of the number of
measurements). Most of the error bars are smaller than the symbols themselves,
except for the measurements of very small $\alpha$ (since very little
attenuation occurs over the entire channel) and large $c$ (since these are
associated with very high driving frequencies and large attenuation). This
shows that for multiple parameter combinations, $\hat{c}$ and
$\hat{\alpha}/\hat{\omega}$ versus $\hat{\omega}\hat{\gamma}$ is still able to
collapse the data, albeit with different scaling forms for different values of
$\hat{P}$.

(a)

(b)

Figure 7: Dispersion and attenuation in 2D with varied pressure for
compression waves. Plots show $\hat{c}$ and $\hat{\alpha}/\hat{\omega}$ versus
$\hat{\omega}\hat{\gamma}$. Solid black lines show the theoretical result from
Eqs. (4) and (5).

The key message from these plots is that the functional forms of $c(f)$ and
$\alpha(f)$ for the low-pressure (i.e., more realistic) packings do not agree
with the continuum picture but instead follow a new, emergent functional form.
For large $\hat{P}$, we recover scaling functions for $\hat{c}$ and
$\hat{\alpha}/\hat{\omega}$ as a function of $\hat{\omega}\hat{\gamma}$ that
agree with the continuum picture (solid black lines), e.g.,
$\hat{\alpha}/\hat{\omega}\propto\hat{\omega}\hat{\gamma}$ at low frequencies.
However, as $\hat{P}$ is decreased, we observe that $\hat{c}$ and
$\hat{\alpha}/\hat{\omega}$ asymptotically approach new scaling forms. In
particular,
$\hat{\alpha}/\hat{\omega}\propto(\hat{\omega}\hat{\gamma})^{1/3}$, which
implies $\alpha\propto f^{4/3}$ for small $\hat{\omega}\hat{\gamma}$ in 2D.
Additionally, rather than $\hat{c}$ constant for
$\hat{\omega}\hat{\gamma}<0.1$, $\hat{c}$ weakly increases over
$10^{-3}<\hat{\omega}\hat{\gamma}<10^{-1}$ before joining the high-pressure
curves at $\hat{\omega}\hat{\gamma}\approx 0.5$ and above. In measurements of
sediment dispersion williams2002comparison , there typically exists a finite
low-frequency asymptote, which is related to the inertial effects of the
saturating fluid, as discussed by WilliamsWilliams2009b . Further model
development using a saturating fluid is key for explaining the low frequency
sound speed behavior.

Figure 8 shows similar measurements but with shear waves instead of
compressional waves. For these measurements, we vibrate the wall as before but
in the orthogonal lateral direction, tracking the lateral motion of particles
down the channel to obtain $\alpha$ and $k$. We then plot the same
dimensionless combinations as before. At high $\hat{P}$, the curves for speed
and attenuation mimic the continuum theory, albeit shifted in both magnitude
as well as the crossover value of $\hat{\omega}\hat{\gamma}$. As $\hat{P}$ is
decreased, the curves shift to lower $\hat{c}$ with less of a plateau at low
$\hat{\omega}\hat{\gamma}$ as well as an emergent scaling of
$\hat{\alpha}/\hat{\omega}$ roughly constant for small
$\hat{\omega}\hat{\gamma}$. This implies that $\alpha\propto f$ for shear
waves in this system.

(a)

(b)

Figure 8: Dispersion and attenuation in 2D with varied pressure for shear
waves. Plots show $\hat{c}$ and $\hat{\alpha}/\hat{\omega}$ versus
$\hat{\omega}\hat{\gamma}$. Solid black lines show the theoretical result from
Eqs. (4) and (5).

These emergent scaling forms likely arise from non-affine motion of grains
during the local deformations associated with the traveling waves. In this
context, “non-affine” refers to motion of individual grains that is not
perfectly aligned with large-scale, average motion utter2008experimental . We
support this claim with two figures. First, in Fig. 9, we show that grains
undergo elliptical motion with semi-major and semi-minor axes $a$ and $b$.
This is in contrast to linear (affine) motion, which would be expected from a
continuum picture with compressive oscillation in the $x$ direction.
Additionally, we find that these ellipses have increased aspect ratio $b/a$
and increased rotation angle $\theta$ as pressure is decreased, with
asymptotic behavior at low pressures (similar to the results in Figs. 7 and
8).

(a)

(b)

(c)

Figure 9: (a) A scatter plot of a single particle’s trajectory during a
simulation (red circle), along with a fitted ellipse (black line) with semi-
major axis $a$, semi-minor axis $b$, and rotation angle $\theta$. This
particle is typical of all particles, which are well-described by elliptical
trajectories. (b,c) Probability distributions $p(b/a)$ and $p(|\theta|)$ are
shown for varying $\hat{P}$ but all other parameters fixed, showing a trend
toward larger aspect ratio and larger rotation angle as pressure is decreased,
with asymptotic behavior for small $\hat{P}$.

Figure 9(a) shows a plot of the deviation of the particles from their initial
position for one grain during one simulation; the data points (red circles)
trace out an elliptical shape to excellent approximation. The black line is a
best-fit ellipse, with semi-major axis $a$, semi-minor axis $b$, and rotation
angle $\theta$, defined as the angle between the $x$ direction and the
direction of the closest semi-major axis. This particular plot is typical of
all grains, i.e., all grains trace out shapes that can be fit to an ellipse
with excellent fidelity.

Figure 9(b) and (c) show probability distributions $p$ of the aspect ratio
$b/a$ and the absolute value of the rotation angle $|\theta|$, since the
statistics of $\theta$ are even about $\theta=0$. For $\hat{P}=0.1$, both
distributions decay rapidly away from 0, suggesting again that at large
$\hat{P}$, the motion of grains is more similar to the continuum picture
(nearly linear along the direction of the applied oscillation at the
boundary). However, as $\hat{P}$ is decreased, the tails of the distribution,
which appear quasi-exponential, decay much more slowly, meaning there are many
more particles tracing out elliptical trajectories with larger values of $b/a$
and $|\theta|$. For $\hat{P}\leq 0.005$, the probability distributions become
independent of $\hat{P}$, similar to the dispersion and attenuation curves
shown previously. Thus, we conclude that the anomalous scaling shown in Figs.
7 and 8 is directly associated with a transition from predominantly linear
motion to elliptical motion.

The second figure to support our claim is Fig. 10. Here we show dispersion and
attenuation as a function of the product of dimensionless damping and
frequency for simulations where we damp only the component of relative
velocity that is along the grain-grain contact, i.e., we do not damp the
tangential motion. This is equivalent to using $\gamma_{n}=\gamma$ and
$\gamma_{t}=0$, as discussed in Sec. III. In such a system, we observe that
the attenuation coefficient mimics the 1D theory, albeit shifted and scaled.
Thus, the anomalous scaling arises when two conditions are met: first,
particle oscillations have significant non-affine components that are not
along the imposed oscillation direction, and, second, these non-affine motions
are damped.

Figure 10: Attenuation in 2D with varied pressure for compressional waves, but
where the tangential motion of grains is not damped. Curves are shifted but
remain parallel to the 1D continuum theory (solid black line, from Eqs. (4)
and (5)).

###  IV.3 3D Simulations

We repeat our 2D simulations for compressional waves in 3D. We find similar
results for the packing statistics as those shown in Fig. 5. We again find
$Z-Z_{c}\propto\hat{P}^{1/2}$ but with $Z_{c}=6$ for frictionless 3D spheres,
as expected van2009jamming . We again find that particles trace out elliptical
paths, as shown in Fig. 11. We again find that, as $\hat{P}$ is decreased, the
ellipses become less aligned with the imposed oscillation at the boundary and
they develop larger aspect ratios, although the summary statistics are not
shown here. We save a complete analysis of particle motion, including spatial
correlations, in 3D for future work.

Figure 11: A scatter plot of a single particle’s trajectory during a 3D
simulation. This particle is typical of all particles, which are well-
described by elliptical trajectories.

We again consider the motion of each grain along the imposed oscillation
direction and fit it to a sine wave, as in Fig. 6. We measure $k$ and $\alpha$
in the same way, and again plot $\hat{c}$ and $\hat{\alpha}/\hat{\omega}$ as a
function of $\hat{\omega}\hat{\gamma}$ for a range of parameter values, as
shown in Fig. 12. We again find that $\hat{\omega}\hat{\gamma}$ is a good
predictor of these two quantities, irrespective of the values of the
individual parameters used to construct this combination, at least for an
individual value of $\hat{P}$.

(a)

(b)

Figure 12: Dispersion and attenuation in 3D with varied pressure for
compression waves. Plots show $\hat{c}$ and $\hat{\alpha}/\hat{\omega}$ versus
$\hat{\omega}\hat{\gamma}$. Solid black lines show the theoretical result from
Eqs. (4) and (5).

The results shown in Fig. 12 are very similar to the 2D case, except for the
fact that $\hat{\alpha}/\hat{\omega}\propto(\hat{\omega}\hat{\gamma})^{\beta}$
with $\beta\approx 1/6$ in 3D, as opposed to $\beta\approx 1/3$ in 2D. This
implies $\alpha\propto f^{7/6}$, which is very similar to $\alpha\propto f$.

##  V Discussion

In this manuscript we have shown that the packing structure of the granular
material may play an important role in controlling the acoustic properties of
marine sediments. By studying frictionless granular packings using only
linearized forces (springs and dashpots), we observe dispersion and
attenuation curves that do not agree with a naive continuum approach, which
predicts $\alpha\propto f^{2}$ for low to moderate frequencies using linear
forces. Instead, we observe that the granular packing structure itself leads
to emergent scaling forms with $\alpha\propto f^{\beta}$, where $\beta\approx
4/3$ for packings of 2D frictionless disks and $\beta\approx 7/6$ for packings
of 3D frictionless spheres. The observed scaling laws are outputs of the DEM
simulations and are not currently explained by any theory.

Thus, the key message of our paper is to challenge a common assumption in
prior theoretical approaches, that linearized forces necessarily lead to
dispersion and attenuation scaling that agrees with classic viscosity and
Biot-Stoll theory, notably $\alpha\propto f^{2}$ at low $f$. Finding a minimal
theoretical framework that would give rise to other scaling laws, e.g.,
$\alpha\propto f$, was a major factor in motivating the GS and VGS models
Buckingham1997 ; buckingham_wave_2000 . In contrast, our results show that
linear forces can give rise to alternative scaling laws that come directly
from the granular packing structure. We observe a transition from
$\alpha\propto f^{2}$ for unphysical, high-pressure packings (with large
overlaps between particles) to $\alpha\propto\beta$ with $\beta\approx 4/3$
for P-waves in 2D, $\beta\approx 1$ for S-waves in 2D, and $\beta\approx 7/6$
for P-waves in 3D. This transition occurs as the packings asymptotically
approach the (more realistic) stiff-particle limit. This transition is also
associated with changes in the force statistics of the packing structure (Fig.
5) as well as a transition from linear to ellipitcal motion of individual
grains (Fig. 9).

Our intent is not necessarily to advocate for any particular set of
experimental data. However, we note that the most direct measurements of
acoustic properties of an idealized granular packing immersed in a viscous
fluid are the laboratory experiments by Hefner and Williams hefner2006sound .
The fact that our simulation results are similar to these measurements is very
encouraging. These and other results were shown above in Fig. 1.

There are a large number of questions that arise from our work. One line of
inquiry might explore how our results would change using explicitly nonlinear
grain-grain forces, like Hertzian contact mechanics or the history-dependent
dashpot terms from the GS or VGS models. Additionally, more realistic granular
packings, like those arising from frictional grain-grain interactions or
irregular grain shapes, might give rise to changes in the scaling laws we
observe. Grain-grain friction, polydispersity (i.e., a larger range of grain
sizes), and irregular grain shape are known to affect the properties of the
intergrain contact network silbert2010jamming ; nguyen2014effect ;
nguyen2015effects .

Future work may also try to further characterize the non-affine grain motions
that are associated with the anomalous scaling behavior. For example, we have
shown distributions of the ellipses traced out by grains in 2D, but have not
considered whether there are spatial correlations in the orientation and
aspect ratio of these ellipses. This may provide additional insight, since it
is not the elliptical motion directly that increases dissipation but the
relative motion between each grain and its neighbors.

Finally, the effect of a saturating fluid beyond its contribution to the
dashpot magnitude was not included in our work. Previous workWilliams2009b
has shown that inertial effects of the fluid and grains must be treated
properly. Including the inertia, as well as compressional and viscous effects
of the saturating fluid is a key future step in relating the work presented
here to marine sediment dispersion.

## Conflict of Interest

The authors declare no conflicts of interest.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

###### Acknowledgements.

The authors thank Dr. Charles Holland for discussions regarding measurements
of sediment dispersion. Funding was proved by the U.S. Office of Naval
Research, for which the authors are grateful. Any opinions, findings, and
conclusions or recommendations expressed in this publication are those of the
author and do not necessarily reflect the views of the U.S. Navy.

## References

  * (1) M. S. Ballard and K. Lee, “The acoustics of marine sediments,” Acoust. Today 13(3), 11–18 (2017).
  * (2) C. Daraio, V. Nesterenko, E. Herbold, and S. Jin, “Strongly nonlinear waves in a chain of teflon beads,” Physical Review E 72(1), 016603 (2005).
  * (3) L. R. Gómez, A. M. Turner, M. van Hecke, and V. Vitelli, “Shocks near jamming,” Physical review letters 108(5), 058001 (2012).
  * (4) J. Shäfer, S. Dippel, and D. Wolf, “Force schemes in simulations of granular materials,” Journal de Physique I 6 (1996) 10.1051/jp1:1996129.
  * (5) T. S. Majmudar and R. P. Behringer, “Contact force measurements and stress-induced anisotropy in granular materials,” nature 435(7045), 1079–1082 (2005).
  * (6) L. E. Silbert, A. J. Liu, and S. R. Nagel, “Vibrations and diverging length scales near the unjamming transition,” Phys. Rev. Lett. 95, 098301 (2005) https://link.aps.org/doi/10.1103/PhysRevLett.95.098301 10.1103/PhysRevLett.95.098301.
  * (7) N. P. Chotiros, _Acoustics of the Seabed as a Poroelastic Medium_ (Springer, 2017).
  * (8) K. L. Williams, “Sand acoustics: The effective density fluid model, pierce/carey expressions, and inferences for porous media modeling,” The Journal of the Acoustical Society of America 125(4), EL164–EL170 (2009) 10.1121/1.3097681.
  * (9) C. W. Holland, P. L. Nielsen, J. Dettmer, and S. Dosso, “Resolving meso-scale seabed variability using reflection measurements from an autonomous underwater vehicle,” The Journal of the Acoustical Society of America 131(2), 1066–1078 (2012) 10.1121/1.3672696.
  * (10) J. Dettmer, S. E. Dosso, and C. W. Holland, “Trans-dimensional geoacoustic inversion,” The Journal of the Acoustical Society of America 128(6), 3393–3405 (2010) 10.1121/1.3500674.
  * (11) J. E. Quijano, S. E. Dosso, J. Dettmer, L. M. Zurk, M. Siderius, and C. H. Harrison, “Bayesian geoacoustic inversion using wind-driven ambient noise,” The Journal of the Acoustical Society of America 131(4), 2658–2667 (2012) 10.1121/1.3688482.
  * (12) C. W. Holland, S. Pinson, C. M. Smith, P. C. Hines, D. R. Olson, S. E. Dosso, and J. Dettmer, “Seabed structure inferences from TREX13 reflection measurements,” IEEE Journal of Oceanic Engineering 42(2), 268–288 (2017) 10.1109/joe.2017.2658418.
  * (13) D. R. Olson, “The effect of seafloor roughness on passive estimates of the seabed reflection coefficient,” The Journal of the Acoustical Society of America 153(1), 586–601 (2023) 10.1121/10.0016846.
  * (14) J. Bonnel, S. E. Dosso, and N. Ross Chapman, “Bayesian geoacoustic inversion of single hydrophone light bulb data using warping dispersion analysis,” The Journal of the Acoustical Society of America 134(1), 120–130 (2013) 10.1121/1.4809678.
  * (15) T. W. Tan and O. A. Godin, “Ambient sound directionality and rapid estimation of empirical Green’s functions in a coastal ocean,” The Journal of the Acoustical Society of America 152(4), A153 (2022) https://doi.org/10.1121/10.0015861 10.1121/10.0015861.
  * (16) C. W. Holland and J. Dettmer, “In situ sediment dispersion estimates in the presence of discrete layers and gradients,” The Journal of the Acoustical Society of America 133(1), 50–61 (2013) 10.1121/1.4765300.
  * (17) M. B. Porter and E. L. Reiss, “A numerical method for bottom interacting ocean acoustic normal modes,” The Journal of the Acoustical Society of America 77(5), 1760–1767 (1985) 10.1121/1.391925.
  * (18) C. W. Holland and P. Neumann, “Sub-bottom scattering: A modeling approach,” J. Acoust. Soc. America 104(3), 1363–1373 (1998) 10.1121/1.424345.
  * (19) K. D. LePage and H. Schmidt, “Spectral integral representations of volume scattering in sediments in layered waveguides,” The Journal of the Acoustical Society of America 108(4), 1557–1567 (2000) 10.1121/1.1289370.
  * (20) K. D. LePage and H. Schmidt, “Spectral integral representations of monostatic backscattering from three-dimensional distributions of sediment volume inhomogeneities,” J. Acoust. Soc. America 113(2), 789–799 (2003) 10.1121/1.1528625.
  * (21) D. Tang and D. Jackson, “A time-domain model for seafloor scattering,” The Journal of the Acoustical Society of America 142(5), 2968–2978 (2017) 10.1121/1.5009932.
  * (22) D. R. Olson and C. W. Holland, “Fast computation of time-domain scattering by an inhomogeneous stratified seafloor,” The Journal of the Acoustical Society of America 147(1), 191–204 (2020) 10.1121/10.0000570.
  * (23) D. Jackson and D. R. Olson, “The small-slope approximation for layered, fluid seafloors,” The Journal of the Acoustical Society of America 147(1), 56–73 (2020) https://doi.org/10.1121/10.0000470 10.1121/10.0000470.
  * (24) D. R. Olson and D. Jackson, “Scattering from layered seafloors: Comparisons between theory and integral equations,” The Journal of the Acoustical Society of America 148(4), 2086–2095 (2020) 10.1121/10.0002164.
  * (25) D. A. Abraham, _Underwater Acoustic Signal Processing_ (Springer International Publishing, 2019).
  * (26) J.-X. Zhou, X.-Z. Zhang, and D. P. Knobles, “Low-frequency geoacoustic model for the effective properties of sandy seabottoms,” The Journal of the Acoustical Society of America 125(5), 2847–2866 (2009) https://doi.org/10.1121/1.3089218 10.1121/1.3089218 publisher: Acoustical Society of America.
  * (27) M. Biot, “Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range,” The Journal of the Acoustical Society of America (1956) https://pubs.aip.org/asa/jasa/article/28/2/168/737238/Theory-of-Propagation-of-Elastic-Waves-in-a-Fluid.
  * (28) M. A. Biot, “Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher frequency range,” The Journal of the Acoustical Society of America 28(2), 179–191 (1955) https://doi.org/10.1121/1.1908241 10.1121/1.1908241.
  * (29) M. A. Biot, “Generalized theory of acoustic propagation in porous dissipative media,” The Journal of the Acoustical Society of America 34(9A), 1254–1264 (1962) https://doi.org/10.1121/1.1918315 10.1121/1.1918315.
  * (30) R. D. Stoll, “Acoustic waves in ocean sediments,” GEOPHYSICS 42(4), 715–725 (1977) https://library.seg.org/doi/10.1190/1.1440741 10.1190/1.1440741 publisher: Society of Exploration Geophysicists.
  * (31) R. D. Stoll, “Theoretical aspects of sound transmission in sediments,” The Journal of the Acoustical Society of America 68(5), 1341–1350 (1980) https://doi.org/10.1121/1.385101 10.1121/1.385101.
  * (32) R. D. Stoll and T. K. Kan, “Reflection of acoustic waves at a water–sediment interface,” The Journal of the Acoustical Society of America 70(1), 149–156 (1981) https://doi.org/10.1121/1.386692 10.1121/1.386692.
  * (33) R. D. Stoll, “Marine sediment acoustics,” The Journal of the Acoustical Society of America 77(5), 1789–1799 (1985) https://doi.org/10.1121/1.391928 10.1121/1.391928.
  * (34) O. A. Godin, “Shear waves and sound attenuation in underwater waveguides,” The Journal of the Acoustical Society of America 149(5), 3586–3598 (2021) 10.1121/10.0004999.
  * (35) O. A. Godin, “Effects of weak shear rigidity of the seabed on sound propagation in a range-dependent ocean,” The Journal of the Acoustical Society of America 153(3_supplement), A86 (2023) https://doi.org/10.1121/10.0018259 10.1121/10.0018259.
  * (36) E. L. Hamilton, “Geoacoustic modeling of the sea floor,” J. Acoust. Soc. America 68(5), 1313–1340 (1980) 10.1121/1.385100.
  * (37) K. L. Williams, D. R. Jackson, E. I. Thorsos, D. Tang, and S. G. Schock, “Comparison of sound speed and attenuation measured in a sandy sediment to predictions based on the biot theory of porous media,” IEEE journal of oceanic engineering 27(3), 413–428 (2002).
  * (38) B. Hefner, D. Jackson, K. Williams, and E. Thorsos, “Mid- to high-frequency acoustic penetration and propagation measurements in a sandy sediment,” IEEE Journal of Oceanic Engineering 34(4), 372–387 (2009) 10.1109/joe.2009.2028410.
  * (39) A. Turgut and T. Yamamoto, “Measurements of acoustic wave velocities and attenuation in marine sediments,” The Journal of the Acoustical Society of America 87(6), 2376–2383 (1990) 10.1121/1.399084.
  * (40) H. J. Simpson, B. H. Houston, S. W. Liskey, P. A. Frank, A. R. Berdoz, L. A. Kraus, C. K. Frederickson, and S. Stanic, “At-sea measurements of sound penetration into sediments using a buried vertical synthetic array,” The Journal of the Acoustical Society of America 114(3), 1281–1290 (2003) 10.1121/1.1594192.
  * (41) D. J. Wingham, “The dispersion of sound in sediment,” The Journal of the Acoustical Society of America 78(5), 1757–1760 (1985) 10.1121/1.392761.
  * (42) B. T. Hefner and K. L. Williams, “Sound speed and attenuation measurements in unconsolidated glass-bead sediments saturated with viscous pore fluids,” The Journal of the Acoustical Society of America 120(5), 2538–2549 (2006).
  * (43) K. I. Lee, V. F. Humphrey, B.-N. Kim, and S. W. Yoon, “Frequency dependencies of phase velocity and attenuation coefficient in a water-saturated sandy sediment from 0.3to1.0mhz,” The Journal of the Acoustical Society of America 121(5), 2553–2558 (2007) 10.1121/1.2713690.
  * (44) J.-P. Sessarego, A. Ivakin, and D. Ferrand, “Frequency dependence of phase speed, group speed, and attenuation in water-saturated sand: Laboratory experiments,” IEEE Journal of Oceanic Engineering 33(4), 359–366 (2008) 10.1109/joe.2008.927584.
  * (45) J. Hare and A. E. Hay, “Phase speed in water-saturated sand and glass beads at MHz frequencies,” The Journal of the Acoustical Society of America 148(4), 2301–2310 (2020) 10.1121/10.0002250.
  * (46) M. J. Buckingham, “Theory of acoustic attenuation, dispersion, and pulse propagation in unconsolidated granular materials including marine sediments,” The Journal of the Acoustical Society of America 102(5), 2579–2596 (1997) 10.1121/1.420313.
  * (47) M. J. Buckingham, “Wave propagation, stress relaxation, and grain-to-grain shearing in saturated, unconsolidated marine sediments,” The Journal of the Acoustical Society of America 108(6), 2796–2815 (2000) https://doi.org/10.1121/1.1322018 10.1121/1.1322018.
  * (48) M. J. Buckingham, “On pore-fluid viscosity and the wave properties of saturated granular materials including marine sediments,” The Journal of the Acoustical Society of America 108(6) (2007) https://pubs.aip.org/asa/jasa/article/122/3/1486/853009/On-pore-fluid-viscosity-and-the-wave-properties-of.
  * (49) H. M. Jaeger, S. R. Nagel, and R. P. Behringer, “Granular solids, liquids, and gases,” Reviews of modern physics 68(4), 1259 (1996).
  * (50) Y. Forterre and O. Pouliquen, “Flows of dense granular media,” Annu. Rev. Fluid Mech. 40, 1–24 (2008).
  * (51) M. van Hecke, “Jamming of soft particles: Geometry, mechanics, scaling and isostaticity,” Journal of Physics: Condensed Matter 22(3), 033101 (2009).
  * (52) A. J. Liu and S. R. Nagel, “The jamming transition and the marginally jammed solid,” Annu. Rev. Condens. Matter Phys. 1(1), 347–369 (2010).
  * (53) R. P. Behringer and B. Chakraborty, “The physics of jamming for granular materials: a review,” Reports on Progress in Physics 82(1), 012601 (2018).
  * (54) A. Amon, P. Born, K. E. Daniels, J. A. Dijksman, K. Huang, D. Parker, M. Schröter, R. Stannarius, and A. Wierschem, “Preface: Focus on imaging methods in granular physics,” Review of Scientific Instruments 88(5) (2017).
  * (55) F. Radjai, J.-N. Roux, and A. Daouadji, “Modeling granular materials: century-long research across scales,” Journal of engineering mechanics 143(4), 04017002 (2017).
  * (56) K. Kamrin, K. M. Hill, D. I. Goldman, and J. E. Andrade, “Advances in modeling dense granular media,” Annual Review of Fluid Mechanics 56, 215–240 (2024).
  * (57) M. Wyart, L. E. Silbert, S. R. Nagel, and T. A. Witten, “Effects of compression on the vibrational modes of marginally jammed solids,” Phys. Rev. E 72, 051306 (2005) https://link.aps.org/doi/10.1103/PhysRevE.72.051306 10.1103/PhysRevE.72.051306.
  * (58) K. Saitoh, K. Taghizadeh, and S. Luding, “Sound characteristics of disordered granular disks: effects of contact damping,” Frontiers in Physics 11, 1192270 (2023).
  * (59) E. T. Owens and K. E. Daniels, “Sound propagation and force chains in granular materials,” Europhysics Letters 94(5), 54005 (2011) https://dx.doi.org/10.1209/0295-5075/94/54005 10.1209/0295-5075/94/54005.
  * (60) C. F. Schreck, T. Bertrand, C. S. O’Hern, and M. D. Shattuck, “Repulsive contact interactions make jammed particulate systems inherently nonharmonic,” Phys. Rev. Lett. 107, 078301 (2011) https://link.aps.org/doi/10.1103/PhysRevLett.107.078301 10.1103/PhysRevLett.107.078301.
  * (61) L. E. Silbert, “Jamming of frictional spheres and random loose packing,” Soft Matter 6(13), 2918–2924 (2010).
  * (62) D.-H. Nguyen, É. Azéma, F. Radjai, and P. Sornay, “Effect of size polydispersity versus particle shape in dense granular media,” Physical Review E 90(1), 012202 (2014).
  * (63) D.-H. Nguyen, E. Azéma, P. Sornay, and F. Radjai, “Effects of shape and size polydispersity on strength properties of granular materials,” Physical Review E 91(3), 032203 (2015).
  * (64) C. Coste and B. Gilles, “On the validity of hertz contact law for granular material acoustics,” The European Physical Journal B-Condensed Matter and Complex Systems 7, 155–168 (1999).
  * (65) S. L. Garrett, _Understanding Acoustics_ (Springer International Publishing, 2020).
  * (66) M. V. Shitikova and A. I. Krusser, “Models of viscoelastic materials: a review on historical development and formulation,” Theoretical Analyses, Computations, and Experiments of Multiscale Materials: A Tribute to Francesco dell’Isola 285–326 (2022).
  * (67) K. Binder, J. Horbach, W. Kob, W. Paul, and F. Varnik, “Molecular dynamics simulations,” Journal of Physics: Condensed Matter 16(5), S429 (2004).
  * (68) M. Montaine, M. Heckel, C. Kruelle, T. Schwager, and T. Pöschel, “Coefficient of restitution as a fluctuating quantity,” Physical Review E 84(4), 041306 (2011) 10.1103/physreve.84.041306.
  * (69) C. S. O’hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, “Jamming at zero temperature and zero applied stress: The epitome of disorder,” Physical Review E 68(1), 011306 (2003).
  * (70) J. D. Thompson and A. H. Clark, “Critical scaling for yield is independent of distance to isostaticity,” Physical Review Research 1(1), 012002 (2019).
  * (71) B. Utter and R. Behringer, “Experimental measures of affine and nonaffine deformation in granular shear,” Physical review letters 100(20), 208302 (2008).