Nonlinear plane waves in saturated porous media with incompressible constituents

We consider the propagation of nonlinear plane waves in porous media within the framework of the Biot–Coussy biphasic mixture theory. The tortuosity effect is included in the model, and both constituents are assumed incompressible (Yeoh-type elastic skeleton, and saturating fluid). In this case, the linear dispersive waves governed by Biot’s theory are either of compression or shear-wave type, and nonlinear waves can be classified in a similar way. In the special case of a neo-Hookean skeleton, we derive the explicit expressions for the characteristic wave speeds, leading to the hyperbolicity condition. The sound speeds for a Yeoh skeleton are estimated using a perturbation approach. Then we arrive at the evolution equation for the amplitude of acceleration waves. In general, it is governed by a Bernoulli equation. With the present constitutive assumptions, we find that longitudinal jump amplitudes follow a nonlinear evolution, while transverse jump amplitudes evolve in an almost linearly degenerate fashion.


Introduction
Originating in the field of geophysics, the theory of porous media has a long history that goes back to the eighteenth century (see the historical review by De Boer [1]). Poroelasticity theories have also been employed in various biomechanical applications involving the deformation of hydrated porous biological tissues. As noted by Ateshian [2], one early application to biological tissues is the modelling of articular cartilage. More recently, multi-phasic models have been used to model the mechanical response of brain tissue, which is known  infinitesimal Biot waves are recalled. The main results are presented in §4, where nonlinear plane wave motion is considered. Prospective future works are discussed in the conclusion ( §5).

Biphasic mixture theory
We consider an unbounded fluid-saturated porous material, to be described within the framework of the theory of porous media (Biot-Coussy biphasic mixture theory). In what follows, we introduce the main equations describing the motion and the deformation of such a fluidsolid mixture. The solid skeleton is assumed elastic, and heat transfer is neglected. Here, several shortcuts are taken for the sake of conciseness. For more details, the reader is referred to various reference textbooks [27][28][29][30][31] and other related works [2,32].

(a) Kinematics
Consider the position x of a particle. Its components are expressed with respect to an orthonormal basis (e 1 , e 2 , e 3 ) of the Euclidean space, and a Cartesian coordinate system is chosen. The Eulerian position vector x is the same for fluid particles and solid particles, but the reference position vectors X f , X s for the fluid and solid phases are independent. Let n α denote the Eulerian volume fraction of the phase α ∈ {s, f}, i.e. the volume fraction of the solid and fluid constituent in the deformed configuration. The saturation condition requires n s + n f = 1. (2.1) The fluid volume fraction n f corresponds to the Eulerian porosity.
In the Eulerian description of motion, spatial differential operators computed with respect to x are written div, grad, etc. In the skeleton-Lagrangian description of motion, the spatial coordinate is the position X s of a skeleton particle in its reference (undeformed) configuration. Spatial differential operators computed with respect to X s are written Div, Grad, etc.
We introduce the deformation gradient tensor F = Grad x of the solid phase, as well as its inverse A = F −1 given by A = grad X s . By introducing the displacement field u s = x − X s of the solid phase, the deformation gradient tensor and its inverse are rewritten F = I + Grad u s and A = I − grad u s , (2.2) where I = [δ ij ] is the metric tensor, here represented by Kronecker delta components δ ij . In related works, the tensor A is called the distorsion tensor [33,34]. Various strain tensors can be defined as functions of F or A, such as the right Cauchy-Green deformation tensor C = F T F and the Green-Lagrange strain tensor E = 1 2 (C − I). To describe the mixture's motion, we introduce the velocity fields v α = x α . The prime with index α denotes the particle time derivative, which is computed while following the motion of the solid (α = s) or of the fluid (α = f). Thus, for any scalar Eulerian field Γ (x, t), we have and similar differentiation operators can be introduced for vectorial and tensorial fields. The velocity fields may be rewritten as v α = (u α ) α , where u α = x − X α denotes the displacement from a particle of α from its reference position X α to its current position x.
We also introduce the Eulerian velocity gradients L α = grad v α and their symmetric part D α = 1 2 (L α + L αT ). In the solid phase, the velocity gradient L s = F s F −1 depends on the deformation gradient and its material derivative, and so does the symmetric part D s = F −T E s F −1 . Finally, to describe fluid motion with respect to the skeleton, we introduce the seepage velocity

