A consistent derivation of soil stiffness from elastic wave speeds Open Access

A

stress constant for M_w

B

stress constant for G_w

b

stress dependency exponent

$Cαβγζ$

homogeneous elastic stiffness tensor

$d¯$

constant for solid fraction dependence in hyperelastic model

e

void ratio

e_min

minimum void ratio

e_max

maximum void ratio

$F(ϕ)$

solid fraction function for hyperelastic model

$G¯$

shear stiffness of hyperelastic model

$G¯hypo$

shear stiffness constant for hypoelastic model

G_s

specific gravity

G_t

tangent shear modulus with respect to $εs$

G_w

shear moduli

$H(ϕ)$

general function of $ϕ$ for hypoelastic model

i

unit imaginary number

$K¯$

bulk stiffness of hyperelastic model

$K¯hypo$

bulk stiffness constant for hypoelastic model

k

wave number

$M¯$

constrained modulus for hyperelastic model

M_t

tangent constrained modulus

M_w

constrained moduli

n

pressure dependence constant for hyperelastic model

p

effective (elastic) pressure

p_atm

atmospheric pressure

q

elastic deviatoric shear stress invariant

$U$

internal energy density

$Uα$

magnitude of Fourier mode

u

displacement

$u¨$

acceleration

V_p

longitudinal (p-) wave velocity

$Vphyper$

longitudinal wave velocity for hyperelastic model

$Vphypo$

longitudinal wave velocity for hypoelastic model

V_s

shear (s-) wave velocity

$Vshyper$

shear wave velocity for hyperelastic model

$Vshypo$

shear wave velocity for hypoelastic model

$δαβ$

Kronecker delta

$εαβ$

strain tensor

$εαβT$

total strain tensor

$εv$

elastic volumetric strain

$εs$

elastic deviatoric shear strain invariant

$εaT$

total axial strain

η

ratio of q and p

ρ

bulk density

$ρs*$

unstressed solid density

$σαβ$

stress tensor

σ_r

radial stress

$ϕ$

solid fraction

$ϕ0$

initial solid fraction

ω

circular frequency

$□*$

quantity at reference state

$□˜$

perturbed quantity

$□dev$

deviatoric tensor

Accurate determination of wave speeds and small-strain stiffness is crucial for a thorough geotechnical evaluation of soil sediments (Atkinson, 2000; Mayne, 2001) and the precise design of geotechnical structures (Clayton, 2011). In geomechanics, the wave speeds are conventionally assumed to be equivalent to the ones derived for linear elastic materials (Hardin & Richart Jr, 1963; Hardin & Blandford, 1989; Wichtmann & Triantafyllidis, 2009; Li et al., 2022). Based on this assumption, empirical relations are obtained to account for the observation that instantaneous stiffness is a non-linear function of the pressure, density, and other state variables. However, assuming the wave speeds are defined for a linear elastic solid is generally an oversimplification for porous media and inconsistent with empirical relationships, necessitating the exploration of more general non-linear elastic models. Here, wave speeds are derived for non-energy-conserving hypoelastic as well as energy-conserving hyperelastic models, assuming a general stiffness dependence on pressure and density while ensuring mathematical consistency.

Conventionally, in geotechnical engineering, wave speeds are assumed to be characterised by a homogeneous, isotropic, linearly elastic solid. Under these assumptions the longitudinal and transverse wave speeds are defined as

where M_w and G_w are the constrained and shear moduli, respectively. After experimentally measuring the wave speeds, the stiffness moduli M_w and G_w are found to depend on the effective pressure p and solid fraction $ϕ$ in soils. Generally, these empirical relations adopt the form (Hardin & Richart Jr, 1963):

where A and B are constants; p_atm = 1 Pa is a reference stress; b is the stress exponent reflecting non-linear pressure dependency; and $H(ϕ)$ is a general function of solid fraction $ϕ$⁠. Note that the void ratio $e=(1−ϕ)/ϕ$ is normally adopted in geotechnical engineering rather than the solid fraction $ϕ$⁠, but the current authors opt for a description in terms of solid fraction for convenience. Typically, $1/3 ≲ b ≲ 1/2$ for particles spanning from perfectly spherical to angular (Santamarina & Cascante, 1998), as validated through effective medium representation (Goddard, 1990; Yimsiri & Soga, 2000; Agnolin & Roux, 2007). Furthermore, more advanced relations have been proposed to include general stress states (Roesler, 1979; Yu & Richart Jr, 1984; Hardin & Blandford, 1989; Bellotti et al., 1996) and other state variables such as grain-size distribution (Wichtmann & Triantafyllidis, 2009, 2010; Payan et al., 2017), particle shape (Liu & Yang, 2018; Tang & Yang, 2021), or fabric (Mital et al., 2020; Li et al., 2022). However, within the current paper, these options are not explored, as the authors seek to solely understand the effect of pressure and density dependence on the wave speed.

