A generalized effective anisotropic poroelastic model for periodically layered media accounting for both Biot's global and interlayer flows

We present a generalized effective poroelastic model for periodically layered media in the mesoscopic scale range, which accounts for both Biot's global and interlayer wave‐induced fluid flow, as well as for the anisotropy associated with the layering. Correspondingly, it correctly predicts the existence of the fast and slow P‐waves as well as quasi and pure S‐waves. The proposed analytical model is validated through comparisons of the P‐wave and S‐wave phase velocity dispersion and attenuation characteristics with those inferred from a one‐dimensional numerical solution of Biot's poroelastic equations of motion. We also compare our model with the classical mesoscopic model of White for a range of scenarios. The results demonstrate that accounting for both wave‐induced fluid flow mechanisms is essential when Biot's global flow prevails at frequencies that are comparable or smaller with respect to those governing interlayer flow. This is likely to be the case in media of high permeability, such as, for example, unconsolidated sediments, clean sandstones, karstic carbonates, or fractured rocks. Conversely, when interlayer flow occurs at smaller frequencies with respect to Biot's global flow, the predictions of this model are in agreement with White's model, which is based on quasi‐static poroelasticity.

the theory of poroelasticity, the porous medium is assumed spatially homogeneous at the macroscopic scale. In the presence of heterogeneities in the mesoscopic-scale range, i.e., at scales much larger than the typical pore size but much smaller than the dominant wavelength, WIFF is caused by fluid pressure gradients related to stiffness contrasts associated with the inhomogeneities (e.g., Müller, Gurevich, and Lebedev 2010). Mesoscopic WIFF can be effectively analysed through Biot's theory by relaxing the inherent assumption of homogeneity at the macroscopic scale (e.g., Rubino, Ravazzoli, and Santos 2009). In close analogy to mesoscopic WIFF, microscopic WIFF (or squirt flow) is associated with the fluid pressure gradients arising between the stiffer bulk part of the pore space and its more compliant components, such as, for example, open grain contacts or microcrack (e.g., Chotiros and Isakson 2004;Gurevich et al. 2010;Müller et al. 2010). The range of validity of common wave propagation theories in response to the topology of the pore space was explored by Sarout (2012).
The interlayer flow in a periodic sequence of porous layers represents one of the simplest forms of mesoscopic WIFF amenable to analytical solutions. This scenario is, however, of particular relevance as layering is arguably the most common heterogeneity in sedimentary environments. In their pioneering work, White, Mikhaylova, and Lyakhovitskiy (1975) derived the effective complex-valued and frequencydependent P-wave modulus for wave propagation normal to an alternating sequence of gas-and water-saturated layers. Important contributions of various authors have confirmed and refined these theoretical results through different methodological approaches (e.g., Norris 1993;Gurevich and Lopatnikov 1995). Gelinsky and Shapiro (1997a) derived the low-and high-frequency limits of the effective stiffness matrix, which accounts for the anisotropy induced by the layering, thus generalizing the well-known Backus average (Backus 1962) to poroelastic media. The frequency dependence linking the low-and high-frequency limits of the stiffness coefficients was then derived by Krzikalla and Müller (2011). The resulting model is valid for any angle of incidence and accounts for the attenuation and velocity dispersion of compressional (P) and for the so-called quasi (SV) and pure (SH) shear waves.
The characteristic frequency ω meso corresponds to the frequency of maximum production of attenuation associated with the mesoscopic WIFF and separates the low-frequency limit, where the fluid pressure in the sample is relaxed, from the high-frequency limit, where the fluid pressure in the sample is unrelaxed. This frequency depends on the diffusivity parameter D and on the characteristic size of the mesoscopic heterogeneities L meso where κ 0 , M, and η denote the permeability, the fluid storage modulus, and the fluid viscosity; and P d and P U are the drained and undrained P-wave moduli, respectively. A common assumption inherent to all of the WIFF models at the mesoscopic scale is that they are valid in the quasi-static regime, which in turn implies that the characteristic frequency ω meso associated with the considered mesoscopic WIFF is much lower than Biot's characteristic or critical frequency, i.e., where ρ f , φ, and τ denote the pore fluid density, the porosity, and the tortuosity of the porous matrix. In the proximity of this frequency, the maximum attenuation due to Biot's global flow is produced. Given that the permeability tends to dominate Biot's characteristic frequency, the assumption of quasistatic conditions is inherently valid for many, if not most, low-permeability formations for the frequencies considered in most practical geophysical applications. However, there are also practically important scenarios involving formations of intermediate to high permeability, such as, for example, unconsolidated sediments, clean sandstones, karstic carbonates, or fractured rocks, which may require accounting for both Biot's global and interlayer flows (e.g., Pride, Berryman, and Harris 2004;Kudarova, van Dalen, and Drijkoningen 2013). Gelinsky and Shapiro (1997b) derived a dynamicequivalent medium approach which captures both WIFF mechanisms, i.e., Biot's global and interlayer flows, as well as scattering effects for a thinly layered medium. This solution is, however, limited to normally incident P-waves and weak contrasts of the material properties. Pride and Berryman (2003a,b) developed an effective theory for heterogeneous porous media, which are assumed composed by a mixture of two porous phases. This model does not impose any restrictions with regard to the geometry of the mesoscopic heterogeneities except that the overall mechanical and hydraulic responses of the composite medium remain isotropic. Despite its fundamentally generic nature, this approach only provides closed solutions for a restricted number of idealized model geometries. Recently, Kudarova et al. (2013) proposed an effective poroelastic model for periodically layered media, which allows for arbitrary contrasts of the material properties but is limited to normally incident P-waves.
The objective of this work is to propose a generalized poroelastic model for layered media, which is capable of accounting for both Biot's global and interlayer flows as well as for the anisotropy induced by the mesoscopic layering. This model therefore naturally predicts the existence of fast and slow P-waves and of SV-and SH-waves. We quantitatively compare the predicted phase velocity and the attenuation of P-wave and S-wave propagating normally to the layering and the corresponding properties obtained through numerical simulations of wave propagation. For this direction of propagation, the phase velocity and attenuation of SV-and SH-waves are equivalent. Such comparisons, which are technically challenging and computationally expensive and correspondingly scarce in the literature, represent an essential step for validating theoretical developments of this kind. Indeed, the reason for limiting the comparison of the proposed model with the numerical benchmark to the one-dimensional case is the very high computational cost of the latter.