(b) Balance principles
Let us neglect mass transfer and external mass supply. We introduce the mass densities ρ α = n α ρ αR , where ρ αR are the real mass densities. Also known as true, intrinsic or effective material density, ρ αR represents the mass of a constituent per volume of that constituent.
In this study, we consider an incompressible skeleton saturated by an incompressible fluid. By definition, the true density ρ αR of an incompressible constituent is invariant. Under this assumption, the Eulerian mass continuity equation for each constituent α ∈ {s, f} reads (n α ) α + n α div v α = 0, (2.4) where the material derivative is defined in equation (2.3). By summation of both continuity equations, the saturation constraint (2.1) yields the condition div n f v f + n s v s = 0, (2.5) which will be used to enforce saturation later on. Introducing the volume dilatation J = det F of the solid phase, the relation J = n s 0 /n s is obtained by integration of the solid mass balance equation of equation (2.4), where n s 0 = 1 − n f 0 denotes the volume fraction of the solid phase in the reference configuration. Therefore, porosity n f = 1 − n s 0 /J is a function of the deformation. Contrary to monophasic incompressible solids that support only isochoric deformations (i.e. J ≡ 1 is prescribed), the volume dilatation of biphasic particles increases with porosity, in a similar fashion to a squeezed sponge.
The balance of linear momentum equation for each constituent reads where the vector b α denotes external body forces per unit mass. The reciprocity conditionp s = −p f follows from the balance of linear momentum applied to the mixture as a whole. Assuming microscopically non-polar constituents, the symmetry of the partial Cauchy stresses σ α = σ αT is deduced from the balance of moment of momentum, where no supply of momentum is included.
In the absence of heat flux (adiabatic case), energy transfer and external energy sources, the local form of the balance of internal energy for each constituent reads where e α denotes the specific internal energy. The colon denotes the double contraction of secondorder tensors. Introducing the density of mechanical energy E α = e α + 1 2 v α 2 , one may express the balance of mechanical energy as for each constituent. Here, σ α v α represents the Poynting vector, and ρ α v α · b α is the work done by the external body force. For consistency with Biot's linear theory of saturated porous media, some adjustments have to be made [28]. Indeed, the local balance of energy (2.7) over the fluid phase α = f does not include the 'tortuosity' effect of Biot's theory, which cannot be captured by the macroscopic mixture approach. Following estimations at the scale of a representative volume, Coussy [28] introduces the tortuosity vector a = (a − 1) w f that modifies the balance of energy (2.7)-(2.8) for the fluid phase as follows: The tortuosity factor a ≥ 1 is defined as the ratio between the seepage energy averaged over an elementary volume, and the corresponding macroscopic quantity ρ f w 2 . It satisfies a → 1 in the single-constituent fluid limit, and a → +∞ in the single-constituent solid limit. As noted by Wilmanski [35], objective relative accelerations may be introduced instead of the relative acceleration w f to account for the tortuosity effect. We will see later on that the tortuosity vector  of equation (2.9) adds −ρ f a to the interaction forcep f , leading to a Cattaneo-type effect on the filtration law. The simple mixture model without tortuosity effect is recovered by setting a ≡ 1.