Although these empirical relations have successfully captured experimental phenomena, it is crucial to note that they were determined by first assuming a linearly elastic isotropic continuum. This introduces a discrepancy between the assumed elastic wave speeds, as given by equations (1) and (2) and the retrieved empirical equations (3) and (4). Fig. 1 highlights this inconsistency, where soils are known to have state-dependent stiffness, yet state-independent linear elasticity is assumed to derive the wave speeds. Despite existing models that incorporate pressure dependence in their constitutive laws (Jiang & Liu, 2003; Einav & Puzrin, 2004; Houlsby et al., 2005; Nguyen & Einav, 2009), wave speeds are usually not derived with consideration for pressure dependence. However, hyperelastic pressure-dependent models assuming specific pressure dependencies have recovered empirical relations under isotropic compression states (Andreotti et al., 2013) and uniaxial compression states that showed shearing alter the evolution of the longitudinal-to-transverse wave speed ratio (Mayer & Liu, 2010).

Herein, analytical wave speeds are derived for both hyperelastic and hypoelastic constitutive models with a general pressure and density dependence. Importantly, the passage of energy is conserved in the hyperelastic model, but this is not the case for the hypoelastic model. Furthermore, this derivation procedure remains accurate and general, irrespective of the model selected. The derived wave speeds reproduce the empirical relations under isotropic stress conditions regardless of whether hyperelasticity or hypoelasticity is considered. However, due to the density dependence, a minor, typically negligible, correctional term arises for the hyperelastic model. Conversely, when the material undergoes shear, the wave speed derived for hyperelasticity is significantly different from the hypoelastic wave speed because of the presence of additional non-negligible stress-dependent terms. These insights shed light on the evolution of the ratio between longitudinal and transverse wave speeds observed experimentally that previously could not be explained with empirical relations when assuming a linearly elastic isotropic continuum.

The analytical expressions for longitudinal and transverse wave speeds will now be derived. The momentum balance with the assumption of negligible body forces is given by:

where ρ is the bulk density; $u¨α$ represents acceleration in the α direction; $σαβ$ is the total stress tensor; and □_,α$=∂□∂xα$⁠. Note that by neglecting body forces it is assumed that the domain of interest is small, as in conventional geotechnical experiments. However, this may not always be the case, and a more robust derivation is required for such more complex scenarios. For a homogeneously deformed solid at equilibrium, this simplifies to:

where $σαβ*$ is the homogeneous reference (or unperturbed) stress. By first-order Taylor series expansion, the stress can be expressed in terms of $σαβ*$ and a stress perturbation $σ˜αβ$⁠, and is given by:

where $Cαβγζ*$ is the homogeneous elastic stiffness tensor at the reference state; $u˜α$ is a displacement perturbation from the reference homogeneous solution of the displacement field $uα*$⁠, such that $uα=uα*+u˜α$ (Stefanou & Alevizos, 2016; Stefanou & Gerolymatou, 2019). The total strain tensor is defined as $εγζT=12(uγ,ζ+uζ,γ)=uγ,ζ$ as small deformations have been assumed, and that there are no rotational contributions, and thus, $ε˜γζT=u˜γ,ζ$ is the perturbed strain tensor. By inserting equation (7) into equation (5) and using a first-order Taylor expansion of the left-hand side, assuming no acceleration in the reference state (⁠$u¨α*=0$⁠) the momentum balance is expressed as:

where $ρ*$ is the homogeneous equilibrium density. The perturbed stress tensor $σ˜αβ$ will be determined later from the constitutive law through first-order Taylor series expansion.

First, consider longitudinal waves that propagate along the z-axis. The perturbed displacement $u˜z$ is decomposed in Fourier modes of the form $u˜z=Uzei(kz−ωt)$⁠, where U_z is the amplitude, i is the unit imaginary number, k is the wave number and ω is the circular frequency. The perturbed strain tensor and its deviatoric component are thus expressed as:

