Hydro-micromechanical modeling of wave propagation in saturated granular media

Biot's theory predicts the wave velocities of a saturated poroelastic granular medium from the elastic properties, density and geometry of its dry solid matrix and the pore fluid, neglecting the interaction between constituent particles and local flow. However, when the frequencies become high and the wavelengths comparable with particle size, the details of the microstructure start to play an important role. Here, a novel hydro-micromechanical numerical model is proposed by coupling the lattice Boltzmann method (LBM) with the discrete element method (DEM. The model allows to investigate the details of the particle-fluid interaction during propagation of elastic waves While the DEM is tracking the translational and rotational motion of each solid particle, the LBM can resolve the pore-scale hydrodynamics. Solid and fluid phases are two-way coupled through momentum exchange. The coupling scheme is benchmarked with the terminal velocity of a single sphere settling in a fluid. To mimic a pressure wave entering a saturated granular medium, an oscillating pressure boundary condition on the fluid is implemented and benchmarked with one-dimensional wave equations. Using a face centered cubic structure, the effects of input waveforms and frequencies on the dispersion relations are investigated. Finally, the wave velocities at various effective confining pressures predicted by the numerical model are compared with with Biot's analytical solution, and a very good agreement is found. In addition to the pressure and shear waves, slow compressional waves are observed in the simulations, as predicted by Biot's theory.


INTRODUCTION
Understanding wave propagation in dry and saturated granular media is vital for nondestructive soil testing, 1-3 seismicity analysis, [4][5][6] and oil exploration, 7,8 among others. Extensive work has been done to study the dispersion and attenuation properties of dry granular media in both laboratory experiments and computer simulations. [9][10][11][12] The kinematics at the microscale associated with wave propagation in dry granular media can be reproduced with the discrete element method (DEM), 11,[13][14][15] given that the parameters are appropriately calibrated. 16,17 To simulate wave propagation in fully/partially saturated granular media, 18 the momentum and energy exchange between fluid, solid, and/or gas phases must be taken into account in a fully coupled manner.
Depending on the spatial and temporal resolution of hydrodynamic interactions in the pore fluid, various fluid-solid coupling methods have been proposed for dense 19 or loose granular materials 20,21 submerged in viscous fluids. While conventional computational fluid dynamics (CFD) methods couple the drag forces on solid particles locally into the Navier-Stokes equations, 22 direct numerical simulation (DNS) techniques such as the lattice Boltzmann method (LBM) 23,24 resolve hydrodynamic interactions at the pore scale and maintain the momentum balance at no-slip fluid-solid interfaces. [25][26][27][28] Conventional CFD methods are sufficient for modeling dilute suspensions of particles in which local pore-scale fluid flow is not essential, whereas DNS techniques are better suited for fully/partially saturated dense granular media. Most DNS techniques (eg, coupled FEM-DEM) require adaptive remeshing around particles that move in the fluid, with high computational costs. Alternatively, the LBM projects particles on a Eulerian lattice to advect and bounce back the fluid and thus avoids adaptive remeshing. The main advantage of the LBM for fluid-solid coupling is the locality of the collision-streaming operations and the explicit time-stepping scheme, 29 which makes LBM simulations highly parallelizable. For geotechnical and geophysical applications in which the saturated granular materials are typically dense, an accurate prediction of the pressure gradients in the pore space is pivotal to the simulation of fluid-solid interaction processes such as consolidation 30,31 and wave propagation. 32,33 In the past decade, coupled LBM-DEM simulations have been applied to solve various geotechnical problems, including multiparticle sedimentation, 29 quicksand, 34 hydraulic fracture, 35 piping erosion, 36 underwater landslide, 37 and the collapse of a sand arch inside a perforation cavity. 38 Recent applications to acoustics have demonstrated the capability of the LBM to capture weak pressure fluctuations in pure fluids with high accuracy. In Viggen, 39 the acoustic behavior was reproduced with a monopole source, which was later extended to multipole in both 2D 40 and 3D. 41 In the present work, for the first time, two-way coupled LBM-DEM modeling is exploited to study the fluid-solid interaction at the pore scale associated with the propagation of pressure and shear waves in saturated poroelastic media. The porous material is reproduced as a regular crystal of identical solid spheres with a viscous fluid that fills the pore space. The LBM resolves the pore-scale fluid flow because of the propagation of elastic waves, while DEM describes the motion and interparticle forces of solid particles. The fluid-solid coupling is handled with the momentum exchange method (MEM), 24,42 which computes the hydrodynamic forces on solid particles and updates local flow at the no-slip fluid-solid interfaces. In order to propagate waves in the fluid-solid system, an oscillating pressure boundary condition that emits acoustic waves coming from the fluid is implemented.
The hydro-micromechanical mechanical model is tested in comparison with the continuum description of compressional (P)-wave and shear (S)-wave in saturated granular media, based on Biot theory, ie, a set of balance equations that macroscopically describe wave propagation in a two-phase poroelastic material. Biot theory predicts that the P-wave typically travels faster in a saturated granular medium than in its dry solid matrix and the S-wave velocity is reduced because of the inertia coupling due to the relative acceleration between the solid and the pore fluid. Given the continuum nature of the theory, both the pore and solid phases are assumed to be homogeneous and the influence of the microstructure on the fluid flow and effective elastic moduli is neglected. However, the effects of the microstructure and local flow, as well as their interaction, should come into play when the wavelength of the acoustic wave is comparable with the size of the constituent solid particles. Our coupled LBM-DEM code and the newly implemented oscillating pressure boundary allow for a comprehensive investigation of the dispersion relations and attenuation in both solid and fluid phases, resulting from the propagation of low-to-high frequency waves. A simplified fluid-solid system, namely, a face-centered-cubic (FCC) granular crystal, is adopted to perform LBM-DEM simulations at various effective confining pressures, and the numerical results are compared with Biot analytical solutions. Other governing mechanisms including the highly attenuated, diffusive slow wave as predicted by the theory are explored as well.
The remainder of this paper is organized as follows. Section 2 explains the fundamentals of the hydro-micromechanical model, including the LBM for simulating fluid in Section 2.1, the DEM for solid particles in Section 2.2, the fluid-solid interactions in Section 2.3, and the acoustic source in Section 2.4. Section 3 benchmarks the fluid-solid coupling scheme and the acoustic source with respect to analytical solutions. The hydro-micromechanical model is applied to simulate elastic wave propagation in fully saturated FCC granular crystals in Section 4, with the model parameters and unit conversion highlighted in Section 4.2. The influence of pore fluid on the dispersion relations is discussed in Section 4.3 and the effect of acoustic source investigated in Section 4.4. Section 5 shows the comparison between the numerical predictions and Biot analytical solutions, as well as the numerical evidence of the slow compressional wave. Conclusions are drawn in Section 6.