(c) Constitutive modelling
In contrast to the above balance principles that are written for each constituent, the postulate of entropy increase is written for the mixture as a whole. In a standard way, we consider a singletemperature mixture such that θ > 0 is the temperature field for all the constituents, and η α are the specific entropies. Thus, the second principle of thermodynamics is expressed by the Clausius- where D is the dissipation in the mixture. The thermodynamic procedure based on the temperature is well-described in the literature [27,31]. To model nearly isentropic processes such as acoustic perturbations [36], one may assume that the biphasic mixture is described by the state variables {η s , η f , E}, where E is the Green-Lagrange strain tensor. Moreover, because we are considering a constrained mixture with incompressible constituents, various simplifications can be performed [31]. Here, phase separation is assumed, which amounts to stipulating that e s is a function of {η s , E}, and that e f is a function of η f only. Thus, according to the Gibbs identity, the total material derivatives of the functions of state e α satisfy The intrinsic incompressibility of the constituents is introduced using the method of Lagrange multipliers. For this purpose, the differential form (2.5) of the saturation constraint is expanded as follows using vector calculus identities: p n s div v s + n f div v f + (grad n f ) · w = 0, (2.12) where the corresponding Lagrange multiplier p has been introduced. Using the conservation of energy (2.7), summation of the above equations (2.10)-(2.12) yields the final expression of the dissipation. Due to the tortuosity effect of equation (2.9), the dissipation becomes (2.13) where θ = ∂ e α /∂η α is required to ensure the positivity of the dissipation for arbitrary transformations. Here, we have used the reciprocity conditionp s = −p f , and we have introduced the strain energy density function W = ρ s 0 e s of the skeleton, with ρ s 0 = n s 0 ρ sR . Following standard arguments, the dissipation inequality (2.10) is satisfied for arbitrary processes if Known as the effective drag force, the quantityp f e entails no dissipation if orthogonal to w. Using the saturation condition (2.1), Terzaghi's effective Cauchy stress reads where σ i = σ s + σ f denotes the inner part of the mixture stress. As explained by Carcione [24], 'the effective-stress concept means that the response of the saturated porous medium is described by the response of the dry porous medium with the applied stress replaced by the effective stress'. The remaining dissipation D in equation (2.14) is ensured positive by settinĝ (2.16) which models the internal friction between solid and fluid. The parameter k f ≥ 0 is the permeability of the fluid, i.e. the ratio of the skeleton's intrinsic permeability and the fluid's dynamic viscosity. It satisfies k f → +∞ in the single-constituent fluid limit, and k f → 0 in the single-constituent solid limit. Injecting the expression ofp f e in the conservation of momentum equation (2.6) for the fluid constituent, one eventually obtains Darcy's filtration law, (2.17) More general forms of Darcy's law may include a permeability tensor instead of the scalar k f .

Remark.
In the case of fluid flow through a rigid porous skeleton, the porosity n f is constant. The interaction force of equations (2.14)-(2.16) becomeŝ where we have used the definition a = (a − 1) w f , and τ a = (a − 1)k f ρ fR /n f is a characteristic time.
Thus, we note that the tortuosity effect of equation (2.9) yields a Cattaneo-type relaxation in the filtration law. Again, one may have replaced the present relative acceleration w f by an objective derivative [35], e.g. in a similar fashion to the so-called Darcy-Jordan-Cattaneo model [37].
We assume that the skeleton's effective mechanical response (2.15) follows from the two-term Yeoh strain energy function where I 1 = tr B is the first principal invariant of the left Cauchy-Green tensor B = FF T . The corresponding constitutive relation reads (2.20) where the positive constants λ, μ are the Lamé parameters of linear elasticity. The Yeoh parameter β ≥ 0 has been introduced for the sake of generality, in view of discussing the influence of the constitutive assumptions. With this choice, the neo-Hookean model β → 0 used by Diebels & Ehlers [18] is recovered as a special case. While the analysis introduced hereinafter is quite general, most of the exact analytical formulae are obtained in the neo-Hookean limit. The more general case β > 0 is addressed in a quasi-analytical fashion. Moreover, we assume that the fluid's permeability k f follows from the formula [38] which is an alternative to the Kozeny-Carman formula of [28,32]. Here, k f 0 represents the fluid's permeability when the porosity n f equals its initial value n f 0 , and κ is a dimensionless parameter. Lastly, we assume that the tortuosity coefficient a satisfies Berryman's formula [24,28,35] but more general expressions could be used [38]. For the purpose of illustration, typical values of the material parameters for a soft biological tissue saturated by an incompressible liquid are specified in table 1. The elastic parameters λ, μ are deduced from Comellas et al. [6], 1 Table 1. Physical parameters of water-saturated brain tissue inferred from [6], where the mass density of water at room temperature is assumed for both constituents α ∈ {s, f}. The potential mismatch between isothermal and isentropic measurements is neglected in the present study. shear stresses are consistent with [6] over a large range of deformations (simple shear strains ranging from −0.7 to 0.7). The tortuosity coefficient (2.22) deduced from the reference value of the porosity n f 0 is a = 3.0.