E F F E C T I V E A N I S O T R O P I C P O R O E L A S T I C M O D E L
Effective models are based on the generic and idealized concept of replacing a heterogeneous medium with an equivalent homogeneous medium. In this work, we consider a periodically layered medium composed of two alternating layers. The period of this system is thus given as H = h 1 + h 2 with h 1 and h 2 denoting the layer thicknesses. Each layer corresponds to a homogeneous and isotropic poroelastic medium obeying Biot's equations. In the long-wavelength limit λ H, it is reasonable to assume that the same governing equations describe the effective model. To account for the overall anisotropy induced by the layering, Biot's equations need to be derived for a transversely isotropic poroelastic medium. The poroelastic coefficients of the effective model can then be evaluated from those of the individual layers through volume averaging. A sketch of the overall model setup and the scale relations is shown in Figure 1.

Biot's equations for effective transversely isotropic media
Assuming time dependence of the form e −iωt , Biot's equations of poroelasticy describing wave propagation in a homogeneous transversely isotropic medium can be written as (e.g., Carcione 2001) where u and w are the average solid and the relative fluid-solid displacements over the period H of the layering, which in turn defines the averaging volume of the sample. The superscript * denotes the effective poroelastic parameters that need to be expressed in terms of volume averages of the corresponding parameters prevailing in the layers. The effective bulk and pore fluid densities are given by (e.g., Pride et al. 2004) where the brackets <> denote the volume average, which is defined for a generic function f (z) as < f >= 1 /H, with f 1 and f 2 denoting the values assumed by the function in the two layers.
The bulk stress tensor σ and the pore fluid pressure p f can be written in Voigt notation as (e.g., Carcione 2001) where e is the vector of the strain components of the solid phase, and C * m is the stiffness tensor of the effective dry matrix. The coefficients α I are the anisotropic extensions of the Biot-Willis constant, which in a transversely isotropic medium are given by (e.g., Carcione 2001) where K s is the bulk modulus of the solid grains, which is assumed homogeneous in this study. The corresponding extension of the fluid storage modulusM is given by (e.g., Carcione 2001; Pride et al. 2004) where φ * =< φ > is the effective porosity, andK = 1 9 (2C * m 11 + 2C * m 12 + 4C * m 13 + C * m 33 ). We assume 1/K * f =< 1/K f > for the effective bulk modulus of the pore fluid (e.g., Gelinsky and Shapiro 1997a).
The stiffness tensor of the effective dry matrix C * m has five independent coefficients in a transversely isotropic medium. These coefficients do not depend on the fluid properties and thus can be obtained by means of Backus averaging applied to the dry frame properties of the layers (e.g., Gelinsky and Shapiro 1997a). Their expression is explicitly given in Appendix A. Analogously, the stiffness tensor of the effective saturated matrix is linked to the dry frame properties through (Gassmann 1951) This expression is inserted into equation (8) for the bulk stress and valid at frequencies that are sufficiently low for the fluid pressure to be equilibrated between the pore space of the layers. Conversely, when considering wave-induced fluid flow between the layers, the frequency-dependence saturated stiffness tensor C * U needs to be accounted for. As shown in Appendix B, we therefore replace the equation (14) with the corresponding analytical equation for C * U (ω) provided by Krzikalla and Müller (2011). Please note that when the dry frame elastic moduli of the layered system are homogeneous, the expressions for the effective dry and saturated stiffness tensors do indeed correctly reduce to their isotropic equivalents (e.g., Johnson 2001).
The remaining parameterρ * (ω) represents the effective relative flow resistance. In a transversely isotropic medium,ρ * (ω) is a tensor with two independent components (Gelinsky and Shapiro 1997a) In analogy with previous works (Gelinsky and Shapiro 1997a;Kudarova et al. 2013), it would be possible to assume the effective relative flow resistance as These coefficients, which control the frequency dependence of Biot's global flow mechanism, are widely used and heuristically justified by their consistency with corresponding hydraulic models based on Darcy's law (e.g., Bear 1988). However, in our model, we propose for the effective components where η * =< η >. These coefficients allow for a better agreement with respect to the numerical benchmark, as it will be shown later. The dynamic permeability κ d (ω) is defined by the model of Johnson et al. (1987) where the parameter n j is taken to be 8 (e.g., Pride 2005) and not critical in the course of our theoretical investigation.