HYDRO-MICROMECHANICAL MODEL
The local hydrodynamics in the pore space of a granular medium is modeled by solving the discrete Boltzmann equation with a Bhatnagar-Gross-Krook (BGK) collision operator. The BGK operator uses a single relaxation time for the probability distribution of local fluid velocities towards an equilibrium distribution. Interactions between contacting solid particles and their microstructural evolution are modeled with the DEM. The fluid-solid coupling scheme is implemented following the link-based MEM. A novel oscillating pressure boundary condition is employed to send acoustic waves from the fluid phase into saturated granular media.

Hydrodynamics
The Boltzmann equation describes the evolution of a collection of particles via their associated probability density functions in space and velocity. 24 Interactions between these particles are accounted by a collision operator, which we take here to have the simple BGK form, such that, Here, f(x, , t) is the probability distribution function associated with particles at position x with velocity ; B is an external force and the relaxation time during which local distributions are relaxed to the equilibrium distribu-tioñ(x, , t). The macroscopic fluid density and velocity are retrieved via = ∫ Δ and u = (1∕ ) ∫ Δ . If the Maxwell-Boltzmann distribution is considered at equilibrium, with v = − u and p the pressure, it can be shown that through the Chapman-Enskog expansion, with the Knudsen number sufficiently small, the Navier-Stokes equations are recovered. The LBM solves the Boltzmann equation discretized in space and velocity. Having its origin in statistical mechanics, the variables in LBM are probability distribution functions on a Cartesian lattice rather than macroscopic pressure and velocity in conventional CFD methods. At each time step, fluid particles on each lattice node are streamed to their immediate neighboring nodes along certain directions. The collision operator is then directly computed at each lattice site from the macroscopic variables.
The discretization of the velocity space requires a finite set of velocity vectors along which the fluid particles can propagate and collide, reducing the continuous distribution function f(x, , t) to f d (x, t) on the discrete velocity directions. In this work, we use the D3Q19 lattice (three dimensions, 19 velocity vectors) with spacing Δx, such that fluid particles can propagate to the nearest and next-nearest neighbors at each time step Δt. The discrete velocity vectors c d (i = 1... 19) are defined as with c = Δx∕Δt, and illustrated in Figure 1. By adopting the discretizations above, the BGK collisional operator (1) can be rewritten as̃d wherẽd(x, t) are the postcollision values to be distributed to the corresponding lattice neighbors in the succeeding streaming step. The dimensionless relaxation time ∈ (0.5, ∞) implicitly determines the kinematic viscosity , together with the lattice spacing Δx and the time step Δt via = ( − Δt∕2)c 2 s . The lattice sound speed c s = Δx∕ √ 3Δt for the D3Q19 lattice is defined by the discretizations. After expansion to second order, the equilibrium distribution where the macroscopic fluid density = f + ′ , with the mean density f and the off-equilibrium density ′ . Note that ′ is proportional to the fluctuation of fluid pressure that is p ′ = c 2 s ′ . Together with the macroscopic fluid velocity, they are calculated via the zeroth and first moments of the distribution functions f d (x, t) as, The lattice weights w d in (6) for a D3Q19 LB model are given as follows.

Micromechanics
The DEM represents granular materials as packings of solid particles with simplified geometries (eg, spheres) and vanishingly small interparticle overlaps. 44 In this way, the particles in contact can interact with their neighbors via repulsive springs and viscous dashpots, resulting in relative motion between particles. Generally, the time step in DEM is very small in order to resolve collisions or oscillations in time sufficiently. Once all the forces acting on each particle, either from interacting neighbors or from the surrounding fluid, are known, the kinematics of each particle in all spatial directions are updated by Newton second law via the explicit time integration scheme. The equations of translational and rotational motion for each solid particle are m where m and I are the mass and the moment of inertia of the particle, with V p and p the resultant linear and angular velocities at the center X p . F c is a contact force applied on the particle, with c summing over the contact points N c p ; X c − X p is the vector connecting the contact point X c to the center position X p ; F f − s and T f − s are the hydrodynamic force and moment acting on the particle; and F g the external body force. The superposed dot represents a time derivative, eg, V p = . X p . The interparticle force between two contacting solid particles can be described by contact-level force-displacement laws in normal and tangential directions, ie, F c = F n + F t . The tangential component of the force is constrained by a Coulomb-type yield criterion. For two contacting spheres with a normal overlap u n and a relative tangential displacement Δu t , the interparticle normal and tangential forces F n and ΔF t can be calculated as where E p and p are the contact-level Young modulus and Poisson ratio, is the interparticle friction coefficient, and R * is the equivalent radius defined as 1∕(1∕R 1 + 1∕R 2 ), where R 1 and R 2 are the radii of the two particles. From the contact force network in a given microstructural configuration, the effective stress tensor is given by the average over the particle assembly where ′ is the bulk stress tensor, N c is the total number of contacts contained within the force network, V is the total volume of the granular assembly including the pore space, and d c is the branch vector connecting the centers of each pair of contacting particles.