(d) Eulerian equations of motion
In the Eulerian specification of motion, 2 spatial differential operators are computed with respect to x. We consider a fluid-saturated poroelastic material with incompressible constituents governed by equations (2.4)-(2.6). In addition, a kinematic relationship between the distorsion tensor A defined in equation (2.2) and the velocity v s of the solid phase is introduced in the first line of equation (2.23) below. The latter can be retrieved by using the equality of mixed partials in equation (8.3) of Godunov & Romenskii [34]. Thus, the equations of motion read as a system of balance laws constrained by the saturation condition of equation (2.1). This system is closed by the constitutive equations for the partial stresses σ α and the interaction forcesp f = −p s in equations (2.14)-(2.16). We therefore end up with a system of 18 scalar equations, which involves the 18 components of {A, n α , v α , p} for α ∈ {s, f}.
Keeping equations (2.4)-(2.6) for the fluid phase and adding the latter to the equations for the solid phase, we may rewrite the above system as which introduces the mixture's inner stress σ i = σ e − pI, see equation (2.15), effective body force ρb = ρ s b s + ρ f b f , and effective density ρ = ρ s + ρ f . We thus end up with a system of 17th scalar equations, which involves the 17 components of {A, n f , v s , w, p}, where w = v f − v s is the seepage velocity. For sake of exhaustiveness, let us mention that appropriate boundary conditions should be provided [28]. Here, plane waves propagating in unbounded domain are considered. Note that standard vector calculus identities can be used to derive alternative forms. In particular, the second line of equation (2.23) might be removed, since the porosity n f is a function of A (consequence of equation (2.4)). Moreover, the last line of equation (2.23) may be rewritten in more compact form by introducing the mixture velocity and the mixture stress tensor, see [2,30]. Contrary to the case a ≡ 1 of simple mixtures, no fully conservative first-order formulation of the equations of motion is known due to the dependency of the tortuosity coefficient (2.22) with porosity. When a is uniformly equal to unity, the above system is analogous to the equations in  [18,19]. Note that equation (2.23) is not straightforwardly linked to the nonlinear Biot theory by Grinfeld & Norris [25], where different inertial terms are proposed.

Biot's theory (a) A linearization
Let us assume that the effective stress in the solid phase satisfies Hooke's Law of linear elasticity σ e = λ tr(ε)I + 2με, where ε = 1 2 (grad u s + grad T u s ) is the infinitesimal strain tensor, and λ, μ are the Lamé constants. When we linearize the equations of motion (2.23) about an undeformed static state by neglecting convection terms, we have where the porosity n f = n f 0 , fluid permeability k f = k f 0 and tortuosity a are constant parameters deduced from the values in table 1. Equation (3.1) corresponds exactly to the low-frequency Biot equations with incompressible constituents [23], for which Biot's effective-stress coefficient β (or α [24]) equals unity and the other Biot parameter M becomes infinite. A more general linearization about arbitrary pre-deformations in a small-on-large fashion would lead to the acousto-elastic equations, see for instance Grinfeld & Norris [25].

(b) Harmonic plane waves
We recall the main dispersion characteristics of this theory hereinafter. To do so, harmonic planewave motion is assumed by setting the space-time dependence of the unknowns to e i(ωt−k ω x) , where ω is the angular frequency, k ω is the wavenumber and i = √ −1 is the imaginary unit. In the absence of body forces b α = 0, non-trivial solutions to equation (3.1) are obtained provided that one of the following dispersion relationships is satisfied: with the coefficients Note that for the simple mixture model a ≡ 1, the above result is the same as that given by De Boer & Liu [15] if fluid compressibility is neglected therein. In equation (3.2), the first family of linear waves with elastic modulus λ + 2μ corresponds to longitudinal compression waves (P) resulting from the interaction of both phases, while the second family with elastic modulus μ corresponds to transverse shear waves (S) mostly supported by the solid skeleton. Recall that although each phase is incompressible, the volume of a particle of mixture can change if the relative quantity of its constituents is modified, i.e. if the porosity is not kept constant (see the comments following equation (2.5)). Figure 1 represents the frequency evolution of the phase velocity v ω = ω/Re k ω and of the attenuation coefficient α ω = −Im k ω , for waves propagating towards increasing x. The horizontal dotted lines mark the respective high-frequency (or inviscid-fluid) phase velocity limits. As noted by Coussy [28], 'the undrained situation is recovered (. . . in) the low-frequency range', where only S-waves propagate.  In the high-frequency range, the analysis reveals that compression and shear waves propagate with strong attenuation resulting from the interaction between the phases. However, this frequency range is not well-described by the present theory, as viscous dissipation inside the fluid phase becomes preponderant [23]; for complements and extensions, see the literature on homogenization theory [41] and enriched continuum models [42]. Thus, this model is valid in the low-frequency range ω ω c , where f c = ω c /(2π ) is a characteristic frequency given in equation (3.3). In practice, the parameter values of table 1 yield f c ≈ 120 MHz (vertical line in figure 1), which means that the present set of parameters provides a valid model up to ultrasonic frequencies.
Of course, the value of the critical circular frequency ω c plays a crucial role (note in passing that its order of magnitude is related to the characteristic time τ a of the Cattaneo-type filtration law (2.18)). If we use the parameter values of the study by Hosseini-Farid et al. [5] instead (see also Forte et al. [4]), i.e. n f 0 = 0.17 and k f 0 = 1.61 × 10 −12 m 2 /(Pa · s), then we find f c ≈ 4.7 MHz, which has similar orders of magnitude to the value obtained previously.
Since the propagation characteristics of both waves depend on the linear Biot parameters, acoustic experimental measurements within an adequate frequency range could provide dynamic estimations of those parameters. Ideally, the frequency range of interest should be chosen high enough for the slow P-wave to propagate, but low enough to limit attenuation [23].