Plane-wave solutions
We consider a general plane-wave solution for homogeneous waves. The displacements vectors are of the form u = u 0 e ik·x and w = w 0 e ik·x , where u 0 and w 0 are complex constant vectors and k is the wavenumber vector. Assuming k = k(l x , l y , l z ), where k is the magnitude of the wavenumber vector and l x , l y , and l z are the direction cosines, and substituting the plane-wave solution into the equations of motion (4) and (5) yields the Kelvin-Christoffel equation where and The solution of the eigenvalue problem of equation (21), whose characteristic equation is det( D −1 − v c 2 I 6 ) = 0, provides the anisotropic and complex-valued velocities v c 2 = ω 2 /k 2 for all the possible wave modes.

V A L I D A T I O N
To validate the proposed model, we compare the analytical solutions for the P-wave and S-wave phase velocities and attenuation for the particular case of wave propagation perpendicular to layering with those obtained from numerical simulations. To do so, we consider a model consisting of the periodic alternation of a harder water-saturated layer and a softer gas-saturated layer ( Table 1). The thicknesses of waterand gas-saturated layers are 4 cm and 1 cm, respectively, which in turn corresponds to an effective gas saturation of S g = 20%. Please note that the primary purpose of this setup is to validate the proposed methodology for the particularly challenging case when the critical frequencies of mesoscopic and macroscopic wave-induced fluid flows (WIFFs) are similar and the attenuation produced by the two mechanisms is of comparable magnitude. As such, the parametrization of this model is somewhat hypothetical and rather generic, which involves a heterogeneous distribution of mechanical and hydraulic properties as well as two different saturating pore fluids. This parameterization does not intend to reproduce any particular geological scenario.

P-and S-wave solutions for vertical direction of propagation
To find the solution for normally incident P-and S-waves, we solve the eigenvalue problem of equation (21) for l x = l y = 0 and l z = 1. The three independent eigenvalues v c 2 we obtain correspond to an S-wave and to the fast and slow P-waves. Please note that, for this direction of propagation, there exists only one S-wave since the SV-and the SH-wave modes are equivalent. The complex velocities for P-waves are given by where the plus and minus signs refer to the slow and fast propagating modes, respectively, and The complex velocity for the S-wave is given by The P-and S-wave phase velocities and inverse quality factors can then be evaluated from the corresponding complex velocities

