Spring-damper equivalents of the fractional, poroelastic, and poroviscoelastic models for elastography

TL;DR

This paper demonstrates that complex fractional, poroelastic, and poroviscoelastic models in elastography can be represented using simple spring-damper systems, simplifying their interpretation and parameter estimation.

Contribution

It shows that various advanced tissue models can be expressed as elementary viscoelastic models, reducing parameter complexity and enabling more physical interpretation.

Findings

Poroelastic shear wave modeled by a three-component Zener model.

Extended poroviscoelastic model equivalent to a four-parameter model.

Complex models can be interpreted through elementary viscoelastic components.

Abstract

In MR elastography it is common to use an elastic model for the tissue's response in order to properly interpret the results. More complex models such as viscoelastic, fractional viscoelastic, poroelastic, or poroviscoelastic ones are also used. These models appear at first sight to be very different, but here it is shown that they all may be expressed in terms of elementary viscoelastic models. For a medium expressed with fractional models, many elementary spring-damper combinations are added, each of them weighted according to a long-tailed distribution, hinting at a fractional distribution of time constants or relaxation frequencies. This may open up for a more physical interpretation of the fractional models. The shear wave component of the poroelastic model is shown to be modeled exactly by a three-component Zener model. The extended poroviscoelastic model is found to be equivalent to what is called a non-standard four-parameter model. Accordingly, the large number of parameters in the porous models can be reduced to the same number as in their viscoelastic equivalents. As long as the individual displacements from the solid and fluid parts cannot be measured individually the main use of the poro(visco)elastic models is therefore as a physics based method for determining parameters in a viscoelastic model.

Tables

Table 1. Table 1: Biot-Stoll parameters and their effect on the shear wave. The pore radius parameter is not included as the flow in the pores is assumed to be laminar. The right-hand column shows whether the parameter affects the shear wave or not

Biot model parameters			Shear?
Bulk properties:
$\varphi$	[0 … 1]	Porosity	Y
${\rho }_s$	[kg/m3]	Solid density	Y
${\rho }_f$	[kg/m3]	Fluid density	Y
$K_s$	[Pa]	Bulk modulus, solid	N
$K_f$	[Pa]	Bulk modulus, fluid	N
Fluid parameters:
$\eta$	[Pa s]	Viscosity	Y
$B$	[m2]	Permeability	Y
$\alpha$	[1 … 3]	Tortuosity	Y
Rigid frame response parameters:
$K_r$	[Pa]	Bulk modulus	N
${\mu }_r$	[Pa]	Shear modulus	Y

Equations

$$\sigma (t)=G(t)\ast \frac{\partial ϵ(t)}{\partial t}.$$

$$\left[1+\sum _{k=1}^p{a}_k\frac{{\partial }^k}{\partial t^k}\right]\sigma (t)=\left[E_e+\sum _{k=1}^q{b}_k\frac{{\partial }^k}{\partial t^k}\right]ϵ(t).$$

$$G(t)=G_e+G_{\tau }+G_{-}\delta (t),G_{\tau }(t)=\sum _{n=0}^{N-1}{E}_n\exp (-t/{\tau }_n).$$

$$E(\omega )=\frac{\sigma (\omega )}{ϵ(\omega )}=\frac{E_e+{\sum }_{k=1}^q{b}_k(i\omega {)}^k}{1+{\sum }_{k=1}^p{a}_k(i\omega {)}^k}$$

$$\sigma (t)=E_eϵ(t)+\eta \frac{\partial ϵ(t)}{\partial t},$$

$$G(t)=E_e+\eta \delta (t).$$

$$E(\omega )=E_e+i\omega \eta =E_e(1+i\omega {\tau }_{ϵ})$$

$$\sigma (t)+a_1\frac{\partial \sigma (t)}{\partial t}=E_eϵ(t)+\eta \frac{\partial ϵ(t)}{\partial t},$$

$$G(t)=E_e+E_e(\frac{{\tau }_{ϵ}}{{\tau }_{\sigma }}-1)e^{-t/{\tau }_{\sigma }}.$$

$${\tau }_{\sigma }=\eta /E\le {\tau }_{ϵ}=\eta /E^{\prime },\frac{1}{E^{\prime }}=\frac{1}{E_e}+\frac{1}{E}.$$