where $ε˜αβdev=ε˜αβ−13ε˜γγδαβ$ is the deviatoric component of the perturbed strain tensor with $δαβ$ being the Kronecker delta.

For a transverse wave, it is assumed that the wave propagates in the z-axis and is polarised in the x-direction, where the perturbed displacement $u˜x=Uxei(kz−ωt)$⁠. The disturbance to the strain field and its deviatoric component for the transverse wave are equal in this case and given by:

Figure 2(a) illustrates these wave perturbations for a general state, where the periodic blue arrow shows the longitudinal wave and the red displays the transverse wave. The periodic gradient shows that the particle motion is in the same direction as the wave propagation for the longitudinal wave. In contrast, the particle motion occurs perpendicular to the wave direction in the transverse wave case.

To determine the wave speed, a constitutive model is required to define the homogeneous elastic stiffness tensor, as this enables an explicit representation of equation (8) following the substitution of the perturbed strains into the perturbed stress. The following section describes two potential constitutive relations to apply these strain perturbations.

Now that a general framework for wave speed derivations has been laid, a constitutive law is required to evaluate the perturbed stress. First, a hyperelastic model for pressure and density-dependent media is introduced. The advantage of hyperelasticity is that it ensures the energy is being conserved as waves travel through the medium, thus keeping faithful to the meaning of a true elastic body. Here, an internal energy density $U$ is considered, that broadens previous work by Alaei et al. (2021), and depends on solid fraction $ϕ$ and elastic strain tensor $εαβ$⁠:

where $K¯$ is the bulk stiffness constant; $G¯$ is the shear stiffness constants, respectively; $εv=εαα$ is the volumetric elastic strain; and $n≥0$ controls the pressure dependence of the stiffness tensor $Cαβγζ=∂2U∂εαβ∂εγζ$⁠. The choice of this energy form will be justified later. Importantly, note that the internal energy is formulated with the elastic strains and not the total strain $εαβT$⁠. The pressure and deviatoric stress are given through differentiation of the internal energy density:

where $σαβ=σαβdev+pδαβ$ is the total stress tensor, and $εs=23εαβdevεαβdev$ is the deviatoric shear strain invariant. It is worth noting the generality of the model above. Typically, $F(ϕ)=ϕd¯$ will be chosen, with $d¯$ as a constant, where this form was motivated by experiments (Rubin & Einav, 2011). For example, linear elasticity is retrieved for n = 0 and $F(ϕ)=1$ (⁠$d¯=0$⁠). Furthermore, if $F(ϕ)=ϕd¯≃(ρρs*)d¯$ is assumed, where $ρs*$ is the unstressed solid density, then depending on the choice of n and $d¯$ the model can recover a variety non-linear elastic models, such as the ones described by Alaei et al. (2021), Riley et al. (2023) and Chen et al. (2023). This versatile structure permits a comprehensive formulation that spans from linear elasticity to pressure-dependent non-linear elastic material, aptly capturing the behaviour typically observed in experiments.

The perturbed stress of equations (12) and (13) can be determined by first-order Taylor series expansion about the homogeneous unperturbed (reference) strain, which gives:

where □* has been dropped from the homogeneous (reference) terms for brevity and the mass balance definition $ϕ˜=ϕε˜v$ was used, assuming no solid density change, as expected in elastic wave experiments.

Consequently, by substituting perturbed strains equation (9) and the perturbed stress equations (14) and (15) into equation (8), and using the definition $ε˜s=23εγζdevε˜γζεs$⁠, the momentum balance for the hyperelastic model with a plane wave perturbation propagating along z can be written as:

where $M¯=K¯+43G¯$ is the constrained modulus, and all state variables are defined in terms of homogeneous unperturbed states. Notably, this result showcases that while the perturbation is along the z-axis, there emerges a dependence on the unperturbed states in the state of other directions, that is $εzβdev$⁠. The wave speed of longitudinal waves is defined as $Vp=ωk$ (Landau et al., 1986; Achenbach, 2012), which gives:

Similarly, by substitution of equation (10) and the perturbed stress equations (14) and (15) into equation (8), the transverse wave speed of a wave propagating along the z-direction and polarised along the x-direction can be found to be