Numerical modelling
Simulations of P-and S-waves propagating perpendicular to the layering are based on a numerical solution of Biot's onedimensional equations in the space-frequency domain. We employ a finite-element technique with standard linear basis functions in conjunction with absorbing boundary conditions at the edges of the numerical domain (Appendix C). The overall size of the considered computational domain is 100 m, which is discretized using elements of size z = 5 · 10 −4 m. The source has the time history of a Ricker wavelet and is located at z 0 = 45 m. We perform several numerical simulations for a wide range of values for the source central frequency f 0 , which then allows for evaluating phase velocities and attenuation of P-wave and S-wave over a broad frequency range. In order to estimate these properties, we apply the crossspectrum and the spectral ratio methods (e.g., Molyneux and Schmitt 2000;Baron and Holliger 2011;Milani et al. 2015) to the synthetic seismograms recorded for all pairwise combinations of receivers, which are located at z r = 50, 52, 54, and 56 m in the computational domain. The primary reason why we limit the validation procedure to normal incidence is that, in conjunction with a very large domain, we need an excessively fine discretization of our numerical model to uniformly achieve the required accuracy over a very broad frequency range covering approximately three orders of magnitude.

Biot's global flow
In the following analyses, we describe the situations in terms of the transition frequencies f Biot and f meso . This allows us to illustrate how these two frequencies are affected by changes in permeability. Before considering the scenario depicted by the material properties of Table 1, we validate the proposed analytical solution in the exclusive presence of Biot's global flow. To realize this scenario, we increase the permeabilities of both layers by two orders of magnitude compared with the values given in Table 1, which shifts the attenuation peaks associated with Biot's global flow and interlayer flow towards significantly lower and higher frequencies, respectively. Figure 2 shows a comparison of the P-wave phase velocity and the inverse quality factor predicted by the method proposed in this study, the model of White et al. (1975), and numerical simulations. There is good agreement between the P-wave phase velocity and inverse quality factor predicted by the proposed effective poroelastic model and those obtained from numerical simulations over a broad frequency range. Conversely, the solution provided by the model of White et al. (1975) only provides a good agreement for the P-wave phase velocity in the low-frequency limit where Gassmann's (1951) relation C * U 33 = C * m 33 + α 3 2M is satisfied by both models. The mismatch over the remainder of the considered frequency range is expected since the model of White et al. (1975) is not designed for describing Biot global flow, which, in this case, occurs at f Biot 100 Hz. Significant discrepancies between the proposed effective poroelastic model and the numerical simulations are then observed at frequencies higher than ∼ 5 kHz. Above this threshold frequency, the long-wavelength approximation, upon which the effective model is based, becomes inappropriate. For visualization purposes, we therefore consider numerical solutions at frequencies lower than 10 kHz. Similarly, Figure 3 shows good agreement between the S-wave phase velocity and the inverse quality factor predicted by the proposed effective poroelastic model and those obtained from numerical simulations over a broad frequency range. Moreover, in this case, the low-frequency limit for the S-wave phase velocity is common to both models and analytically described by Gassmann's (1951) relation C * U 44 = C * m 44 . Figures 2 and 3 show that using equations (16) and (17) leads to an identical solution for the P-wave (dashed blue curve), whereas for the S-wave, it produces a solution  Figure 1 and wave propagation perpendicular to the layering. Model parameters correspond to those of Table 1 apart from the permeability, which has been increased by two orders of magnitude. Green and red curves correspond to the analytical solution of White et al. (1975) and the effective poroelastic model proposed in this study, respectively. The dashed blue curve corresponds to this effective poroelastic model using the averaging rule of equation (17) for the relative flow resistance. The black curve represents the results obtained from numerical simulations.  Figure 1 and wave propagation perpendicular to the layering. Model parameters correspond to those of Table 1 apart from the permeability, which has been increased by two orders of magnitude. Green and red curves correspond to the analytical solution of White et al. (1975) and the effective poroelastic model proposed in this study, respectively. The dashed blue curve corresponds to this effective poroelastic model using the averaging rule of equation (16) for the relative flow resistance. The black curve represents the results obtained from numerical simulations. that significantly deviates from the numerical simulations. Conversely, the coefficients of equations (18) and (19) proposed in this work allow for a reconciliation of the S-wave phase velocity and inverse quality factor with respect to the numerical simulations. Furthermore, the proposed coefficients are proportional to the arithmetic and harmonic averages of the permeability for fluid-solid displacements tangential or normal to the layering, respectively, thus naturally retaining  Figure 1 and wave propagation perpendicular to the layering. Model parameters correspond to those of Table 1 apart from the permeability, which has been lowered by two orders of magnitude. Green and red curves correspond to the analytical solutions of White et al. (1975) and the effective poroelastic model proposed in this study, respectively. The black curve denotes the results obtained from numerical simulations. a generic characteristic of fluid flow in layered media. In this context, it is also important to note that the coefficients that we propose reduce to those of equations (16) and (17) for periodically layered media saturated by a single fluid.