Fluid-solid coupling scheme
Boundary conditions in the LBM are implemented by constructing the distribution functions from physical boundary constraints such as pressure and velocity. When coming to the fluid-solid coupling, two popular coupling schemes are the MEM and Noble-Torczynski method. The MEM was first presented by Ladd and Verberg 45 and Aidun et al 46 for handling moving boundary conditions between fluid and solid. The main idea is to project solid obstacles such as solid particles onto the lattice, marking the nodes contained by the particles as solid and the rest as fluid. A solid node has at least one link with neighboring fluid nodes. Along the links, as shown in Figure 2, momentum is exchanged and no-slip boundary conditions, ie, the macroscopic flow velocity matching the particle's velocity, are applied on the solid surface. From the exchanged momentum, the hydrodynamic force acting on the particle can be obtained as a consequence of the no-slip condition applied to fluid advecting towards solid nodes. As the solid particle moves across the lattice, fluid nodes will be converted to solid nodes and vice versa. Figure 2 shows an illustration of the coupling scheme. Alternatively, the coupling can be achieved by modifying the local collisional operator using the solid volume fractions of fluid cells occupied by solid particles. 47 Comparative studies 48,49 have shown that the Noble-Torczynski method causes numerical dissipation, with the accuracy depending on the relaxation time, whereas the MEM is prone to high-frequency fluctuation because of the binary transition between solid and fluid nodes. Nevertheless, the noise appears at scales much smaller than particle sizes and thus is smeared out after local averaging. [50][51][52][53] For this reason, the MEM has been selected as the fluid-solid coupling scheme for this work. Fluid advecting onto the surface of a solid particle is simply bounced back half way between the solid and the fluid nodes by setting where d * denotes the direction opposite to d such that is the velocity at the boundary position x b located midway between the fluid and the solid nodes. The second term of Equation 14, in addition to the standard bounce-back solid boundary, is contributed by the motion of the particle moving across the lattice. Following the momentum exchange idea, the force F d −s transferred to the solid particle at position x b over a fluid-solid link f-s from every bounce-back collision is It can be seen from Equation 15 that the hydrodynamic forces coupled to the solid particle are a direct consequence of the unbalanced momenta along the fluid-solid links.
As the particle moves and occupies neighboring fluid nodes, the fluid nodes at x f→s are converted to solid ones. During the conversion, the momentum on the fluid nodes has to be transferred to the particle, so that the global momentum is conserved. Therefore, a correction term is added to the fluid-solid coupling, The lattice nodes regarded as solid nodes do not participate in the LBM calculation cycles, namely, Equations 4 to 6. However, as the particle moves away from its solid nodes, the uncovered solid nodes at x s→f have to be converted back to the fluid phase. In order to ensure a consistent reconstruction of the probability distribution on each new fluid node, the velocity of the solid particle therein is used to calculate the equilibrium distributions in (6), and the new density corresponds to the average of neighboring nodes. Similar to Equation 16, the final correction that takes the additional momentum away from the particle is Therefore, the total force F f-s and torque T f − s acting on a solid particle are obtained from (1) the forces given by the momentum exchange, with the sum running over the fluid-solid links and then the boundary nodes, and (2) the forces resulting from the destruction and reconstruction of the relevant fluid nodes, namely,

Oscillating pressure boundary condition
A novel acoustic source is implemented to mimic pressure waves entering a fluid. The numerical approach consists in enforcing the boundary nodes to emit the correct numbers of fluid particles, which macroscopically leads to oscillations in the velocity component along the propagation direction. A pressure boundary condition is first needed. To model fluid flow in granular media, the pressure boundary conditions are typically set with constant local densities. The flow is then driven by an imposed pressure gradient between two opposite boundaries. This method was first suggested by Zou and He for the D2Q9 lattice 54 and later extended to the D3Q15 and D3Q19 lattice by Kutay et al 55 and Hecht and Harting. 56 At a pressure boundary, the unknown velocity component in the flow direction can be determined from the two known velocity components and a given density . For a D3Q19 LB model, this is done by enforcing the bounce-back rule to the nonequilibrium part of the distribution d − eq d at the boundary, namely, d − eq d = d * − eq d * . Therefore, for the pressure boundary located at the bottom (x z = 0) of the D3Q19 lattice in Figure 1, the unknown values of the distribution function f 5 , f 9 , f 13 , f 15 , and f 17 can be calculated as where N z are the transverse momentum corrections along the tangential directions that vanish at equilibrium but become nonzero when boundary-induced stresses are present.
To send out plane waves from the pressure boundaries, the local densities at the boundary nodes are locked to periodic functions that oscillate around f for finite time steps. A similar acoustic source was proposed in Viggen 39 and verified against the analytic solutions of viscously damped sinusoidal waves. For the pressure boundary nodes, which agitate sinusoidal plane waves, the local density varies according to where is the input angular frequency and t n the number of time steps during which (t) ≠ f . Because elastic wave propagation is a linear phenomenon and the Navier-Stokes equation can only be assumed to be linear for small disturbances, it is important to ensure ′ ≪ f so that nonlinear wave effects are avoided. Moreover, the computational Mach number M = |u max |∕c has to be much smaller than 1.0, because the macroscopic quantities given by the discrete-velocity Boltzmann equation converge to the solution of the incompressible Navier-Stokes equation with order M 2 .

VERIFICATION OF THE HYDRO-MICROMECHANICAL MODEL
The fluid-solid coupling scheme is benchmarked with single-particle sedimentation, while the acoustic source with a one-dimensional wave propagating in a pure fluid. The terminal velocities of a single sphere and the wave propagation/attenuation in a fluid are obtained from the simulations and then compared with the respective analytical solutions.
To verify the acoustic source for saturated granular systems, a granular chain of spherical particles is inserted into the same fluid domain. One-dimensional waves are agitated by the same acoustic source as in the second benchmark. Three additional benchmarks, including the propagation of an elastic wave in the same system in dry condition, are used to ascertain that the acoustic response due to the presence of the particles (and their motions) is reproduced correctly. In all simulations below, lattice units are used consistently for convenience, ie, the lattice spacing Δx = 1 and the time step Δt = 1 (for both LBM and DEM parts), which results in the lattice sound speed c s = 1∕ √ 3. The particles are monodisperse spheres with the radius R = R 0 + ||u n ||∕2, where R 0 = 5 and ||u n || = 0.01.

Terminal velocity of a single sphere
Although the LBM-DEM code used is well established and has been verified and validated with multiparticle systems, 28,43 a simple benchmark, ie, a single-particle settling in a fluid is considered here to showcase the robustness of the fluid-solid coupling scheme 42 as a function of the input parameters. Figure 3A shows the simulation setup to obtain the terminal velocity of a single-sphere settling in a viscous fluid. The sphere is released at the center of the cubic fluid domain (160 × 160 × 160). Periodic boundary conditions are adopted at all sides of the cubic domain. The sedimentation of the sphere is driven by a constant body force that accelerates and thus gives rise to a drag force on the sphere. The drag force F f − s , This simulation setup is well suited for the verification of the fluid-solid coupling, because the geometry is simple and an analytical solution for the evolution and the terminal velocity of a sphere is available at various Reynolds numbers. 21,46 From the particle Reynolds number Re = 2|V p −ū|R∕ with the kinetic viscosity andū the macroscopic flow velocity averaged over the whole lattice, the external body force F g balancing the drag force on the sphere can be computed analytically via where V t is the terminal velocity and is only related to the relaxation time once the Δx and Δt are set. In what follows, the influence of the relaxation time on the terminal velocity V t of a single sphere will be investigated by varying from 0.6 to 2.0, which corresponds to ∈ (0.033, 0.5) in lattice unit.
To ensure that the terminal velocity is reached, the simulations are run for sufficient numbers of time steps. The number of time steps needed increases dramatically with the decrease of the relaxation time and the corresponding increase of the particle Reynolds number Re. After running for sufficient time, the saturated values of |V p −ū| approximate the terminal velocities V t , which enters Equation 26 for the calculation of the drag force acting on the sphere. The relative deviation of the drag force depending on the choice of the relaxation time is shown in Figure 3B. The deviation given by the momentum exchange method is lower than 2% for ∈ (0.8, 1.8). The predicted drag force is seemingly independent of the relaxation time, especially in the intermediate Reynolds number regime. Interested readers should refer to Krüger et al 24