Nonlinear plane waves
Now we analyse the characteristics of various types of nonlinear wave solutions to equation (2.23) where body forces b α are neglected. We restrict the study to a one-dimensional configuration, assuming invariance along the y-and z-directions.
Because the motion does not depend on y and z, the deformation gradient and distortion tensors of equation (2.2) are of the form where J = det F denotes the volume dilatation given by J = F 11 = (A 11 ) −1 . Therefore, the motion includes possibly a volume-changing compressive deformation (11-component) and a volumepreserving shear deformation (21-and 31-component).
With the above invariances in mind, the system (2.23) can be rewritten as a first-order quasilinear system of partial differential equations, where the coefficients of M ν for ν ∈ {t, x} and R are specified in appendix A. Note in passing that these arrays do not depend on v s 2 , v s 3 and p. Moreover, they are not symmetric, and the matrix M t is not invertible.

(a) Smooth travelling waves
We consider plane wave solutions propagating with constant speed c, such that the field variables are smooth functions of ξ = x − ct only. Hereinafter, primes denote differentiation with respect to ξ , so that ∂ x = (·) and ∂ t = −c (·) according to the chain rule. Hence, our system (4.2) reduces to the ordinary differential system Attempts to exhibit such solutions numerically have been unsuccessful up to now, suggesting that travelling waves may not propagate in such a material. Viscous dissipation, compressibility or compaction might be needed for this peculiar nonlinear phenomenon to emerge [22,37,43].

(b) Characteristic wave speeds
Let us consider particular wave solutions for which the matrix M(q) is singular. To do so, let us focus on the homogeneous system M(q) q = 0 by setting R equal to zero. According to equation (A 1) of the appendix, this amounts to assuming inviscid flow for which k f → +∞. Nontrivial solutions for q can be obtained if det M(q) vanishes, restricting the value of the wave speed c to one of the generalized eigenvalues of M x and M t .
In general, it is a difficult task to compute these characteristic wave velocities analytically. If the material has a neo-Hookean behaviour (β = 0), then we find that the wave speed c equals one of the following values:   Acoustic waves propagate, i.e. hyperbolicity is ensured, if the sound speeds in equation (4.4) are real. For this purpose, the radical's argument in the expression of c ± P and c ± S must be nonnegative. In both cases, one notes that this quantity is of the form a + bw 2 1 . The propagation condition a + bw 2 1 ≥ 0 can be simplified if the coefficient b is non-negative, in which case imposing a ≥ 0 will be sufficient to ensure hyperbolicity.
While the analysis of hyperbolicity is straightforward for shear waves with speed c ± S , in which case b is non-negative, the analysis is less obvious for compression waves with speed c ± P . In the case of simple mixtures (a ≡ 1), the coefficient b in the expression of c ± P has the same sign as (n f − 1); hence, it is always negative. If Berryman's formula (2.22) is used instead (a ≡ 1), this coefficient has the same sign as (n f − 1)(n f + 1 2 ) + , where the constant = 1 16 ρ fR /ρ sR depends on the ratio of the reference densities. In the present low-porosity material with ρ fR ρ sR (see the values in table 1), the coefficient b in the expression of c ± P is negative. Thus, the propagation of compression waves requires that the seepage velocity has a moderate amplitude |w 1 | along the direction of propagation.   Using Berryman's formula (2.22), the coefficients in equation (4.6) satisfy ϑ → +∞ and θ → 1 at zero porosity. At unit porosity, they satisfy ϑ, θ → 0. Thus, as shown in the figure (solid lines), compression waves do not propagate at zero porosity, and both waves do not propagate at unit porosity. The first remark relates to the fact that the monophasic solid limit n f → 0 is an incompressible solid in which shear waves propagate, but not poroelastic compression waves. The second remark expresses the fact that the monophasic fluid limit n f → 1 does not support shear stresses or poroelastic compression. Figure 2 compares the evolution obtained for variable tortuosity (2.22) with that of the simple mixture model where the tortuosity coefficient a ≡ 1 is not porosity-dependent. One notes that the tortuosity effect has a significant influence on the P-wave velocity at low porosities, while the shear-wave velocity does not seem to be significantly affected by this feature, see also Wilmanski [35] where comparisons between Biot's theory and the simple mixture model are proposed. In a different context, the fact that the tortuosity factor has a major influence on the speed of the slow P-wave is a well-known feature [44].