Interlayer flow
In analogy to the case of Biot's global flow considered earlier, we can validate our model for mesoscopic WIFF by reducing the permeabilities of both layers by two orders of magnitude with respect to the values given in Table 1. This results in a clear separation of the interlayer and global flows by shifting the corresponding attenuation peaks towards significantly lower and higher frequencies, respectively. For the interlayer flow case, we only consider the P-wave solution because the S-wave does not produce WIFF at normal incidence to layering. Figure 4 shows good agreement over a broad frequency range between the P-wave phase velocity and inverse quality factor predicted by the proposed effective poroelastic model (red curve) and the model of White et al. (1975) (green curve), as well as with corresponding values obtained from the numerical simulations (black curves) over a broad frequency range. As expected, when f meso 200 Hz is much smaller than f Biot 1 MHz, the solution by White et al. (1975) provides a proper description of the P-wave attenuation and velocity dispersion due to interlayer WIFF. Compared with the previously considered scenario where only Biot's global flow prevailed (Figure 2), scattering effects become visible at a higher threshold frequency of ∼ 10 kHz. The reason for this shift is that scattering effects now prevails in the highfrequency limit of the mesoscopic WIFF where the stiffness contrast between the layers is reduced (e.g., Johnson 2001).

Biot's global and interlayer flows
Finally, we consider the combined effects of the two WIFF mechanisms acting at similar frequencies by using the permeability values given in Table 1. Figure 5 shows that, for frequencies smaller than ∼ 4 kHz, the proposed effective poroelastic model agrees well with the numerical results, whereas the model of White et al. (1975) systematically underestimates the attenuation at all considered frequencies. In the frequency range between 4 kHz and 10 kHz, some discrepancies between the proposed model and the numerical solutions can be noticed. These discrepancies are likely be related to the fact that our analytical model is based on a superposition of the global and interlayer flow mechanisms and hence does not account for inertial effects arising close to Biot's critical frequency in the mesoscopic flow mechanism. However, in this  Figure 1 and wave propagation perpendicular to the layering. Model parameters correspond to those of Table 1. Green and red curves correspond to the analytical solutions of White et al. (1975) and the effective poroelastic model proposed in this study, respectively. The black curve denotes the results obtained from numerical simulations.
context, it is important to note that the long-wavelength assumption, upon which our model is based, looses its validity in this frequency range, and hence scattering effects are also likely to contribute to the observed discrepancies with regard to the numerical solutions.