Wave propagation in fluid and dry/saturated granular chains
The propagation of a plane wave in a viscous fluid is chosen as the first benchmark for the verification of the acoustic source, following Viggen. 39,40 The 2D cross section at x 2 = 6 of the cuboidal fluid domain (12 × 12 × 620) is shown in Figure 4A. The acoustic source is located on the left-hand side (x 3 = 1), where a P-wave is agitated by changing the local density according to a sinusoidal waveform (see Equation 25). The pressure at the opposite side (x 3 = 620) is kept constant and equal to the mean of the sinusoidal waveform. The lateral boundaries are periodic to represent an infinite fluid between the two pressure boundaries. The boundary nodes have an off-equilibrium density magnitude of ′ = 1 × 10 −4 , with the mean density f = 1, the angular frequency = 0.03, and the number of the simulation steps t n = 5000.
The second benchmark considers the response of a dry, one-dimensional, granular chain, subjected to a mechanical impulse. The system is simulated using DEM only, as shown in Figure 4B. This benchmark does not involve the fluid phase, but it is preliminary to the following saturated systems. A total of 60 monodisperse spheres are inserted one after another in the x 3 direction at x 3 = i(2R − n ), i = 1, 2, … , 60, with x 1 and x 2 equal to 6. Confined by the leftmost and rightmost spheres fixed in space, the granular chain has an overlap of n = 0.01. To propagate a mechanical impulse in the dry granular chain, the first sphere is slightly shifted to the right at x 3 = 10.001, producing a small perturbation between the first and second spheres. The impulse is preferred over an oscillatory displacement/velocity, because the frequency domain response obtained with the impulse is clean and convenient for fitting, as shown in Figure 5F.
The third and fourth benchmarks consist in adding the granular chain into the representative fluid volume identical to the fluid domain of the first benchmark, and the same acoustic source and boundary conditions as in the pure fluid are applied, as shown in Figure 4C,D. In the third benchmark, all solid spheres except for the leftmost and rightmost ones are allowed to move and interact, whereas in the fourth benchmark, all spheres are fixed in space. With the particle motion completely constrained, the propagation of mechanical waves caused by interparticle collisions is excluded. By doing so, the influence of hydrodynamic interactions ( Figure 4D), accounted for by the MEM, and the interplay between mechanical and hydrodynamic interactions ( Figure 4C) on the acoustic behavior of the pore fluid can be studied separately. In the cases where the spheres are allowed to move and interact, the interparticle forces in normal direction are governed by Equation 11 with E p = 10 and p = 0.2. The same relaxation time = 1.0 is used in all benchmark cases. Note that the mechanical impulse is not used here, because an oscillatory displacement to particles in LBM-DEM simulations poses a significant numerical challenge. Because particles are discretized in the regular grid used by the fluid, in order to accurately resolve the small displacements over time, the resolution of the grid would have to be extremely fine, making the simulations too computationally expensive.
The propagation of waves at different speeds is primarily observed from the snapshots at t = 800, as shown in Figure 4 for all physical systems. Note that the spatial variations in the middle columns are truncated at x 3 = 450 approximately, because the waves have not yet reached further at t = 800. It can be seen that the acoustic attenuation in the x 3 direction differs depending on the presence and the mobility of the solid spheres (see Figure 4A,C,D). Figure 4D shows some specific features: the signal quickly dissipates, because the solid phase is completely rigid and the wave travels only through the pore fluid; however, the oscillatory nature does not vanish entirely, as also shown in 5K. Because the total travel distance becomes longer, the saturated medium as a whole appears to be more viscous than the pure fluid. To quantify the propagation speed and attenuation in the saturated systems, the off-equilibrium density ′ within each cross section of the fluid domain perpendicular to the x 3 axis is averaged at each time step. The resulting evolution of the averaged ′ in space and time is plotted for the pure fluid in Figure 5A, the saturated chain of mobile particles in Figure 5G, and the saturated chain of fixed particles in Figure 5J. For the dry granular chain, the elastic wave propagation can be inferred from the evolution of the particle velocity in space and time, as shown in Figure 5D.
The propagation of the averaged ′ in the x 3 direction at t = 800 is extracted from Figure 5A,D,G,J. Those are plotted in Figure 5B,E,H,K and compared with the analytical solution derived as follows. Thanks to the acoustic sources that emit the sinusoidal waves, the acoustic responses in the first, third, and fourth benchmark cases can be described analytically by the stationary solution of the one-dimensional lossy wave, namely, where A is a constant, k is the angular wavenumber, and s is the relaxation time for acoustic attenuation s = 2 ∕c 2 s ; s is the spatial absorption coefficient, which is negligible in low-frequency ranges but becomes more pronounced as the frequency increases. Note that s is related to the relaxation time of the lattice BGK model (4) 27 is A, which is determined from the variation of cross-section averaged ′ from the source at x 3 = 1 to x 3 = 450, as shown in Figure 5B. The solid curve shows that the behavior predicted by the LBM simulation correctly matches the analytical solution with ∠A Re ∕ ′ = −0.023 and A Im ∕ ′ = −1.003.
With the acoustic source verified in the pure fluid, we now let c s and s be the additional free parameters, in order to quantify the propagation speed and attenuation in the saturated granular chains. By fitting Equation 27 with the acoustic responses, we obtain A Re ∕ ′ = −0.578, A Im ∕ ′ = −0.117, s = 6.627 × 10 −3 , and c s = 0.698 for the saturated granular chain of mobile particles, and A Re ∕ ′ = −1.22, A Im ∕ ′ = −0.236, s = 4.819 × 10 −2 , and c s = 0.391 for the saturated granular chain of fixed particles. Figure 5H,K shows the comparison between LBM simulations and analytical solutions, with a good agreement in both cases. It is noteworthy that when the constituent particles are allowed to interact elastically, the wave undergoes much less attenuation and the sound speed almost doubles.
The discrete Fourier transform (DFT) of the space-time domain data gives the dispersion relations of the single/two-phase systems, as shown in Figure 5C,F,I,L. The sound speeds can be inferred from a line connecting the inserted frequency and the most significantly agitated wavenumbers (white broken lines), eg, the sound speed in Figure 5C that is c s ≈ 1∕