$$E(\omega )=E_e\frac{1+i\omega {\tau }_{ϵ}}{1+i\omega {\tau }_{\sigma }}$$

$${\left(\frac{k}{\omega }\right)}^2=\rho \kappa (\omega )=\frac{\rho }{E(\omega )}$$

$$\kappa (\omega )=\frac{1}{E(\omega )}=\frac{1}{E_e}\frac{1+{\omega }^2{\tau }_{ϵ}{\tau }_{\sigma }-i\omega ({\tau }_{ϵ}-{\tau }_{\sigma })}{1+{\omega }^2{\tau }_{ϵ}^2}$$

$$G_{\tau }(t)=b{\int }_0^{\infty }{R}_{\sigma }(\tau )e^{-t/\tau }d\tau ,$$

$$S_{\sigma }(\Omega )=b\frac{R_{\sigma }(1/\Omega )}{{\Omega }^2}$$

$$G_{\tau }(t)={\int }_0^{\infty }{S}_{\sigma }(\Omega )e^{-\Omega t}d\Omega .$$

$$J_{\tau }(t)=a{\int }_0^{\infty }{R}_{ϵ}(\tau )(1-e^{-t/\tau })d\tau$$

$$S_{ϵ}(\Omega )=a\frac{R_{ϵ}(1/\Omega )}{{\Omega }^2}$$

$$J(t)=-{\int }_0^{\infty }{R}_{\sigma }(ϵ)(1-e^{-t/\tau })d\tau$$

$$\sigma (t)+{\tau }_{\sigma }^{\alpha }\frac{{\partial }^{\alpha }\sigma (t)}{\partial t^{\beta }}=E\left[ϵ(t)+{\tau }_{ϵ}^{\alpha }\frac{{\partial }^{\alpha }ϵ(t)}{\partial t^{\alpha }}\right]$$

$$E(\omega )=E_e\frac{1+(i\omega {\tau }_{ϵ}{)}^{\alpha }}{1+(i\omega {\tau }_{\sigma }{)}^{\alpha }}$$

$$\sigma (t)=E\left[ϵ(t)+{\tau }^{\alpha }\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}ϵ(t)\right].$$

$$E(\omega )=E_e\left(1+(i\omega {\tau }_{ϵ}{)}^{\alpha }\right)$$

$$R_{ϵ}(\tau )=\frac{1}{\pi \tau }\frac{\sin \alpha \pi }{(\tau /{\tau }_{ϵ}{)}^{\alpha }+(\tau /{\tau }_{ϵ}{)}^{-\alpha }+2\cos \alpha \pi }$$