A N G L E D E P E N D E N C Y O F P -, S V -, A N D S H -W A V E S
In order to compute the analytical solutions valid for any angle of propagation, we consider wave propagation in the xz plane and evaluate the eigenvalues of equation (21) by setting l y = 0. The resulting linear system is associated with a 6 × 6 matrix that can be decomposed into two uncoupled linear subsystems of dimensions 4 × 4 and 2 × 2. This mathematical decomposition reflects the fact that the motion of the SHwave is decoupled from the other wave modes. The resulting solutions for the phase velocities and inverse quality factors of P-, SV-, and SH-waves depend on the angular frequency ω and on the angle of incidence. This is illustrated in Figure 6 where, for a frequency of 4 kHz, we compare the angle-dependent solutions of the proposed model with those of the model of Krzikalla and Müller (2011) for the scenario accounting for both Biot's global and interlayer flows (Table 1). At normal incidence, the agreement between the proposed model and the numerical benchmark is excellent at this frequency ( Figure 5).
As expected, we see that, for both models, the phase velocities of the three wave modes exhibit a pronounced anisotropy due to the elastic contrast between the two layers for both models (e.g., Carcione 2001). That is, the P-and SH-wave velocities are increasing as the direction of propagation tends to be parallel to the layering, whereas the SV-wave velocity reaches its maximum for an angle of propagation of 45 • . Accounting for both Biot's global and interlayer wave-induced fluid flow (WIFF) increases the phase velocity with respect to the model of Krzikalla and Müller (2011), which only considers interlayer WIFF, at all angles. This discrepancy is expected to vanish in the low-frequency limit and to reach its maximum in the high-frequency limit ( Figure 5). Again, the attenuation predicted by the proposed model is significantly higher than that of Krzikalla and Müller (2011). The reason for the differing P-and S-wave attenuation characteristics is that the former is affected by Biot's global and interlayer WIFF, whereas the latter is entirely governed by Biot's global flow. This also explains the observation that the S-wave attenuation level predicted by mesoscopic model of Krzikalla and Müller (2011) is largely negligible, whereas it is very significant for the model proposed in this study. The moderate anisotropy of the inverse P-wave quality factor for the effective poroelastic model is entirely caused by mesoscopic WIFF, whereas the Biot mechanism is purely isotropic as for the SV-and SH-waves attenuation ( Figure 6). Figure 6 P-, SV-, and SH-wave phase velocities and inverse quality factors at frequency of 4 kHz as functions of the angle of incidence predicted by the proposed effective poroelastic model (red curves) and by the model of Krzikalla and Müller (2011) (green curves). The setup of the model and its parametrization are given in Figure 1 and Table 1, respectively. Therefore, both Biot's global and interlayer flows are considered in this situation ( Figure 5).
The main cause for the isotropic attenuation related to Biot's global flow is that only the water-saturated layer contributes significantly to the corresponding WIFF; therefore, the two coefficients for the relative flow resistance of equations (18) and (19) are very similar. For completeness, we have also verified that the two solutions for the SV-and SH-waves are equivalent for propagation perpendicular to the layering.

D I S C U S S I O N A N D C O N C L U S I O N S
We have developed a generalized effective poroelastic model which superimposes Biot's global and interlayer flows for a medium composed of periodically distributed mesoscopic layers, which also accounts for the inherent anisotropy related to the layering. The proposed model is useful whenever at least one of the layers is characterized by high permeability, thus shifting the characteristic frequency of Biot's macroscopic mechanism into the sonic to seismic frequency range, including all cases producing a negligible amount of attenuation due to interlayer WIFF. To include interlayer flow, we replace the effective stiffness tensor for the undrained sample in the low-frequency limit with a corresponding frequencydependent tensor, which analytically describes WIFF in the mesoscopic scale range as well as the anisotropy associated with the layering. The general solution resulting from this model correctly predicts the existence of the fast and slow P-wave modes as well as of quasi and pure S-wave modes. The proposed model was validated quantitatively through comparisons with a numerical solution of Biot's poroelastic equations of motion for wave propagation perpendicular to the layering and qualitatively by assessing its angle dependence and comparing it with the mesoscopic model of Krzikalla and Müller (2011). The reason for the limitation of the quantitative validation to the one-dimensional case is the high computational cost of the numerical solution, due to the very fine spatial discretization that is required in conjunction with the large model size and broad frequency range of approximately three orders of magnitude.
Comparisons of the P-and S-wave phase velocities and inverse quality factors for normal incidence with the corresponding numerical results show very good agreement as long as the considered frequencies fulfill the inherent longwavelength assumption with respect to the layer thicknesses. When interlayer flow occurs at frequencies much lower than those governing Biot's global flow, the generalized model proposed in this study is also in good agreement with the classic mesoscopic model of White et al. (1975). For the scenarios f Biot f meso and f meso f Biot , we were able to quantitatively corroborate this validation only for the WIFF mechanism prevailing at lower frequencies because, at higher frequencies, scattering effects were increasingly contaminating the numerical solution (Figs. 2 and 4). Attempts to extend this validation procedure to the respective higher frequent WIFF mechanism were hampered by the virtual absence of geologically realistic models satisfying f Biot f meso f scattering or f meso f Biot f scattering . A reasonably strong indication for the fundamental validity of the proposed model is, however, provided by the fact that we observe good agreement with its numerical benchmark for f meso f Biot at frequencies close to 10 kHz , where the model by White et al. (1975) diverges in a systematic and significant manner (Figure 4).
A limitation of the proposed model is that, for f Biot f meso , it simply superimposes Biot's global and interlayer WIFF without being able to account for inertial effects on interlayer flow. The reason for this is that the analytical solutions of White et al. (1975) and Krzikalla and Müller (2011) for mesoscopic WIFF, which are employed in our model, are exact only for frequencies lower than f Biot and only provide approximate solutions for higher frequencies. The magnitude of the associated inaccuracies is as of yet unknown and inherently difficult to quantify through comparisons with numerical solutions as scattering effects become important in the same frequency range. In the view of the fact that the corresponding error is caused by the influence of the unaccounted dynamic terms associated with the global flow mechanism on WIFF in the mesoscopic scale range, it is reasonable to assume that its magnitude is likely to be related to that of the corresponding attenuation. In this context, it is important to note that, despite the rather high attenuation level, the proposed model provides a solution that is not too far from its numerical benchmark in the frequency range where both inertial and scattering effects are expected to be at play ( Figure 5).
Finally, it is important to note that the proposed methodology can be further generalized to account for the effects of squirt flow. To this end, one needs to replace the current dry frame elastic moduli by those resulting from the considered squirt flow model (e.g., Gurevich et al. 2010). A corresponding approach has, for example, been taken by Milani et al. (2015) when investigating the potential importance of this mechanism in unconsolidated sediments. Given the methodological character of this work, we chose to ignore this mechanism as its inclusion would make a rigorous numerical benchmarking impossible.