√
3, which proves that the acoustic wave propagates with accurate speed in the pure fluid. Note that the slopes agree well with the sound speeds fitted with Equation 27. Interestingly, the dispersion relation of the dry granular chain in Figure 5F is highly nonlinear because of the interparticle collisions that give rise to frequency dependence. It appears that the nonlinear dispersion relation of the dry granular chain is completely suppressed by the coupled motion of the pore fluid and the hydrodynamic interactions, as shown in Figure 5I. The sound speed therein is higher than both in the pure fluid and the dry granular chain. We will discuss this point further in Section 5.
Finally, by constraining the degrees of freedom of all solid particles (see Figure 5J,K,L) while allowing a flow in the pore space, the sound speed significantly decreases, as can be seen by comparing Figure 5G and 5J. This is because the fluid path along which the acoustic wave propagates is always longer than the travel distance. In addition, when the solid phase is rigid, the only attenuation mechanism is the pore-scale fluid flow. When the solid particles are allowed for elastic interactions, the attenuation becomes less pronounced because of the in-phase motion of the solid particles and the pore fluid. Nevertheless, the acoustic attenuation in the pure fluid is much weaker than in the dry/saturated granular chains, as indicated by the spatial absorption coefficients and shown in Figure 5D,G,J.

MODELING ELASTIC WAVE PROPAGATION IN SATURATED GRANULAR CRYSTALS
As envisaged by Biot, 32 the simplest setup to evaluate the so-called viscodynamic operator is an FCC packing of equally sized spheres, pushed by an alternating motion from the fluid. The ordered FCC packing is apparently more complicated than the particle chains in Section 3.2. However, the elastic wave speeds of such packings can be computed analytically in dry conditions by means of the effective medium theory. 57 Such configuration is then selected for further analysis, including the influence of the acoustic sources and effective confining pressure on the wave propagation. In order to simulate the propagation of waves in saturated granular media, the discretization parameters are carefully selected so as to meet the true sound speed in water.

FCC packing of frictionless spheres
The FCC packing of spheres is shown in Figure 6. In a similar fashion as the benchmarks in Section 3.2, the oscillating pressure boundary condition is applied on the left-hand side (x 3 = 1Δx) of the fluid domain (40Δx × 40Δx × 1280Δx), while a constant pressure equal to the mean of the oscillatory pressure is maintained on the right-hand side. The other boundary conditions are assigned to be periodic in order to mimic an infinitely large, homogeneous and fully saturated packing. Within the cuboidal fluid domain, a 4 × 4 × 200 FCC packing of equally sized frictionless spheres is inserted. The leftmost and rightmost layers of solid spheres are fixed in space, so that effective stresses on the FCC packing are kept constant. The radius of each sphere is set to R = R 0 + 0.5|u n |; R 0 = 5Δx is the radius that leaves zero overlap between the spheres in the FCC packing and |u n | is the overlap that gives rise to the effective stress.

Model parameters and numerical aspects
Setting the model parameters of LBM-DEM simulations can be difficult, especially in the case of wave propagation. Relevant physical parameters need to be properly selected to reproduce both the flow and acoustic behaviors. The kinematic viscosity depends on the choice of the spatial and temporal resolutions and a single model parameter . For the computational Mach number M to be sufficiently small, the lattice sound speed c s is tuned to be much larger than the maximum flow velocity |u max | via x and t. For LBM simulations of fluid flow, M is typically smaller than 0.1, meaning that c s is not necessarily the true sound speed of the fluid. One can choose a preferred spatial resolution x and then for the target viscosity, which in turn gives the corresponding value for t. However, this cannot be the case for LBM-DEM modeling of acoustic waves. First, c s has to be chosen equal to the true sound speed in the fluid (eg, c s = 1500 m/s for water), leaving dependent on and x, that is, = c 2 s − √ 3c s Δx∕6. Second, for the LBM-DEM simulations to be numerically stable, the BGK collision operator requires to be larger than 0.5. As a consequence, to match the kinematic viscosity of water, = 10 −6 m 2 /s, with the smallest allowable relaxation time = 0.5, the spatial resolution x would have to be 2.598 km, which is obviously beyond the scale of micromechanics. If one wants to use the length scale of a typical solid particle, for example, x = 10 −3 m, the resultant then is as small as 1.92 × 10 −7 . Therefore, in this work, we choose to only match the sound speed of water rather than the viscosity. Ongoing work involves the implementation of a regularized LBM to reduce the lower limit of the relaxation time (not shown here).
For the LBM-DEM simulations of wave propagation in the saturated FCC crystal, we use the parameters listed in Table 1 in both dimensionless and dimensional units. The LBM and DEM calculation cycles share the same time step t. In order to investigate the effect of the acoustic source, a variety of input waveforms and frequencies are considered (see Table 2). In contrast to the continuous signal in Section 3.2, the saturated FCC systems in this section are agitated by a single-period pulse. Different values for the overlap |u n | between solid spheres are selected, such that the effective confining pressure on the packing varies from 0.1 to 30 MPa. Figure 7 shows snapshots of a typical LBM-DEM simulation of wave propagation in an FCC packing of spheres. The pressure wave is agitated by a cosine signal with the magnitude of the off-equilibrium density ′ = 0.33 kg/m 3 at an input frequency of 1.30 MHz. Similar to Section 3.2, the macroscopic fluid velocityū = ∑ u∕ ∑ is averaged within each cross section (within the solid spheres, both |u| and are zero). Applying the DFT to the evolution ofū 3 in time and the propagation direction x 3 gives the P-wave dispersion relation of the pore fluid. Similarly, the dispersion relation of the solid phase can be obtained from the particle velocity vectors V p . In particular, the components of V p along x 3 and perpendicular to x 3 reflects the P-and S-wave propagations, respectively.