(c) Perturbation approach
Let us investigate the influence of the Yeoh parameter β on the characteristic speeds by using a perturbation method [45]. For this purpose, we introduce pairs of left and right generalized eigenvectors l, r deduced from the condition M(q) q = 0. The vectors r form a basis of the right null space of M, while the vectors l form a basis of the left null space of M, i.e. they belong to the right null space of M T .
As shown in the matrices' expression (appendix A), the matrix M x is linear in β, but M t does not depend on β. This can be rewritten as a perturbation of the form M x = M x0 + βM x1 , where the zeroth-order matrix M x0 corresponds to the neo-Hookean case discussed in the previous section. Thus, we seek generalized eigenvalues and eigenvectors as power series of β: c = c 0 + βc 1 + · · · , l = l 0 + βl 1 + · · · , r = r 0 + βr 1 + · · · , (4.7) where the zeroth-order quantities c 0 , l 0 , r 0 correspond to the case of neo-Hookean behaviour at first order in β. One observes that the increment of the speed of sound is linear with respect to the (presumably small) perturbation βM x1 of the matrix M x . In practice, the pairs of vectors l 0 , r 0 deduced from the previous section are normalized in such a way that l 0T M t r 0 = 1, which greatly simplifies equation (4.9). Compression waves. Let us go back to the zeroth-order neo-Hookean case. Using a computer algebra system, one pair of vectors l 0 , r 0 is deduced from M(q) by solving the generalized eigenvalue problem corresponding to the characteristic speed c 0 P = c + P of equation (4.4). The components of these vectors lead to the perturbation (4.9) of the speed of compression waves. The quantity Q 1 11 denotes the first-order increment of the coefficient Q 11 given in appendix A. Note that the speed of sound is no longer exclusively a function of volume-changing strain A 11 , and that the above perturbation has a non-zero value in the undeformed state.
of the shear wave speed with polarization along y, where the coefficient Q 1 22 is deduced from the appendix A. One notes that the speed of sound is no longer independent of the shear deformation A 21 , A 31 and that this dependency is quadratic, which is coherent with related studies [13]. Here, the sound speed in an undeformed state (4.6) is unchanged. An expression similar to equation (4.11) is found for shear waves polarized along z, where the increment Q 1 22 needs to be replaced by a coefficient Q 1 33 obtained in a similar fashion from the expressions in the appendix. Figure 3 illustrates the validity of the above perturbations. Figure 3a compares the perturbation (4.10) of the Yeoh P-wave speed (blue dashed line) with the same value obtained by numerical resolution of the generalized eigenvalue problem of M x and M t (blue solid line). The value of the perturbation parameter β = 2.2 is taken from table 1, as well as the values of other parameters. Here, a static state under pure dilatation is considered, i.e. the shear strain components A 21 , A 31 are set to zero while A 11 is varied. Thus, the porosity n f = 1 − 0.8 A 11 is not constant. The agreement between both curves is very good in the vicinity of the static undeformed state A 11 1. Similarly, figure 3b illustrates the effect of the perturbation (4.11) on the shear wave speed. Here, a static state under simple shear is considered, i.e. A 11 = 1 and A 31 = 0 are imposed while A 21 is varied, and the porosity n f = 0.2 is constant.
In both figures, the black solid line corresponds to the neo-Hookean case (4.4) where β = 0, and the vertical dotted line marks the undeformed state. In figure 3b, one observes that the situation is almost symmetric with respect to the undeformed state, and that the neo-Hookean case yields a strain-independent shear wave speed. The picture is different in figure 3a, where the curves are not symmetric with respect to the undeformed state, and where the neo-Hookean model yields already strain-dependent sound velocities. These observations suggest that the nonlinearity of poroelastic P-wave propagation is of a very different nature to that of shear wave propagation.