A C K N O W L E D G E M E N T S
This work has been supported by a grant from the Swiss National Science Foundation. L.B. Monachesi gratefully acknowledges an extended visit to the University of Lausanne financed by the Fondation Herbette. L.B. Monachesi and J.I. Sabbione would like to thank Dr. Claudia Ravazzoli and Dr. Fabio Zyserman for their comments on the implementation of the numerical modelling code.

A P P E N D I X A : S T I F F N E S S T E N S O R O F T H E E F F E C T I V E D R Y M A T R I X
The elastic stiffness tensor for the dry matrix of a layered medium in the long-wavelength limit is effectively transversely isotropic and given by the so-called Backus average (Backus 1962). For the symmetry axis in the z-direction, i.e., horizontal layering, the stiffness tensor C * m can be written in Voigt notation: Effective anisotropic poroelastic model 1147 with P d = λ m + 2μ and λ m denoting the dry frame P-wave modulus and Lamé coefficient.

A P P E N D I X B : S T I F F N E S S T E N S O R O F T H E E F F E C T I V E S A T U R A T E D M A T E R I A L
The poroelastic extension of the Backus average (Backus 1962) provides the relaxed and unrelaxed limits for fluidsaturated porous layered media (Gelinsky and Shapiro 1997a). At intermediate frequencies, interlayer flow, which causes velocity dispersion and attenuation, may occur due to stiffness contrasts between the layers. The complex-valued and frequency-dependent effective undrained stiffness tensor C * U (ω) can be written in Voigt notation (Krzikalla and Müller 2011) (ω) C * U 12 (ω) C * is the scalar complex-valued relaxation function. White et al. (1975) derived the equivalent P-wave modulus C 33 (ω) based on Biot's (1956aBiot's ( , 1956b) theory for periodically layered media with two alternating porous layers where the subscripts refer to the properties of the two layers.
In this equation, α is the so-called Biot-Willis constant, P U = P d + α 2 M is the undrained P-wave modulus, and M is the fluid storage modulus. These parameters are given by where K m denotes the bulk modulus of the dry matrix. The remaining parameter k s = √ iω/D is the wavenumber of Biot's slow wave. The pressure diffusivity D is given by The set of equations for the relaxed moduli C relaxed I J is given by (Gelinsky and Shapiro 1997a)   where λ U = λ d + α 2 M is the undrained Lamé coefficient.