Comparison of wave propagation in dry and saturated FCC packings 4.3.1 Mechanically agitated impulse
The acoustic source at the pressure boundary aims to simulate the propagation of acoustic waves from a viscous fluid to a saturated granular medium. Alternatively, the wave can be agitated by perturbing the solid phase at a given position, 11,58,59 with an interparticle force slightly bigger than elsewhere, eg, between the boundary particles fixed in space and their neighbors. The resulting unbalanced forces on the neighboring spheres, in turn, induce a mechanical impulse propagating into the granular packing. For FCC packings, if the unbalanced force is aligned with the propagation direction, a P-wave will be agitated; otherwise, an S-wave will be triggered. The mechanically agitated impulse is well suited for DEM modeling of wave propagation in dry granular materials. In this section, we aim to compare the propagation of waves in dry and saturated granular systems. Therefore, the mechanically agitated source is chosen in order to keep the same type of acoustic source between the dry and saturated cases.

Time-domain and frequency-domain responses
To understand the effect of pore fluids on the dispersion relations of granular crystals, the acoustic responses of the dry and saturated FCC packings under an effective confining pressure of 5 MPa are compared. The elastic waves in both cases are agitated by a perturbation applied in the leftmost layer of solid spheres (x 3 = 2Δx), pointing towards the propagation direction x 3 . Following the same averaging scheme as in Section 4.2, the time-and frequency-domain responses of the saturated and dry FCC packings are obtained and plotted in Figures 8 and 9. It can be observed from the space-time signal in Figure 8A that the main pulse broadens as the wave travels through the saturated FCC packing, before reflection. Interestingly, there exists another wave, which is located almost close to the source (x 3 = 2Δx). This wave is slow and highly dissipative and has features typical of the slow P-wave predicted by Biot. 32 Figure 8B shows the amplitude spectrum of the time-domain signal of the averaged fluid velocity. The dispersion relation, relating the frequency and the wavenumber k, is almost linear. Figure 9A shows that the P-wave dispersion relation of the dry FCC packing is highly nonlinear, ie, the wave velocity decreases with the increase of the wavenumber and the frequency. 13,14 The group velocity approaches zero in the very high-frequency regime, as shown in Figure 9A. The amplitude spectrum of the particle velocity components along x 3 in Figure 9C highlights that the P-wave dispersion relation of the solid phase is identical to the one for the fluid phase of the saturated FCC packing (see Figure 8B). This proves that the in-phase motions of the pore fluid and the submerged solid spheres indeed cause the propagation of P-waves in saturated granular media. By comparing Figure 9A and 9C, one can find that the fluid-solid coupling qualitatively changes the P-wave dispersion relation of the solid phase. It is known that the P-and S-waves are decoupled for regular crystalline structures like the FCC configuration, ie, no shear wave is induced by P-wave propagation and vice versa. Therefore, it is not surprising that no S-wave branches are present in Figure 9B. In the case of the saturated packing, however, the S-wave branch appears ( Figures 9D) and is nearly linear in the low-frequency regime. Nevertheless, the S-wave is significantly more dissipative and dispersive than the P-wave, as shown by the small Fourier amplitude in Figure 9D. A few inclined branches that have the same slope as the P-wave dispersion relation also enter the frequency domain in Figure 9D. These weak branches are caused by fluctuations due to frequently updated fluid-solid links (see Figure 2) due to strong longitudinal motions of the solid particles.

Effect of acoustic source
Unlike typical ultrasonic experiments in which wave velocities are inferred from the time evolution of signals, numerical simulations allow for direct measurements of wave velocities from dispersion branches in the frequency domain. 13,14 For dry granular materials, it is known that input frequencies affect the identification of travel time and travel distance from received signals and, in turn, wave velocities. Therefore, the wave velocities derived from dispersion branches, obtained from simulations, can be treated as the "ground truth." In this section, we investigate the effect of input waveforms and frequencies on the wave velocities in saturated granular materials.

Effect of input waveform
Three waveforms, namely, single-period step, cosine, and sine functions, are tested in order to study the effect of input waveforms on the acoustic response of saturated granular media, as described in Table 2 and Figure 10.  Figure 11A,D,G shows the evolution of the cross-section averaged fluid velocityū 3 in time and space, in response to the step, cosine, and sine signals. In the case of the square signal, the local fluid velocity shows a sharp increase near the source that stays almost constant in the remaining time (see Figure 11A). Differently, the cosine signal guarantees a smooth transition of local densities at the source from equilibrium to nonequilibrium and then back to equilibrium during t ≤ t n Δt. The sine signal, however, inserts a peak and a trough of equal magnitudes around the equilibrium density, which leads to more significant dissipation while the signal is propagating, as shown in Figure 11G. Similar to Figure 9C,D, Figure 11B,E,H and 11C,F,I shows the P-and S-wave dispersion branches for the respective input waveforms. Compared with the sine signal, both the square and the cosine signals activate the P-wave dispersion relations in the low-wavenumber range (see Figure 11B,E). The P-wave dispersion relation in Figure 11H is more pronounced in the high-frequency range ( ∕2 > 2 MHz), which is associated with less broadened signal in Figure 11G.
Regardless of the difference in the amplitude spectrum, the slopes of all dispersion branches are seemingly identical. Interestingly, a few frequency bands seem to be more or less active for the different inserted signals. While frequency bands caused by the sine and cosine signals are consistent, the square signal creates low-intensity frequency bands with equal intervals but activates higher frequencies as well.
Similar features can be observed in the S-wave dispersion relations in Figure 11C-I. Additionally, the transverse motions of the solid particles induced by the square and cosine signals of the pressure wave are stronger than those induced by the sine signal. The intensity on the S-wave dispersion branch in Figure 11F is comparable with the noise induced by the P-wave, making it difficult to identify the slopes in Figure 11I accurately. Overall, the P-and S-wave dispersion branches in Figure 11E,F are less fragmented by low-intensity frequency bands, which leads to more accurate estimates of the wave velocities. Therefore, we use the cosine waveform in the following sections, although the dispersion relations in the low-frequency ranges are clear for both the square and cosine signals.