$$R_{ϵ}(\tau )\propto \{\begin{array}{ll}{\tau }^{\alpha -1}& \text{for}\tau /{\tau }_{ϵ}\ll 1,\\ {\tau }^{-\alpha -1}& \text{for}\tau /{\tau }_{ϵ}\gg 1\end{array}$$

$$S_{ϵ}(\Omega )=\frac{1}{\pi E}\frac{({\tau }_{ϵ}^{\alpha }-{\tau }_{\sigma }^{\alpha }){\Omega }^{\alpha -1}\sin \alpha \pi }{(\Omega {\tau }_{ϵ}{)}^{2\alpha }+1+2(\Omega {\tau }_{ϵ}{)}^{\alpha }\cos \alpha \pi }.$$

$$S_{ϵ}(\Omega )\propto \{\begin{array}{ll}{\Omega }^{\alpha -1}& \text{for}\Omega {\tau }_{ϵ}\ll 1,\\ {\Omega }^{-\alpha -1}& \text{for}\Omega {\tau }_{ϵ}\gg 1\end{array}$$

$$E(\omega )=K+iH$$

$${\left(\frac{k}{\omega }\right)}^2=\rho ({\kappa }_r-i{\kappa }_i),$$

$${\kappa }_r=\frac{1}{\mu }\frac{1+{\gamma }_{22}\frac{{\gamma }_{11}{\gamma }_{22}-{\gamma }_{12}^2}{({\gamma }_{12}+{\gamma }_{22}{)}^2}{\left(\frac{f}{f_c}\right)}^2}{1+{\left(\frac{{\gamma }_{22}}{{\gamma }_{12}+{\gamma }_{22}}\right)}^2{\left(\frac{f}{f_c}\right)}^2},{\kappa }_i=\frac{1}{\mu }\frac{f}{f_c}\frac{{\gamma }_{12}+{\gamma }_{22}}{1+{\left(\frac{{\gamma }_{22}}{{\gamma }_{12}+{\gamma }_{22}}\right)}^2{\left(\frac{f}{f_c}\right)}^2}$$

Full text

For a medium expressed with fractional models, many elementary spring–damper combinations are added, each of them weighted according to a long-tailed distribution, hinting at a fractional distribution of time constants or relaxation frequencies. This may open up for a more physical interpretation of the fractional models.

The shear wave component of the poroelastic model is shown to be modeled exactly by a three-component Zener model. The extended poroviscoelastic model is found to be equivalent to what is called a non-standard four-parameter model. Accordingly, the large number of parameters in the porous models can be reduced to the same number as in their viscoelastic equivalents. As long as the individual displacements from the solid and fluid parts cannot be measured individually the main use of the poro(visco)elastic models is therefore as a physics based method for determining parameters in a viscoelastic model.

1 Introduction

In modeling of data from MR elastography, it is common to use a simple elastic model for the medium. This is the case even for ultrasound elastography. For more accurate modeling, there are three families of models that are used. These are the linear viscoelastic models such as the Kelvin-Voigt and Zener models, the fractional extensions of these models, and poroelastic models based on the theory of Biot.

The linear viscoelastic models are among those that have been used for fitting frequency dependency of shear wave data from MR elastography in the brain [1]. The fractional Kelvin-Voigt model was fitted to breast MR elastography data in [2] and also analyzed in [3] and compared to other models for elastography in [4].

It may be argued that these single-phase models are too simplistic and that in tissue a bi-phasic model which distinguishes between the solid and liquid phases would be more accurate. This potential has already been demonstrated in models of the quasi-static biomechanics of hydrocephalus [5] and of infusion-induced swelling in the brain [6]. The poroelastic model has also been used for elastography. In [7] ultrasound elastography was modeled with either an elastic model (i.e. without viscosity) or a simplified poroelastic model which depended on porosity, composite density and fluid density. Similarly [8] evaluated the full poroelastic model vs. an elastic one for MR elastography. In [9] they went one step further and compared the poroleastic model with a viscoelastic one with a complex shear modulus. They also tried to use these models for inversion and observed that the viscoelastic model produced better reconstructions at 50 Hz while the poroelastic model was superior at 1 Hz.

The challenge in making a reconstruction algorithm based on the poroelastic model is the need for capturing the displacement of the solid and the fluid independently. Chapter 5 in [10] states that due to the voxel size of MR elastography, it cannot detect properties of individual pores, but only sees a homogeneous effective medium. Further, since it is sensitive to signals from hydrogen atoms in the voxel, one cannot separate the signal from the solid and the fluid, even if the solid should contain fluid which is considered not to be in the pores. Likewise ultrasound is scattered from the soft tissue which is mainly part of the solid matrix. Ultrasound will therefore produce a strain estimate in the solid matrix which is only indirectly influenced by the wave in the fluid [7]. In this paper, the goal is however not to develop a poroelastic model for reconstruction of the MR elastography image, but only to discuss it in the framework of explaining variations in parameters due to physiological and physical parameters.

The claim of this paper is that these models are much more similar than they appear to be. The fractional models can be developed as sums of ordinary viscoelastic elements weighted in a particular way. The fractional model has its strength in that it gives a parsimonious description of the phenomenon, i.e. one with a minimal number of parameters, especially when power-law variation in frequency is observed. But it is not fundamentally different. Likewise the poroelastic model for the wave mode of interest for elastography, the shear wave, can be described in terms of standard viscoelastic models. In this case it is the viscoelastic model which requires the smallest number of parameters. The strength of the poroelastic formulation is that it gives a way of finding how these parameters depend on physical and even physiological parameters.

The linear viscoelastic model is therefore first generalized in Sec. 2 from the simple two- and three-component models to chains of spring-damper elements described by time- and frequency-spectral functions. These functions are used to show that the fractional viscoelastic models in Sec. 3 can be described as sums of ordinary viscoelastic models. A surprising result is that the weighting in the sum follows a long-tailed distribution reminiscent of a fractal.

The shear wave solution of the poroelastic theory is then developed in Sec. 4 both from the original formulation of Biot [11] and also from that of Stoll [12]. Somewhat unexpectedly, it is found that there is a one to one correspondence between the poroelastic shear wave response and that of a simple spring-dashpot network.

2 Linear viscoelastic models

The linear viscoelastic model is expressed in three different ways. The different ways are needed in order to generalize from the simple two and three-component models (e.g. Kelvin-Voigt and Zener) to higher order models.

In order to illustrate the similarities between models it is sufficient here to express the models in one dimension, although three dimensions are really needed for a complete shear wave description. This section, as well as Sec. 3 builds to a large degree on [13].

2.1 Three descriptions

In linear viscoelasticity there are three different ways of describing the medium’s response. The first is the hereditary model of Boltzmann [14, 15] where the constitutive equation is a convolution integral [13] (ch. 2):

[TABLE]

The kernel $G(t)$ is called the relaxation modulus and represents a fading memory, i.e. one where changes in the past have less effect now than more recent changes. In order to ensure causality, the kernel, $G(t)$, has to be zero for non-negative time. The relaxation modulus is the strain response of the system to a step excitation in strain.

The second description is a linear differential equation between stress and strain with constant coefficients:

[TABLE]

The third description is in the form of a relaxation spectrum or Prony expansion:

[TABLE]

A physically realizable model is obtained if $G(t)$ is modeled by a parallel network of pairs of springs and dashpots in series where the spring constants are $E_{n}$ and the viscosity of the dashpots are $\eta_{n}=\tau_{n}E_{n}$. Furthermore, both the spring constants and the viscosities are non-negative. This is the Maxwell-Wiechert model of Fig. 1. The left-hand spring leads to the constant $G_{e}=E_{e}$, the equilibrium modulus, and if there is a dashpot directly across the terminals (e.g. if $E_{1}=0$), the response will have an impulse at time zero as well given by $G_{-}$. This is the case in the Kelvin-Voigt model which will soon be discussed.

Not all the coefficient sets of (2) will lead to a fading memory model, so the linear differential equation is the most general model. Also, a model described as a relaxation spectrum, (3), can always be described with a linear differential equation, but not vice versa. In addition, a relaxation spectrum model, (3), always results in a fading memory model as it consists of sums of positively weighted falling exponentials.

The frequency domain equivalent of (2) is also useful:

[TABLE]

The dynamic modulus, $E(\omega)$, is often called $G(\omega)$ in elastography, but since $G$ here means the relaxation modulus in this paper, $E$ is used instead.

2.2 Elementary linear viscoelastic models

The two simplest viscoelastic models are the Kelvin-Voigt and the Zener models as shown in Fig. 2. The Kelvin-Voigt model is found by just keeping $E_{e}$ and $\eta_{1}$ in Fig. 1. The linear differential equation is

[TABLE]

and the relaxation modulus is:

[TABLE]

Thus $p=0$, $q=1$, $a_{1}=0$, and $b_{1}=\eta$ in (2) and (4). Further $G_{e}=E_{e}$, $G_{-}=\eta$, and $G_{\tau}(t)=0$ in (3). The dynamic modulus is:

[TABLE]

where $\tau_{\epsilon}=\eta/E_{e}$. The even simpler elastic model, which is often the reference in elastography as noted in the Introduction, is obtained by setting the viscosity, $\eta$, to 0.

The Zener model adds one more term to the linear differential equation of the Kelvin-Voigt model:

[TABLE]

and the relaxation modulus is:

[TABLE]

Compared to the Kelvin-Voigt model, the effect of the second spring is to ‘soften’ the impulse in (6) into a falling exponential. In this model $p=q=1$ and $b_{1}=\eta$ in (2) and (4). The time constants are $\tau_{\epsilon}$ as in the Kelvin-Voigt model and $\tau_{\sigma}=a_{1}$. The parameters of the model relate to the physical parameters of Fig. 1 as follows:

[TABLE]

For this model the dynamic modulus is:

[TABLE]

The addition of the extra spring leads to frequency dependent denominator in the dynamic modulus which gives more degrees of freedom in fitting this model to data.

The dynamic modulus can also be included in a dispersion relation when propagating waves are involved. In [16], Eqs. (13-15), it is shown that the complex wave number in that case is:

[TABLE]

where $\rho$ is the density and $\kappa(\omega)$ is the dynamic compressibility, the inverse of the dynamic modulus. This result is important in the subsequent analysis of the poroelastic model. Note that the roles of $\tau_{\epsilon}$ and $\tau_{\sigma}$ are reversed here compared to [16]. Here the convention of [13] is followed instead.

A slightly rewritten version of the compressibility based on factoring of (11) will be needed when analyzing the poroelastic model:

[TABLE]

Judging from Fig. 1, a Zener model is characterized by two springs and a damper. These three parameters may be expressed in different ways, and one common way is via the low frequency asymptote of the propagation speed, $c_{0}$, and the two cross-over frequencies, $\omega_{\epsilon}$ and $\omega_{\sigma}$. The two frequencies express approximately where the phase velocity starts rising and where it approaches its asymptotic value $c_{\infty}=c_{0}\sqrt{\omega_{\sigma}/\omega_{\epsilon}}$ [16]. Taking into consideration that $c_{0}^{2}=E_{e}/\rho$ and thus depends on both the shear modulus and the density, then a Zener medium will depend on four independent parameters: $E_{e}$, $\rho$, $\omega_{\epsilon}$ and $\omega_{\sigma}$.

2.3 Time- and frequency-spectral functions

The continuous generalization of (3) is given in [17, 18, 13]:

[TABLE]

where $b=G_{\tau}(0)$ is a non-negative constant which for the Zener model is $b=E(\tau_{\epsilon}/\tau_{\sigma}-1)$ [18]. $R_{\sigma}(\tau)$ is a non-negative relaxation spectrum.

There is a corresponding frequency-spectral function which is

[TABLE]

By substituting $\Omega=1/\tau$ so $d\Omega=-d\tau/\tau^{2}$ this gives

[TABLE]

The relaxation modulus is, as noted, the stress response of the system to a step excitation in strain. The relaxation spectrum could also have been expressed with the creep response, which is the strain response to a step excitation in strain. It is denoted by $J(t)$ and plays the same role as $G(t)$ in (1) if $\sigma(t)$ and $\epsilon(t)$ are interchanged. The time-spectral function is:

[TABLE]

where $a=J_{\tau}(\infty)$ is a non-negative constant which for the Zener model is $a=(1/E)(1-\tau_{\sigma}/\tau_{\epsilon})$ [18]. This corresponds to a decomposition in a sum of elementary models in series as shown in Fig. 3. This is the conjugate of the Maxwell-Wiechert model of Fig. 1, i.e. a different configuration of springs and dash-pots which has the same characteristics [19]. It follows that the conjugate of the Zener model in the right-hand part of Fig. 2 is the top three elements of Fig. 3 also.

The frequency-spectral function for the creep is found by a transformation which is analogous to that for the relaxation:

[TABLE]

and the frequency-spectral function is:

[TABLE]

This frequency spectral function will be used to find the spring-damper equivalent of the fractional models.

3 Fractional models

Both the Kelvin-Voigt model and the Zener model can be generalized by introducing non-integer, fractional derivatives of order $\alpha$ in the constitutive equation (2). For the Zener model this is:

[TABLE]

The dynamic modulus is now:

[TABLE]

while in the simpler Kelvin-Voigt model shown in Fig. 4, the constant $\tau_{\sigma}=0$:

[TABLE]

The dynamic modulus is:

[TABLE]

In the limit as the frequency and/or viscosity is very large the dynamic modulus approaches $E(\omega)\rightarrow(i\omega\eta)^{\alpha}$. This is equivalent to the case where the spring of Fig. 4 can be neglected. That seems to be the case often for elastography data [20, 21].

The relaxation modulus of the fractional Kelvin-Voigt model is expressed by a power law function of time, while for the fractional Zener model it follows a Mittag-Leffler function. That function is a generalization of the exponential function with power-like behavior. The creep compliance of both models also follows a Mittag-Leffler function. That means that the creep time- and -frequency-spectral functions of the two models will be similar also.

3.1 Fractal time- and frequency-spectral function

In [17, 18, 13] it has been shown that for the fractional Zener model, the creep time-spectral function of (17) is:

[TABLE]

It was plotted in [18, 13] with linear axes and is a decreasing function of $\tau$ for small $\alpha$ and gets a more and more pronounced peak as $\alpha$ approaches 1. For the non-fractional case, $\alpha=1$ it is a delta function at $\tau/\tau_{\epsilon}=1$ showing that it is equivalent to a single relaxation process in that case.

Here the function is plotted on a log log scale in Fig. 5. That brings out the properties of this particular function in a different way than if it is plotted in a linear plot. The log log plot fixes the attention on the asymptotes, rather than the peak, and they are:

[TABLE]

Both of them are power laws. As long as the exponent of the high-frequency tail falls off slower than $\tau^{-2}$, i.e. for $\alpha<1$, the variance of the distribution will not exist. Such long-tailed distributions are scale-invariant and therefore may indicate fractal properties. This indicates that the distribution of elementary spring-damper models in the medium, as given in Figs. 1 and 3, may be a multi-fractal with two different fractal orders.

The references above didn’t derive the frequency-spectral function, but it can be found by substituting (18) and using the normalizing constant $a=J_{\tau}(\infty)=(1/E)(1-(\tau_{\sigma}/\tau_{\epsilon})^{\alpha})$ for the fractional Zener model. The creep frequency-spectral function of (19) is then:

[TABLE]

This result was first given in [22] in the form above. There it was derived with a starting point in the multiple relaxation theory of [23]. It should be noted that the frequency-spectral function is exactly the same as the time-spectral function except for a scaling factor.

The plot for $\alpha=0.8$ and $\alpha=0.01$ in Fig. 6 shares many properties with the plot of the time-spectral function. The asymptotes are also very similar:

[TABLE]

Interestingly, the relaxation spectrum approaches a single fractal for the limiting case $\alpha\rightarrow 0$, where both the low- and high-frequency parts will approach $\Omega^{-1}$. See the dash-dotted line for $\alpha=\text{\char 15\relax}=0.01$ in Fig. (6) with asymptotes $\Omega^{-1+\text{\char 15\relax}}$ for low frequencies and $\Omega^{-1-\text{\char 15\relax}}$ for high frequencies.

3.2 Relationship with hysteresis

The results found here have some resemblance with those derived from hysteresis. Hysteresis is a hypothetical loss element with a constant phase lag between stress and strain at all frequencies [24] resulting in attenuation that increases linearly with frequency [25]. The dynamic modulus for hysteresis is:

[TABLE]

where $K$ and $H$ are constants. The model leads to a noncausal model and therefore [24] is concerned with making a bandlimited approximation where the constants are allowed to vary with frequency, $K(\omega)$ and $H(\omega)$, in order to ensure causality.

By comparing the dynamic moduli, (28) and (23), it is evident that hysteresis can be viewed as the limiting case as $\alpha\rightarrow 0$ for the fractional Kelvin-Voigt model. Thus an alternative way to ensure causality is to employ a fractional Kelvin-Voigt model with parameter $\alpha=\text{\char 15\relax}=0^{+}$.

Hysteresis is particularly interesting in the context of fractional operators because the result for the limiting case $\alpha\rightarrow 0$ of the previous section is similar to that recently found for hysteresis in [26]. There it was shown that the hysteresis model is the same as a sum of relaxation processes weighted with a long-tailed power law $\Omega^{-1+\text{\char 15\relax}}$. This parallels our discussion of the asymptotic result of (27) for $\alpha=\text{\char 15\relax}$ and Parker’s result can be interpreted as a special case of (27).

4 Poroelastic model

The Biot poroelastic model deals with a saturated porous medium with a solid phase and a fluid phase. Wave propagation in such a medium is described by a set of coupled vector wave equations as given in Eq. (4.2) in [11]. The variables are the displacement vectors $\mathbf{u}$ and $\mathbf{U}$ for the displacement in the solid and the fluid respectively. That paper assumes that as the fluid moves in the pores, the flow is laminar and that losses are given by Darcy’s law and proportional to the relative velocity, $\partial(\mathbf{u-U})/\partial t$. The theory predicts three solutions, two compressional waves and one shear wave. In a biological porous medium like cancellous bone (bone with a low volume fraction of solid, less than 70%), all three waves have been detected [27].

This model is much more complex than the viscoelastic models of the previous sections. It also seems more appropriate for a complex medium like brain, liver, or bone. Here a surprising exact relationship between the poroelastic and viscoelastic models will be shown for the shear wave solution.

4.1 Biot’s original formulation

As elastography is only concerned with shear waves, we restrict the analysis here to that mode. Then the following dispersion relation was derived in [11], Eqs. (7.5)–(7.6):

[TABLE]

where

[TABLE]

This expression depends on three normalized densities:

[TABLE]

and the aggregate or composite density which is

[TABLE]

These formulas depend on the porosity, $\phi$, and the fluid and solid densities. The parameter $\rho_{12}$ represents a negative mass coupling density between fluid and solid, $\rho_{11}$ represents the total effective mass of the solid moving in the fluid, and $\rho_{22}$ is the total effective mass of the fluid. There is also a characteristic frequency

[TABLE]

where $\eta$ is the fluid viscosity, and $B$ is the permeability. The effective characteristic frequency of (30) is however a lower frequency:

[TABLE]

since $\rho_{22}=\alpha\phi\rho_{f}$ [28] where $\alpha$ is a structure constant related to tortuosity.

Comparison with (12) and (13) shows that (29) and (30) are exactly the same as those of the Zener model. This is a remarkable result which shows that for shear waves, the poroelastic model is also a linear viscoelastic model. What distinguishes it from the ordinary linear viscoelastic models is that it provides a sophisticated way of determining the parameters from physical considerations.

4.2 Biot-Stoll formulation

The parameters of the original Biot theory are often considered to be hard to find in practice and the theory has been rewritten in terms of the relative displacement between fluid and solid, $\mathbf{u-U}$, or the volume of fluid that has flowed in or out of an element, $\zeta=\phi\nabla(\mathbf{u-U})$, in combination with the solid displacement $\mathbf{u}$. The material parameters are then transformed to a new set of parameters.

In that case Eq. (16.77) in [29] gives the shear dispersion of the low frequency Biot model. It can also be found from Eq. (12) of [12]:

[TABLE]

where the sign of $\omega$ has been reversed compared to the original articles due to a different definition of the Fourier transform than here (see also Appendix A of [30]). The mass coupling density is

[TABLE]

where $\alpha$ is the tortuosity. Comparing (35) to (12) shows that the dynamic modulus is:

[TABLE]

As expected this is also equivalent to that of a Zener model as can be seen by comparison with (11). The constants when $\omega_{\epsilon}=1/\tau_{\epsilon}$ and $\omega_{\sigma}=1/\tau_{\sigma}$ are:

[TABLE]

Insertion of the expression for the mass coupling density, (36), in the equation for the effective characteristic frequency, (34), shows that this frequency also is the same as $\omega_{\epsilon}/(2\pi)$.

4.3 Redundant parameters in the Biot shear wave model

The low frequency Biot model is characterized by the ten parameters shown in Table 1. Seven of them affect shear wave propagation as indicated in the right-hand column, but several of them are connected as (38) indicates.

The parameter $\omega_{\epsilon}$ depends on the mass coupling density $\rho_{c}$. In addition it depends on the ratio of the viscosity, $\eta$, and the permeability, $B$. They do not influence any of the other parameters $\omega_{\sigma}$ and $c_{0}$ so it is clear that it is their ratio which matters. Likewise, the sound velocity, $c_{0}$, depends on the ratio of the shear modulus, $\mu_{r}$ and the aggregate density, $\rho$.

Furthermore for the third parameter, $\omega_{\sigma}$, the denominator in the expression is usually larger than $0.8$ so the ratio of $\omega_{\sigma}$ and $\omega_{\epsilon}$ is not very sensitive to changes in the densities, $\rho$, $\rho_{c}$, or $\rho_{f}$. That means that neither is it very sensitive to changes in tortuosity, and it will mainly be $\omega_{\epsilon}$ which depends on $\alpha$ in the ratio $\eta/(\alpha B)$.

In this way three of the seven parameters can be said to be redundant and the seven parameters from Table 1 ($\eta$, $B$, $\alpha$, $\rho_{f}$, $\rho_{s}$, $\phi$, and $\mu_{r}$) may be reduced to four by combining the three first ones into one, $\eta/(\alpha B)$. Alternatively, the four parameters may be stated as $\eta/\alpha$, $\rho_{c}$, $\rho$, and $\mu_{r}$. In case one wants to compute $\omega_{\sigma}$ exactly, a fifth parameter, $\rho_{f}$, also needs to be included.

In this way the number of parameters in the Biot shear model with some approximation is four as in the Zener medium model, and five in the exact case.

4.4 Extension to poroviscoelasticity

In the sediment acoustics field, the Biot model has been amended in order to extend the model from a rigid porous frame, like a porous rock or bone, to one where the grains are allowed to move. This is a model for a saturated sediment and it is not unlikely that it may be more appropriate for tissue than the rigid frame implied in the Biot model also.

In this case viscosity is introduced in the solid frame in addition to in the flowing liquid. The model is called the Biot squirt flow and viscous drag (BICSQS) model [31]. This viscosity is added by allowing a relaxation model for the shear modulus, thus allowing for a frequency dependent complex shear modulus or a dynamic shear modulus, described by a characteristic frequency $\omega_{\mu}$:

[TABLE]

When the relaxation is included in (37) and (11), the result is

[TABLE]

In the limit as frequency goes to zero the dynamic modulus approaches $E(0)=\mu$ and as the frequency approaches infinity it approaches $E(\infty)=\infty$. This model is called the non-standard Four-parameter model in [19] (Chap. 3.4.2). Its spring–damper realization is shown in Fig. 7 with the equivalent conjugate model to the right. Therefore, even in this case, an equivalent viscoelastic model may be found.

5 Discussion and Conclusion

The concept of the time- and frequency-spectral decompositions of viscoelastic systems has been developed and applied to the fractal Zener and Kelvin-Voigt models. This summation of multiple elementary models is an idea which is independent of the fractional models. In fact it is the idea behind the relaxation spectrum models of Kelvin and Wiechert dating from 1888 and 1893 respectively [15]. In addition to the coverage in [13], the book [19] devotes several sections to it (Chap. 3.5-3.6). Similar ideas of model fitting have been used in biomechanics and elastography also. In the white matter model of [32] the shear relaxation modulus was for instance modeled with three very slow relaxation terms (0.01 Hz and slower) and a similar model with two terms was used in [33].

The idea behind the particular weightings implied in the time-spectral and frequency-spectral functions found here is however to make the sum approximate the behavior of the fractional models, i.e. in terms of power laws in the frequency domain. Therefore it parallels the modeling of arbitrary power law attenuation in medical ultrasound over a limited frequency range with a few terms in [34].

The surprising result is that the weighting has a long-tailed distribution. This is a property which is associated with a fractal geometry. Here it characterizes the distribution of individual relaxation processes both in the time and the frequency domains. An interesting special case is the hysteresis model recently proposed by [26] and this could possibly shed some light on the interpretation of the distribution. The fractal distribution is an explanation that seemingly is based on a completely different physical mechanism than the link shown between the fractional damper and a linearly time-varying viscosity in [35]. But perhaps the link is even deeper and suggests that such a time-varying viscosity in some way is associated with a medium with a fractal distribution of relaxation processes?

The connection between the fractal distribution and the fractional constitutive equation of an absorbing medium is also a result that complements recent results for the alternative attenuation mechanism due to scattering rather than absorption. An example is [36, 37] where it is shown that a fractal distribution of random scatterers leads to attenuation that varies with frequency according to a power law. The fractal – fractional derivative connection can also be compared to the geometrical and physical interpretation of fractional integration and differentiation in [38], but the result here is more specific to one particular application, namely viscoelasticity.

The poroelastic model depends on ten independent parameters, of which seven influence the shear wave solution, and seems to build on an entirely different foundation than the viscoelastic models. Despite that it has been shown to be equivalent to a Zener model. This was done by comparing dispersion relations, and not unsurprisingly the result is the same whether the dispersion relations from the original Biot theory or those from the Biot-Stoll theory are analyzed. This result is an extension of [39] where it was shown that in the low loss/low frequency approximation, an equivalent can be found between the poroelastic model and a Kelvin-Voigt model. That was done by comparing approximate expressions for the velocity and attenuation. As the Zener model in the low loss/low frequency case is equivalent to the Kelvin-Voigt model [16], the result found here also agrees with that result.

It is also shown that the seven parameters of the poroelastic model can be reduced to five, and even four as in the Zener model if a small approximation is allowed. When the Biot model is extended to include viscosity in the frame as in [31], an extra damper has to be added to the Zener model making it into what is called the non-standard four-parameter model in [19].