Notably, the longitudinal wave speed in equation (17) and transverse wave speed in equation (18) both depend on strain (or stress states) that cannot be represented solely by invariants, which are a result of the chosen wave perturbation direction and model. Thus, if different propagation or polarisation directions other than the ones illustrated in Fig. 2 were used on the same unperturbed state, the wave speed would naturally depend on different strain components. However, as shown later, only the longitudinal wave speed maintains this dependence on directionality for triaxial conditions.

The particular choice of the internal energy density equation (11) is motivated here. To this end, it will be shown that the tangent (instantaneous) moduli for the simplified scenario of zero shear strain are equivalent to equations (3) and (4). To begin, the constrained tangent modulus is determined by $Mt=Czzzz=∂σzz∂εzz$⁠, which is found to depend on the volumetric strain and solid fraction as follows:

Importantly, the right-hand side of the longitudinal wave speed in equation (17) is equivalent to this modulus in the absence of deviatoric shear components $εxzdev=εyzdev=0$⁠.

In the absence of shear strain, which happens under isotropic stress states whereby the shear stresses are zero, it is readily shown by substituting equation (12) that the constrained modulus simplifies to:

where an additional linear dependence on pressure is noted, a structure that seems to differ greatly from observed empirical trends.

However, by recognising that $p∝F(ϕ)M¯εv11−b$⁠, it is apparent that $F(ϕ)M¯εvb1−b∝pεv$⁠. Furthermore, considering the reasonable and general choice (Rubin & Einav, 2011) of $F(ϕ)=ϕd¯$ the second term reduces to $d¯p$⁠. It is now clear that when assuming small elastic volumetric deformations $εv ≪ 1$ that:

Therefore, the approximate constrained tangent modulus can be expressed through substitution of equation (12):

where p_atm was introduced such that A has units of stress. Similarly, the shear tangent modulus G_t with respect to triaxial shear strain with substitution of equation (12) is:

where again p_atm was introduced such that B has units of stress. Thus, under the assumption of isotropic conditions and small elastic strains (⁠$εs=0,εv ≪ 1$⁠), the tangent moduli reduce to the experimental forms found previously, equations (3) and (4) in the hyperelastic model.

Hypoelastic models tend to adopt the experimentally developed empirical stiffness in equations (3) and (4) without the requirement of being determined from internal energy. Thus, hypoelastic models exhibit path-dependency and do not conform to the thermodynamic requirements of energy conservation (Zytynski et al., 1978; Einav & Puzrin, 2004). The stress rates determined directly from the empirical relations (equations (3) and (4)) are:

where $K¯hypo$ and $G¯hypo$ are the bulk and shear stiffness constants for the hypoelastic model, which have a stress dimension. The superscript emphasises that although these are elastic constants, their values are not inherently equal to those of the hyperelastic model.

Following a procedure similar to the derivation of the hyperelastic wave speeds, the perturbed stress is substituted into equation (8). However, hypoelasticity does not have an explicit expression for total stress $σαβ$ and as such, it was assumed that $σ˜αβ=σ˙αβ$ as this was the result of the first-order Taylor expansion of equations (12) and (13). The resulting longitudinal and transverse wave velocities are

which are equivalent to equations (1) and (2) when equations (3) and (4) are substituted, and thus the empirical relations are recovered. However, the hypoelastic wave speeds lack the additional stress-dependent terms that emerged in the hyperelastic derivation due to the energy conservation of the latter model. These differences are critical to the comparison against experimental results that will be shown between the hyperelastic and the hypoelastic model in triaxial compression. Notably, A and B in equations (26b) and (27b) depend on the elastic constants in a different way than their counterparts in equations (22b) and (23b).

For both the hyperelastic and hypoelastic models, linear elasticity is recovered when b = 0 and $H(ϕ)=1$⁠. In this limit, the longitudinal wave speed and transverse wave speed for both hyperelastic and hypoelastic models reduce to:

where the superscript notation is dropped as the stiffness constants for the hyperelastic and hypoelastic models are equivalent in this case. Therefore, equations (1) and (2) are recovered in the linear elastic limit.

In the isotropic compression scenario (i.e. $εs=εαβdev=0$⁠), it is shown that the hyperelastic wave speeds equations (17) and (18) reduce to the same form as the wave speeds of the hypoelastic model equations (26a) and (27a). If equation (12) while $εs=0$ (isotropic compression) is substituted into equation (17), then the longitudinal wave speed is:

where A is defined in Equation (22b) and the density correction term, that is $ϕpF(ϕ)∂F(ϕ)∂ϕ$⁠, is ignored as discussed previously (cf. equation (21)). Similarly, by substitution of equation (12) into equation (18), the corresponding transverse wave speed for the hyperelastic model is:

where B is defined in equation (23b). The derived wave speeds are equivalent to the hypoelastic model (equations (26a) and (27a)). However, it is crucial to note that A and B are combined parameters of the bulk and shear stiffness constants, $K¯$ and $G¯$⁠, respectively. Therefore, to accurately determine the stiffness constants, both V_p and V_s must be measured for the hyperelastic and hypoelastic models.

Both models are fit to an experimental dataset of wave velocities to elucidate better the equivalence of these derived wave speeds and the empirical relations equations (3) and (4), as shown in Fig. 3. The data were acquired by digitising results presented in Dutta et al. (2021). The wave measurement was performed on Toyoura sand with a specific gravity of $Gs=2·64$⁠, a minimum void ratio of e_min = 0·569, and a maximum void ratio of e_max = 0·936. Notably, V_p and V_s were fit to experimental measurements simultaneously. Despite the sparsity of data, both models fit the experimental results well, and their equivalency in this scenario is clear. The resulting coefficients are shown in Table 1, which will also be used for the following section. The stiffness constants for the hyperelastic and the hypoelastic model are shown in Table 2. Notably, the models have several order of magnitude differences between the elastic coefficients, but it is important to recall that these are from two different models and, therefore, cannot be directly compared with each other. Moreover, the rather large elastic stiffness values of the hyperelastic model occur because once equations (22b) and (23b) are substituted into equations (30) and (31) it becomes clear that the pressure is scaled by $K¯$ as opposed to p_atm as in the hypoelastic model.

While the wave speeds of both models are shown to be mathematically equivalent for isotropic compression, this is not the case for triaxial states. For the hypoelastic model, wave speeds do not change from equations (26a) and (27a) because they solely depend on $ϕ$ and p. However, the hyperelastic model shows that stress dependence is more involved in triaxial states due to being derived from an energy potential. Fig. 2(b) shows the wave perturbations and relevant directions of propagation and particle motion along principal directions in this scenario. For the hyperelastic model, the longitudinal wave speed in triaxial conditions is equivalent to equation (17), but since the reference state is in triaxial conditions, $εzxdev=εzydev=0$ and only $εzzdev≠0$⁠, while the transverse wave speed equation (18), by substitution of equation (12), becomes:

where $η=qp$⁠, and $q=32σαβdevσαβdev$ is the deviatoric shear stress invariant. Notably, the transverse wave speed depends only on the stress invariants; thus, for triaxial states, the wave perturbation direction plays no role in the wave speed when the wave propagates along a principal direction as illustrated in Fig. 2.

Interestingly, both $Vphyper$ and $Vshyper$ can become imaginary when $η=qp=3G¯2bK¯$⁠, which indicates material instability (Rice, 1976; Zhang et al., 2012; Stefanou & Gerolymatou, 2019). Such instabilities have also been documented in other hyperelastic models by evaluating the convexity of the energy (Jiang & Liu, 2003; Agnolin & Roux, 2007; Nguyen & Einav, 2009), where it is important to emphasise that imaginary wave speeds indicate the same loss of convexity. The generation of these imaginary wave speeds corresponds to the development of a positive Lyapunov exponent (Rice, 1976; Sulem & Vardoulakis, 1995), indicating a bifurcation or generation of a localised deformation region (Sulem & Vardoulakis, 1995). Thus, the model directly indicates that the experimental use of elastic wave speeds to determine material properties is no longer valid, as this is often accompanied by assumptions of a homogeneously deformed continuum. This does not preclude the models used to determine analytical wave speeds following this bifurcation but would require a more complicated derivation in which a finite localised region is embedded within the domain of interest as in Zhang et al. (2012) and, thus, would also require knowledge of both density gradients and orientations of the localised region. However, the development of this analytical methodology also implies that experiments would need to accurately determine wave travel times both within and outside the localised region, so that the material parameters could be calibrated appropriately. This approach would be similar to the principles applied for wave speed measurement of layered elastic media wherein it is acknowledged that the domain of interest is non-uniform (Achenbach, 2012). However, the exploration of such scenarios is outside the scope of the current paper.