Effect of input frequency
Using the cosine waveform and the input frequencies given in Table 2, the P-wave and resulting S-wave dispersion relations are obtained for three different values of the input frequency. The space-time evolution of the longitudinal particle velocity is plotted in Figure 12A-G, and the P-wave and S-wave dispersion branches in Figure 12B-H and 12C-I, respectively. Pulse broadening can be observed in Figure 12A, ie, the high-frequency signals decay rapidly as the wave propagates in time and space. With decreasing input frequency, the pulse broadening phenomenon becomes less pronounced as shown in Figure 12D,G, which reflects the dispersive nature of the waves in granular media. As the input frequency decreases, the dispersion relation becomes clearer at the frequencies smaller than approximately twice as large as the input value (1.30, 3.25, and 13.0 MHz in Figure 12B, 12E, and 12H, respectively). In contrast to the dispersion relation of dry granular media in Figure 9A, the largest Fourier coefficients always appear in the low-frequency ranges, which are associated with the frequency bands activated by the input frequency. Figure 12C,F shows very noisy patterns, with additional branches caused by the longitudinal motion of the particles. Using low input frequencies, eg, f = 1.30 MHz, the noisy signals are reduced, as shown in Figure 12I. Note that the S-wave dispersion relations here are obtained with the transverse velocity components of the solid particles. Although not shown here, only the P-wave dispersion branch is obtained from the DFT of the cross-section averaged fluid velocity in the transverse directionū 3 . This is reasonable because the pore fluid cannot sustain shear.
Irrespective of the input frequencies, clear dispersion relations can be identified, which are consistently linear for both the P-and S-waves. In order to accurately fit each dispersion relation with a straight line, a large number of relevant data points in the frequency domain are needed. Therefore, we select an input frequency of 13.0 MHz for the following LBM-DEM simulations at various effective confining pressures, while keeping a sufficient number of time steps during which the acoustic source remains active for a single period. Biot 32 derived theoretical equations for the wave velocities in a fully saturated poroelastic medium, given the elastic constants of the dry solid matrix and the characteristics of the fluid. Biot theory incorporates the governing mechanisms such as the viscous and inertial interaction between the pore fluid and solid matrix, neglecting the pore-scale interaction between the constituent particles and local flow. The present hydro-micromechanical model is capable of numerically reproducing the relevant mechanisms in the poroelastic medium, ie, granular interactions, fluid-solid coupling, and pore-scale fluid flow in detail. Therefore, the wave velocities as predicted by Biot theory are well suited for the validation of the present numerical model at long wavelengths (much larger than pore/particle size). The goal here is twofold. First, the numerical and analytical predictions are compared. Second, results from the direct simulations are analyzed to investigate the characteristic features as predicted by Biot theory.

WAVE VELOCITIES FROM THE HYDRO-MICROMECHANICAL MODEL AND BIOT THEORY
The system chosen is the same FCC packing of monodisperse spheres as in Section 4. In order to highlight the dependence of wave velocities on pressure, the overlap and thus interparticle forces between the solid spheres are controlled, such that the macroscopic effective confining pressure ranges from 0.1 to 30 MPa. The same model parameters as in Table 1 are adopted here. Given the discussions in Section 4.4, a single-period cosine signal at an input frequency of 13.0 MHz is consistently used at various levels of effective stress.

Dispersion relations at various effective confining pressures
In Figure 13, we show the variation of the P-and S-wave dispersion relations with increasing effective stress. As expected, the slopes of the tangent lines, ie, the wave velocities, increase with the effective confining pressure, for the P-waves in both the solid and fluid phases (Figure 13A,D,G and Figure 13B,E,H) and for the S-waves in the solid phase ( Figure 13C,F,I). Interestingly, the P-wave dispersion relations in the solid and in the fluid are identical for each level of effective confining pressure, as indicated by the slopes (green lines) in Figure 13. This confirms again that pressure waves travel in the pore fluid and in the solid skeleton with the same velocity. Wide horizontal bands appear in the frequency domain response of the fluid phase ( Figure 13A,D,G), approximately around the input frequency (f = 13.0 MHz), and merge with the P-wave branches. However, the horizontal bands are not present in the spectra of the solid phase ( Figure 13B,E,H). The bands are possibly associated with the input frequency inserted from the acoustic source in the fluid domain. Figure 13C,F,I shows the S-wave dispersion branches induced by the P-waves at three effective pressure levels. The pressure dependence of the S-wave velocities is more pronounced than for the P-waves. Because of the presence of interstitial fluid in the pores, the pressure pulse induces also a shear wave in the granular crystal. However, the effect of the pore fluid on the S-wave velocity is almost negligible. Therefore, the behavior of the S-wave velocity with increasing effective stress in saturated granular crystals strongly resembles that in dry condition, irrespective of the acoustic source in the fluid or solid phase, as predicted by Biot theory.

Biot equations of wave velocities in saturated poroelastic media
Biot theory predicts the wave velocities in saturated poroelastic materials, depending on the pore geometry, densities and elastic properties of the solid skeleton, constituent grains, and pore fluid. The material properties are summarized in Table 1. The pressure-dependent bulk and shear moduli of the dry FCC packing K dry and G dry can be calculated following the effective medium theory 57 when the Young modulus E p and Poisson ratio p of the solid particles (see Table 1) and the particle configuration are known. We refer the reader to Pagano et al 60 for detailed expressions of K dry and G dry .
Knowing K dry and G dry , the porosity , the bulk modulus of the solid particles K p from E p and p , the bulk modulus of the pore fluid K f from the sound speed in the fluid c s , and the solid and fluid densities p and f , the wave velocities in the saturated FCC packing are given by 11 and 22 are the effective mass of the solid matrix and the fluid, incorporating the induced mass 12 caused by inertia drag that arises from the relative acceleration between the solid particles and the pore fluid,̂is the total apparent mass of the saturated FCC packing with the tortuosity parameter = 1 − r(1 − 1∕ ) (r = 0.5 for spheres), P, Q, and R are the elastic coefficients defined as The two solutions of v p given in Equation 28 correspond to the fast and slow P-waves. The fast P-wave can be easily measured in both the laboratory and the field, and it reflects the in-phase motion between the solid matrix and the pore fluid. As shown in Equation 28, Biot theory predicts a secondary slow P-wave, which is highly dissipative resulting from the overall out-of-phase motion between the solid and the fluid. The S-wave velocity v s given in Equation 29 incorporates an inertia coupling term which modifies G dry with the apparent mass and the geometry of the packing.
From the P-and S-wave dispersion relations in Figure 13 (see the Supporting Information for the amplitude spectra at other effective confining pressures), the wave velocities from the hydro-micromechanical model are obtained. The numerical predictions and the analytical solutions are compared in Figure 14. As predicted by Biot theory, the in-phase motion of the pore fluid flow and the FCC packing of solid spheres leads to higher P-wave velocities compared with those of the dry packing. On the other hand, the inertia drag between the particles and the surrounding pore fluid increases the effective mass of the solid phase and thereby slightly reduces the S-wave velocities. Both phenomena are qualitatively reproduced including the pressure dependence, as shown in Figure 14A,B. As the effective confining pressure increases from 0.1 to 30 MPa, the relative deviation between the analytical solution and the LBM-DEM prediction of the P-wave velocity decreases from 19.5% to 10.5%, whereas the relative deviation for the S-wave velocity drops from 3.5% to 0.15%. It is not surprising that the analytical solutions of Biot do not match perfectly with the simulation results because some mechanisms such as viscous losses within the pore fluid and pore-scale squirt flow are not taken into account. In addition, the FCC packings of solid spheres are weakly anisotropic rather than isotropic because the left and right boundaries are not periodic (see Section 4.1).