(d) Acceleration waves
Similarly to [16,21,35], let us analyse the speed and evolution of acceleration waves. For such wave solutions, the primary field q is continuous across the surface ξ (x, t) = 0 with ξ = x − s(t), but its normal derivative ∂ ξ q may be discontinuous. Typically, such solutions represent situations in which the field variables experience a brutal change of slope; for instance, an initial-value problem with piecewise linear initial data. As proposed by Müller & Ruggeri [46], we assume that the wave propagates into a domain where the primary field q is a constant equilibrium stateq of equation (4.2), for which the seepage velocityw equals zero; more general cases are discussed in the literature [46]. The jumps [[ · ]] of the partial derivatives across the moving surface are related to those of the normal derivative ∂ ξ q according to [ If r is scaled in such a way that it has the same dimension as q componentwise, then the wave amplitude Π is expressed in m −1 . Now we derive Bernoulli's evolution equation satisfied by the wave amplitude following §8.4 of [46]; a similar result was obtained by Ciarletta et al. [21] for general Biot-like models. For this purpose, we consider a vector l belonging to the left null space of M(q), and such that l T M t r = 1. As shown in the literature [46], a Bernoulli differential equation is obtained where d/dt denotes the directional derivative ∂ t + c ∂ x along the curve that follows the position of the surface. A well-known analytical solution to equation (4.13) yields the time-evolution, , (4.14) of the jump amplitude as the wave propagates, in terms of the coefficients (4.15) evaluated at the constant equilibrium stateq. Assuming positive coefficients Ω 1 , Ω 2 in the expression (4.14) of the wave amplitude, we observe that the denominator vanishes at some finite time, (4.16) if the initial jump Π (0) is smaller than −Ω 1 /Ω 2 . Conversely, such acceleration wave solutions are stable for positive times under the condition Π (0) > −Ω 1 /Ω 2 . Note that Ω 1 vanishes in the case of inviscid flow k f → +∞. Since Ω 1 depends on R, it accounts for attenuation. As can be seen from equation (4.14), the constant Ω 1 is responsible for the decay of the jump amplitude, and therefore provides a smoothing effect on wave solutions. The constant Ω 2 vanishes when the characteristic speed c corresponds to a linearly degenerate eigenspace [47]. Thus, this constant expresses the nonlinearity of wave propagation, and therefore may yield a competing wavefront steepening effect.
Compression waves. Assume thatq = [1, 0, 0, n f 0 , 0, 0, 0, 0, 0, 0,p] T corresponds to a motionless undeformed equilibrium state, and that the material's behaviour is neo-Hookean (β = 0). Using a computer algebra system, one pair of vectors l, r is deduced from M(q) by solving the generalized eigenvalue problem corresponding to the characteristic speed c = c + P ≈ 4.99 m s −1 , which is the value displayed in figure 2 and Shear waves. Again,q = [1, 0, 0, n f 0 , 0, 0, 0, 0, 0, 0,p] T corresponds to a motionless undeformed equilibrium state, and the material's behaviour is assumed neo-Hookean (β = 0). For the characteristic speed c = c + S ≈ 2.71 m s −1 , two distinct pairs of vectors l, r are found, corresponding to shear waves polarized along y or z. For both polarizations (i ∈ {2, 3}), we find that the acceleration jumps are linked through where the tortuosity coefficient a is deduced from the porosity in the stateq. This relationship suggests that acceleration S-waves propagate both in the fluid and the solid phase if a = 1.
, Ω 2 = 0 (4.20) evaluated atq. Therefore, transverse acceleration waves decay exponentially, consistent with the study by De Boer & Liu [17] where a ≡ 1. The characteristic distance of decay deduced from table 1 is c + S /Ω 1 ≈ 101 nm. Note that the order of magnitude of this characteristic distance relates to that of the high-frequency attenuation distance α −1 ω deduced from the dispersion analysis (figure 1). Figure 4 displays the time evolution of the amplitude Π deduced from equation (4.14), for nonlinear poroelastic acceleration P-waves and S-waves propagating in a neo-Hookean material (β = 0). While shear wave amplitudes decay exponentially (smoothing effect), compression waves are subject to a critical amplitude −Ω 1 /Ω 2 deduced from equation (4.18) (horizontal dashed line in figure 4a), below which the solution becomes infinite in finite time. As stated in Müller & Ruggeri [46, p. 183], 'if the initial discontinuity in the derivatives is too strong, it cannot be damped; instead it grows to infinity and thus the acceleration wave develops into a shock wave'. Beyond the critical amplitude, nonlinearity overpowers attenuation effects, leading to the formation of shock waves.
Of course, the addition of Yeoh behaviour with β > 0 modifies the picture slightly (see figure 3, and the modified sound velocities in equation (4.9)), since the shear sound velocities are no longer independent on the shear deformation. Therefore, the linear degeneracy property for shear waves will potentially be lost. Nevertheless, since the shear wave speed is quadratic in the shear strains, the values found for Ω 1 , Ω 2 should not be greatly affected by this modification, at least about an undeformed state where the sound velocity is nearly constant. In a small deformation range, the significant difference of magnitude for the coefficients Ω 1 , Ω 2 in compression and shear waves will remain a predominant feature of the material. Therefore, shock waves will still develop more easily in compression than in shear, while poroelastic P-waves are subject to faster smoothing than S-waves due to stronger attenuation.

Conclusion
A mixture-theoretic Biot model for large deformations in incompressible media has been presented, in view of future biomechanical applications. Here, saturated Yeoh-type porous solids were considered. The main features are the existence of shear waves and slow compression waves, which linear dispersive properties follow from the Biot theory. The computation of plane-wave solutions with discontinuous gradients shows that shear jumps decay exponentially (in a similar fashion to the linear theory), while the compression jumps are governed by a nonlinear Bernoulli equation. Thus, in the neo-Hookean limit, large compressive jumps can lead to the formation of shocks, which is not the case for shear jumps.
These results can be used for the validation of numerical methods [18,20]. Moreover, the modelling framework and the methodology are applicable to other fields, for instance where wave propagation problems in compressible or triphasic mixtures arise. As discussed above, shock waves may emerge. In this regard, a first difficulty lies in the quasi-linear (non-conservative) form of the equations of motion, for which the definition of shock wave solutions is not straightforward [46]. While this problem can be circumvented in the case of simple mixtures a ≡ 1, the derivation of a conservative form is less obvious in the case of porosity-dependent tortuosity coefficients a ≡ 1. Nevertheless, shock wave solutions can still be investigated numerically. A possible strategy would be to rely on shock-capturing finite volume methods [47], e.g. based on an 'artificial compressibility' approach to account for the saturation constraint.
As far as the present problem is concerned, several improvements need to be mentioned. First, one should be aware that the generality of the results is limited by the constitutive assumptions, but that the same approach could be used for variations of this model. Second, the use of poroelasticity in applications requires the experimental determination of relevant model parameters (table 1). In practice, the brain mechanics literature suffers from a lack of experimental data in dynamic configurations, which would be representative of head trauma configurations. Lastly, several modelling refinements could be introduced in potential fine-tuning steps, such as viscoelastic behaviour [5,6], non-Darcy flow [38] or objective derivatives [35], to name a few. A conclusive experimental or computational assessment of multi-phasic effects in TBI is still needed.
Data accessibility. This work does not have any external supporting data. Competing interests. The author declares that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Funding. This work was supported by the Irish Research Council (project ID GOIPD/2019/328). Acknowledgements. The author is grateful to Michel Destrade (NUI Galway) for fruitful discussions, careful reading and support.

Appendix A. System matrices
Using the invariance assumption along y, z, the deformation gradient and distorsion tensors can be simplified, see equation (4.1). Moreover, the system (2.23) becomes for indices i ranging from one to three, where Einstein's notation for repeated indices was used.
Here, we have used the notation ρv = ρ s v s + ρ f v f for the mixture momentum, where ρ = ρ s + ρ f denotes the mixture density. Let us rewrite this system in quasi-linear form (A 1). To do so, we expand the spatial derivatives by using the product rule. By virtue of the chain rule, the spatial derivatives of the stress components σ e become ∂ x σ e i1 = −Q ij ∂ x A j1 with the coefficients Q ij = −∂σ e i1 /∂A j1 . Considering the unidimensional deformation defined in equation ( the quasi-linear first-order system (4.2), with the 11 × 11 matrices M ν and the vector R specified below.