To understand the differences between the models, wave speed predictions are compared with wave measurements obtained from experimental data, which were digitised from results presented by Dutta et al. (2021). The triaxial experiment conducted was a drained compression test on dry sand, ensuring that the wave speeds would not be influenced by the presence of water (Cho & Santamarina, 2001; Barrière et al., 2012; Leong & Cheng, 2016). The experiments were performed under a radial stress $σr=100$ kPa for three different initial solid fractions: $ϕ0≈0·52,ϕ0≈0·55$ and $ϕ0≈0·60$⁠. It is crucial to highlight that the densest experiment (⁠$ϕ0≈0·60$⁠) experienced shear localisation (Dutta et al., 2021). For lack of information regarding the emergence of the shear localisation, experimental results are shown only up to the peak stress (⁠$εaT≈0·04$⁠) for this experiment, as beyond this point homogeneous deformation cannot be assumed. The prediction has no free parameters at this stage since the stiffnesses obtained from the previous section in isotropic experiments are kept the same, as shown in Table 2.

Figure 4 plots the transverse and longitudinal velocities against the axial strain. The hyperelastic model under-predicts the experimental data for the transverse wave velocities, whereas the hypoelastic model under- or over-predicts experimental data depending on the initial density. However, the hyperelastic model generally offers a more accurate prediction of the data for longitudinal wave speed than the hypoelastic model. Notably, for $ϕ0≈0·60$⁠, the hyperelastic model could not produce a prediction beyond the red dotted line in Fig. 4(c) as these stress states resulted in an imaginary wave speed, suggesting the emergence of an instability or localisation, which in this experiment did occur according to Dutta et al. (2021). However, the hyperelastic model predicted this shear localisation prior to the peak stress of the experiment, where the black arrows denote experimental data points that, according to the model, correspond to localised samples. It is also important to note that the current authors assumed homogeneity until the peak stress due to a lack of sufficient information. However, advanced imaging techniques on triaxial experiments suggest that samples for which localisation occurs experience a loss of a homogeneous deformation prior to this peak stress (Desrues et al., 2018; Amirrahmat et al., 2019). It is also encouraging to emphasise that the model was not calibrated for this experiment but rather used isotropic compression, yet given the material states, it accurately predicted the occurrence of localisation in this experiment. Moreover, while it is indeed always possible to measure wave speeds experimentally, the authors would like to stress that in the context of localised deformation, the standard experimental procedure to extract elastic properties from those measurements is no longer valid.

The observations are further substantiated by Fig. 5, which compares the measured values with those predicted by both models. Encouragingly, both models’ predicted wave speeds are within a range of ± 10% of the measured values, although the hyperelastic model tends to yield better results overall. Both models accurately reproduce shear wave velocities, although the hyperelastic model has an absolute error greater than 5% for $ϕ0≈0·52$⁠, whereas the hypoelastic model maintains an absolute error of less than 5%. Yet, both models perform roughly equivalently for $ϕ0≈0·55$⁠. However, the hyperelastic model outperforms predicting longitudinal wave speeds and essentially stays within the ± 10% error range for $ϕ0≈0·55$ and $ϕ0≈0·60$⁠. The hypoelastic model exceeds an absolute error of 10% in these cases. These results suggest that the additional stress-dependent terms in the hyperelastic model provide correction factors missing in the hypoelastic models. In geotechnical engineering, measurements of V_s are typically more common than V_p due to their empirical correlations with liquefaction properties (Tokimatsu & Uchida, 1990) or in situ measurement techniques such as standard cone penetration (Anbazhagan et al., 2012) or cone penetration tests (Mayne & Rix, 1995), among various other correlations outside the determination of elastic constants (Hussien & Karray, 2015). However, the measurements of both V_p and V_s are equally important for the complete determination of the elastic stiffness constants, as shown by equations (17), (18), (26a) and (27a). While this paper focuses on dry granular media, such media are often accompanied with some level of fluid saturation, resulting in a mixture whose constituents have differing stiffnesses. In such scenarios, V_p and V_s would naturally be affected by the combination of stiffness or compressibilities of the relevant materials in a likely non-linear manner. However, when the fluid is Newtonian (e.g. water or air), its contribution to V_s is likely to be negligible due to the lack of elastic shear stiffness. Future work could investigate multi-phase granular media by incorporating the internal energy corresponding to the compressibility of the constituent media, similar to the approach in Jiang et al. (2017) and Chen et al. (2023).