Numerical evidence of the slow compressional wave
With the help of the coupled LBM-DEM, not only the pressure-dependent dispersion relations can be obtained, but mechanisms like slow-wave diffusion can be investigated as well. The weakly inclined dispersion branches in Figures 9C, 11B-E, and 12B-H suggest the existence of the slow, dissipative P-wave that propagate slowly from the acoustic source, irrespective of the difference in the acoustic source (eg, perturbation in fluid/solid, waveforms, and frequencies). The same feature is highlighted in the insets of Figures 8A, 11A-D, and 12A-G. The fact that the slow P-waves exist irrespective of the acoustic sources on the solid or fluid phase (see Figures 9C, 8A, and 11D,E) also suggests that the slow waves are not numerical artifacts.
To better illustrate the slow waves, the propagation of the locally averaged fluid momentum fluctuationū 3 along the x 3 direction at various times are plotted for the square, cosine, and sine waveforms in Figure 15. The lags between theū 3 signals at various times can be utilized to estimate the P-wave velocities. All subplots of Figure 15 show the propagation of the fast P-waves as indicated by the wavefronts on the right-hand side. From Figure 15A-C, the peaks relevant to the slow P-waves can be observed near the source, most clearly from the ones given by the cosine input signals. Although the signals agitated with the square waveform are noisy and contain high-frequency contributions (coda), as shown in Figure 15A, the peaks and troughs therein are seemingly correlated in time, which suggests the presence of the slow P-wave. With the sinusoidal waveform, the attenuation is so significant that the signals near the source are hardly noticeable. From Figure 15A-C, it is confirmed again that the cosine waveform is best suited to agitate pressure waves from the fluid phase.
The effect of input frequencies on the slow P-waves is investigated in Figure 15D-F. With decreasing input frequency, the existence of the slow P-wave becomes more evident. However, the travel distances and the wave velocities of the slow P-waves at different times cannot be measured because of the slowness and high diffusivity of the waves. The propagation ofū 3 at different times along the x 3 direction associated with the influence of effective confining pressure is plotted in Figure 15G-I, respectively. In addition to the increasing fast P-wave velocity, the slow P-wave becomes more dissipative as the effective confining pressure increases, which can be seen from the increasingly vanishing peaks relevant to the slow P-waves at t = 10 −3.8 ms. The enhanced dissipation may be attributed to the increasing natural frequencies of the FCC packing and the slightly enlarged overlap needed for the elevated effective confining pressure.
To accurately calculate the slow P-wave velocities, two numerical issues need to be tackled: (1) one has to overcome the lower limit of the relaxation time to allow for realistic fluid viscosities in LBM simulations of wave propagation; (2) the size of the numerical sample, at least in the propagation direction, has to be large enough so that the slow P-wave can propagate sufficiently far, before the reflection of the fast P-wave interferes with the slow one. Ongoing work aims to overcome the first issue with a regularized version of the LBM.

CONCLUSIONS
A hydro-micromechanical model has been proposed to simulate wave propagation in saturated poroelastic granular crystals. The hydrodynamics in the pore fluid is resolved with the lattice Boltzmann method, and the translational and rotational motion of solid particles with the discrete element method. The novelty of the work lies in the application of the two-way coupled LBM-DEM method to study the fluid-solid interaction at local level, induced by the small amplitude wave traveling in a saturated medium. As a new implementation, the oscillating pressure boundary can emit the correct acoustic waves from the fluid phase into the fluid-solid mixture, as inspired by Biot original model. The fluid-solid coupling scheme is benchmarked against the terminal velocities of a single sphere settling in viscous fluids at various Reynolds numbers. The relative deviations are less than 2% for intermediate Reynolds numbers, and the accuracy is found to be independent of the choice of the relaxation time. The oscillating pressure boundary is applied to a pure fluid and two saturated granular chains either with the particles fixed or mobile. The acoustic responses are compared with one-dimensional lossy wave equations to infer the spatial absorption coefficients and the wave velocities, which agree excellently with those derived from the dispersion relations.
Furthermore, the FCC granular crystals of equally sized frictionless spheres are inserted into the fluid domain for studying acoustics in saturated granular materials. By activating/deactivating the hydrodynamics, it is proved that the presence of the pore fluid changes the dispersion relations from nonlinear to almost linear. The compressional-and shear-wave dispersion relations, which are uncoupled in dry FCC configurations, become coupled under the influence of the hydrodynamic interactions with the solid particles. Among various input waveforms, the cosine is the most smooth and leads to the cleanest signals in both the time and frequency domains.
Finally, using the cosine as the acoustic source at the boundary, the compressional-and shear-wave dispersion relations are obtained under a wide range of effective confining pressures. The compressional-and shear-wave velocities qualitatively agree with the analytical solutions given by Biot theory, with relative deviations between 19.5% and 10.5% for the former and 3.5% to 0.15% for the latter, at an effective pressure ranging from 0.1 to 30 MPa. In addition to the fast compressional waves, the slow waves as predicted by Biot theory are numerically reproduced in the LBM-DEM simulations. Again, the cosine input waveform yields the cleanest signals that show clear correlations in time that indicate the existence of the slow waves.
In conclusion, our enhanced hydro-micromechanical model is perfectly suited to study wave propagation in poroelastic granular media. It allows exploring physical coupling mechanisms neglected in the classical theory of poroelasticity and the influence of frequencies ranging from medium to very high. As next step, detailed investigations of pore-scale squirt flow will be performed for which modified Biot equations are needed to take the squirt flow into account. Moreover, the numerical tool will be extended to the simulation of random, frictional packings in dry and saturated configurations. Finally, the coupled LBM-DEM code has the potential to simulate a wide range of saturation degrees. Although currently, it is only possible to simulate mixtures of fluids with significantly different densities and viscosities, yet still far from the ratio between common fluid and gases, the LBM-DEM simulations could still provide a better understanding of acoustics in partially saturated granular media. On the methodological side, ongoing work aims at overcoming the lower limit of the relaxation time, which will make it possible to simulate wave propagation in fluids with realistic viscosities.