Lastly, the $Vp/Vs$ ratio is examined to explore shearing-induced differences in the evolution of V_p and V_s. Fig. 6 illustrates the variation of $Vp/Vs$ against $εaT$⁠. Notably, the data reveal an initial increase of $Vp/Vs$ before plateauing, with $ϕ0≈0·60$ not exhibiting a plateau as the experimental results are not shown beyond $εaT≈0·04$ (peak stress). However, the hypoelastic model maintains a constant ratio because the power law coefficients are assumed to be identical for both V_p and V_s, as is required for isotropic materials. One might argue that the power coefficients should differ for empirical relations equations (3) and (4) or include fabric dependency (Dutta et al., 2020; Li et al., 2021), leading to some evolution of $Vp/Vs$⁠, and therefore, extending the empirical models to anisotropic materials. However, the hyperelastic model captures this evolution without additional assumptions (or extra parameters) thanks to the additional stress-dependent terms that arise when deriving the wave speeds. The prediction of evolving $Vp/Vs$ suggests that, at least to a first order, incorporating a more sophisticated stress dependency or fabric might not be essential, as discussed in Mayer & Liu (2010). This offers a compelling methodology to accurately account for stress-induced anisotropy and inherent material anisotropy, enhancing the understanding of prevailing mechanisms within soils.

In this paper, wave speeds were derived for hyperelastic and hypoelastic models that depend on pressure and density in a general manner. The hypoelastic model’s wave speed is equivalent to the empirical relations from Hardin & Richart Jr (1963). The hyperelastic model is also equivalent to those same empirical relations for isotropic compression.

However, for triaxial states, the wave speeds of the hyperelastic model become more elaborate, revealing previously unobserved stress dependencies, which are not present in the hypoelastic model. The hypoelastic model, equivalent to the empirical relations, generally mirrors the overall trends in longitudinal and transverse wave speed. Yet it does not accurately capture the evolution of their ratios $Vp/Vs$⁠. Conversely, the hyperelastic model predicts an evolution in the ratios of the longitudinal and transverse waves due to these additional stress dependencies arising from terms ensuring energy conservation. This suggests that material anisotropy may not be required for a first-order approximation of the evolution of $Vp/Vs$⁠. Furthermore, these intricate stress dependencies might illuminate the observed evolution of V_p under high stresses in oedometric compression (Jia et al., 1999; Makse et al., 2004). More specifically, these findings emerge from not losing energy while waves travel in hyperelastic media, unlike the case of hypoelastic media, and from the consistently derived wave speeds for pressure-dependent elastic media.

Furthermore, the analytical derivation shows that the often used empirical formulations of the type based on Hardin & Richart Jr (1963) for compression and shear waves give the constrained stiffness and shear stiffness, respectively. This indicates that both wave speeds need to be measured to appropriately determine the independent stiffness parameters, that is bulk stiffness and shear stiffness. This distinction is important, as elastic stiffness often serves as a benchmark for design methodologies (Atkinson, 2000; Clayton, 2011) and model calibration (Fahey & Carter, 1993; Ayala et al., 2022; Huang et al., 2023).

In addition, the derived transverse wave speed for hyperelasticity indicates that the wave propagation direction does not affect triaxial states when propagating along principal directions. This suggests that there is a necessity for either fabric or an anisotropic constitutive model to reproduce experimental results for materials that produce this effect (Kuwano & Jardine, 2002; Zamanian et al., 2020; Liu et al., 2022). Furthermore, this procedure could be extended to incorporate other state variables, such as fabric, grain size and shape, as well as saturation. Thus, it provides a methodological and consistent manner in which the effect of additional state variables on wave speed can be considered.

Finally, the hyperelastic model predicts the onset of localisation with the emergence of imaginary wave speeds. This was found to correspond to the emergence of shear localisation in triaxial compression experiments and showcases the predictive capabilities of the model. Notably, while the experiments still produce real-valued wave speeds, the model suggests that those measurements cannot be used directly for the measurement of elastic properties, due to the breakdown of the conventional assumption made in typical experimental procedures of a homogeneously deformed specimen. Thus, the model provides an indicator for when such assumptions are violated and the necessity to consider more advanced interpretation of the experimental conditions by considering layered media and augmenting the analytical derivation as in Zhang et al. (2012).

The authors would like to thank the Australian Research Council for funding through DP190103487.