Fractional Viscoelastic Wave Propagation

• Fractional viscoelastic wave propagation models use fractional derivatives and distributed-order kernels to represent material memory and frequency-dependent attenuation.
• The framework couples linearized momentum balance with fractional constitutive laws, linking microstructural physics to macroscopic dispersion and causality.
• Numerical techniques like diffusive representations and memory compression enable efficient simulation of broadband attenuation in fields such as geophysics and medical imaging.

A fractional viscoelastic wave propagation model is a wave model in which constitutive memory is represented by fractional operators, distributed-order kernels, or related nonlocal operators, so that attenuation, dispersion, and relaxation are governed by power-law or Mittag–Leffler memory rather than by a single exponential timescale. In its canonical form, the framework couples linearized momentum balance and strain–displacement relations to a fractional constitutive law, yielding causal wave equations with frequency-dependent attenuation and dispersion. The same mathematical structure appears in several distinct settings: fractional Zener and Kelvin–Voigt rheologies, distributed-order and Burgers media, poroviscoelastic and constant-$Q$ formulations, and grain-shearing models for marine sediments in which the fractional order is tied to identifiable microphysics rather than introduced ad hoc (Nasholm et al., 2012, Pandey et al., 2015).

1. Constitutive foundation and memory operators

The constitutive core of the subject is the replacement of Hookean or standard linear-solid stress–strain laws by hereditary relations containing fractional derivatives. A widely used form is the fractional Zener law

$$\sigma(t)+\tau_\epsilon^{\beta}\frac{\partial^\beta \sigma(t)}{\partial t^\beta} = E_0\left[\epsilon(t)+\tau_\sigma^{\alpha}\frac{\partial^\alpha \epsilon(t)}{\partial t^\alpha}\right],$$

with elastic modulus $E_0$, characteristic times $\tau_\epsilon,\tau_\sigma$, and fractional orders $\alpha,\beta$. In the physically preferred Cole–Cole-type case one imposes $\alpha=\beta$, which the literature associates with thermodynamic admissibility and monotone relaxation. When this law is combined with $\epsilon=\nabla u$ and $\nabla \sigma=\rho_0 u_{tt}$, the material memory passes directly into the wave operator through the generalized compressibility $\kappa(\omega)=\epsilon(\omega)/\sigma(\omega)$ and the dispersion relation $k^2(\omega)=\omega^2\rho_0\kappa(\omega)$ (Nasholm et al., 2012).

The memory mechanism is explicit at the operator level. In Caputo form, the fractional derivative is a convolution with a weakly singular kernel,

$$\sigma(t)+\tau_\epsilon^{\beta}\frac{\partial^\beta \sigma(t)}{\partial t^\beta} = E_0\left[\epsilon(t)+\tau_\sigma^{\alpha}\frac{\partial^\alpha \epsilon(t)}{\partial t^\alpha}\right],$$0

so the constitutive law is non-Markovian even when the balance laws remain local. This operator viewpoint is central in formulations that reformulate the stress law into an integro-differential equation with a Mittag–Leffler kernel. For example, one can invert a Caputo constitutive equation by Laplace transform and obtain a stress history integral with kernel $$\sigma(t)+\tau_\epsilon^{\beta}\frac{\partial^\beta \sigma(t)}{\partial t^\beta} = E_0\left[\epsilon(t)+\tau_\sigma^{\alpha}\frac{\partial^\alpha \epsilon(t)}{\partial t^\alpha}\right],$$1, which is then approximated numerically by a sum of exponentials without changing the underlying hereditary interpretation (Yuan et al., 2024).

Fractional constitutive laws also admit an equivalent relaxation-spectrum interpretation. The fractional Zener model can be written as a continuum Maxwell–Wiechert or Nachman–Smith–Waag representation, in which a distributed relaxation spectrum reproduces the same complex compressibility as the fractional law. This equivalence is conceptually important because it places fractional rheology alongside conventional springs-and-dashpots models rather than outside them. A plausible implication is that the adjective “fractional” describes the closure used for a broadband relaxation spectrum, not a departure from continuum viscoelasticity itself (Nasholm et al., 2012).

2. Canonical fractional wave equations

The most common time-fractional elastic wave equation derived from the fractional Zener constitutive law is

$$\sigma(t)+\tau_\epsilon^{\beta}\frac{\partial^\beta \sigma(t)}{\partial t^\beta} = E_0\left[\epsilon(t)+\tau_\sigma^{\alpha}\frac{\partial^\alpha \epsilon(t)}{\partial t^\alpha}\right],$$2

This four-parameter model is causal, implies frequency power-law attenuation, and reduces to the fractional Kelvin–Voigt model when $$\sigma(t)+\tau_\epsilon^{\beta}\frac{\partial^\beta \sigma(t)}{\partial t^\beta} = E_0\left[\epsilon(t)+\tau_\sigma^{\alpha}\frac{\partial^\alpha \epsilon(t)}{\partial t^\alpha}\right],$$3. In one-dimensional form it also appears as a convolutional wave equation

$$\sigma(t)+\tau_\epsilon^{\beta}\frac{\partial^\beta \sigma(t)}{\partial t^\beta} = E_0\left[\epsilon(t)+\tau_\sigma^{\alpha}\frac{\partial^\alpha \epsilon(t)}{\partial t^\alpha}\right],$$4

which is the formulation used in Cauchy-problem analyses of fractional Zener rods (Nasholm et al., 2012, Konjik et al., 2011).

A second canonical equation is the time-fractional diffusion-wave equation,

$$\sigma(t)+\tau_\epsilon^{\beta}\frac{\partial^\beta \sigma(t)}{\partial t^\beta} = E_0\left[\epsilon(t)+\tau_\sigma^{\alpha}\frac{\partial^\alpha \epsilon(t)}{\partial t^\alpha}\right],$$5

which interpolates continuously between diffusion $$\sigma(t)+\tau_\epsilon^{\beta}\frac{\partial^\beta \sigma(t)}{\partial t^\beta} = E_0\left[\epsilon(t)+\tau_\sigma^{\alpha}\frac{\partial^\alpha \epsilon(t)}{\partial t^\alpha}\right],$$6 and the wave equation $$\sigma(t)+\tau_\epsilon^{\beta}\frac{\partial^\beta \sigma(t)}{\partial t^\beta} = E_0\left[\epsilon(t)+\tau_\sigma^{\alpha}\frac{\partial^\alpha \epsilon(t)}{\partial t^\alpha}\right],$$7. In viscoelastic wave propagation this equation arises when the creep law is a power law and the constitutive relation becomes

$$\sigma(t)+\tau_\epsilon^{\beta}\frac{\partial^\beta \sigma(t)}{\partial t^\beta} = E_0\left[\epsilon(t)+\tau_\sigma^{\alpha}\frac{\partial^\alpha \epsilon(t)}{\partial t^\alpha}\right],$$8

For $$\sigma(t)+\tau_\epsilon^{\beta}\frac{\partial^\beta \sigma(t)}{\partial t^\beta} = E_0\left[\epsilon(t)+\tau_\sigma^{\alpha}\frac{\partial^\alpha \epsilon(t)}{\partial t^\alpha}\right],$$9, the model describes “fractional diffusive waves,” an intermediate regime that is neither purely parabolic nor purely hyperbolic (Mainardi, 2012).

Several important generalizations enlarge this canonical set. The space-time fractional Zener equation replaces the local strain $E_0$0 by a symmetrized fractional nonlocal strain $E_0$1 and yields

$E_0$2

where $E_0$3 is a convolution operator with Laplace symbol $E_0$4. The model includes finite-speed waves for local or weakly nonlocal regimes and a non-propagating disturbance when $E_0$5, since then $E_0$6 and the solution degenerates to $E_0$7 (Atanackovic et al., 2014).

Other constitutive classes produce different wave equations and support properties. Thermodynamically consistent fractional Burgers models divide into two classes. One class produces kernels with full spatial support and therefore infinite propagation speed; the other yields compact support $E_0$8 with

$E_0$9

Distributed-order constitutive laws likewise lead to a fractional wave equation of the form

$\tau_\epsilon,\tau_\sigma$0

with the wavefront speed tied to the glass modulus by $\tau_\epsilon,\tau_\sigma$1 when $\tau_\epsilon,\tau_\sigma$2 is finite (Oparnica et al., 2019, Konjik et al., 2017).

3. Physical derivations and microstructural interpretations

A major theme in the literature is that fractional wave equations need not be interpreted as phenomenological curve-fits only. In grain-shearing models for saturated, unconsolidated marine sediments, the medium consists of grains that develop micro-asperities under overburden pressure and then undergo stick-slip sliding through the pore fluid. The viscous resistance increases with shearing time, a behavior identified as rheopecty, and the time-dependent dashpot viscosity is modeled as

$\tau_\epsilon,\tau_\sigma$3

The resulting material impulse response behaves asymptotically as

$\tau_\epsilon,\tau_\sigma$4

which is precisely the power-law kernel required by fractional calculus. Under the GS simplification $\tau_\epsilon,\tau_\sigma$5, the compressional equation becomes a Kelvin–Voigt fractional-derivative wave equation and the shear equation becomes a time-fractional diffusion-wave equation. In this interpretation the order

$\tau_\epsilon,\tau_\sigma$6

measures the competition between elastic stiffness and strain-hardening rate, so the fractional order is a microstructural descriptor of grain-shearing rather than an arbitrary fitting parameter (Pandey et al., 2015, Pandey et al., 2016).

A related physical derivation appears in high-frequency poroelasticity. In the Biot–JKD model, the dynamic permeability correction has square-root frequency dependence, and in the time domain the viscous drag becomes a shifted half-order derivative,

$\tau_\epsilon,\tau_\sigma$7

The operator has a convolution kernel proportional to $\tau_\epsilon,\tau_\sigma$8 and admits a diffusive representation as a continuum of exponentially decaying memory variables satisfying local ODEs. This places fractional wave propagation inside the established Biot framework for fluid-saturated porous media rather than outside it (Blanc et al., 2012).

Constant-$\tau_\epsilon,\tau_\sigma$9 seismology provides a third interpretation. In Kjartansson-type fractional models the constitutive law is written directly in terms of Caputo derivatives of order $\alpha,\beta$0, with

$\alpha,\beta$1

Because realistic earth materials may have very small $\alpha,\beta$2, the corresponding Caputo kernels are very long-memory. The same papers show that a sum-of-exponentials approximation of the fractional kernel is equivalent to a generalized Maxwell body with finitely many relaxation times. This suggests that constant-$\alpha,\beta$3 fractional rheology and generalized Maxwell representations are not competing descriptions but two levels of resolution of the same hereditary mechanism (Guo et al., 2023).

4. Attenuation, dispersion, causality, and propagation speed

Fractional viscoelastic wave models are primarily used because they reproduce coupled attenuation and dispersion over broad bands. In the fractional Zener case with $\alpha,\beta$4, the attenuation coefficient derived from $\alpha,\beta$5 exhibits three asymptotic power-law regimes,

$\alpha,\beta$6

This multiregime structure is one reason the model is used in acoustics, ultrasound, photoacoustics, and elastography, where a single exponent rarely describes the entire band. The same formulation is causal, and its frequency-domain structure is compatible with Kramers–Kronig relations (Nasholm et al., 2012).

Propagation speed is more model-dependent than the phrase “fractional wave equation” might suggest. Some fractional constitutive classes retain finite wavefronts, while others do not. In distributed-order media and in certain fractional Zener-type formulations, the support of the fundamental solution is conic and the wave speed is controlled by the instantaneous material response, specifically the glass modulus $\alpha,\beta$7, with $\alpha,\beta$8. In the second class of fractional Burgers models the same principle appears through the finite glass modulus $\alpha,\beta$9, giving $\alpha=\beta$0. By contrast, the first Burgers class has zero glass compliance and a kernel with full support on $\alpha=\beta$1, implying infinite propagation speed (Konjik et al., 2017, Oparnica et al., 2019).

These distinctions are mirrored in broader comparative studies of Zener, Maxwell, and Voigt models, both classical and fractional. Newton’s model becomes parabolic after substitution and does not produce waves, while Zener, Maxwell, and Voigt formulations are dissipative wave models; their fractional counterparts preserve this qualitative distinction but alter dissipation rate and dispersion. Numerical comparisons report that decreasing the fractional order slows dissipation, fractional Maxwell shows slower creep than classical Maxwell, and Voigt models show strong dispersion (Brown et al., 2018).

Admissibility is likewise nontrivial. For the fractional Zener law, the $\alpha=\beta$2 case is singled out as the thermodynamically admissible one because the corresponding relaxation spectrum can be kept nonnegative. In finite-deformation incompressible settings, several fractional viscous stresses that appear natural are nevertheless ruled out by failure of frame-indifference or thermodynamic consistency; only certain Kelvin–Voigt-type constructions survive both tests. A plausible implication is that model selection cannot be based on attenuation fits alone, especially when one departs from small-strain linear theory (Nasholm et al., 2012, Berjamin et al., 2023).

5. Fundamental solutions, Cauchy problems, and analytical structure

A substantial part of the theory concerns the generalized Cauchy problem for nonlocal wave operators. In the one-dimensional fractional Zener rod, the problem is formulated distributionally as

$\alpha=\beta$3

with

$\alpha=\beta$4

For $\alpha=\beta$5, the model admits a unique solution represented as convolution with a fundamental solution $\alpha=\beta$6. The support of $\alpha=\beta$7 lies in the cone $\alpha=\beta$8, which provides a finite-speed interpretation in scaled variables (Konjik et al., 2011).

Existence and uniqueness have also been established for more general nonlocal operators. For the space-time fractional Zener equation, generalized solutions are constructed in $\alpha=\beta$9 with support in $\epsilon=\nabla u$0, and the solution is written as

$\epsilon=\nabla u$1

where $\epsilon=\nabla u$2 is a transform-defined kernel obtained as a distributional limit of regularized kernels. The analysis identifies the symbol

$\epsilon=\nabla u$3

whose zeros govern inversion and residue calculations (Atanackovic et al., 2014).

Distributed-order and Burgers-type models yield analogous transform formulas but different support properties. In distributed-order viscoelasticity the Laplace-domain kernel has the form

$\epsilon=\nabla u$4

and explicit inversion shows finite wavefronts whenever the glass modulus is finite. In fractional Burgers media, the analytical form of the kernel depends on the branching structure of the constitutive symbol $\epsilon=\nabla u$5; the support can be either full or conic, and Dirac-initial-data profiles may exhibit secondary peaks and endpoint behavior that the authors describe as not usually expected in standard wave propagation (Konjik et al., 2017, Oparnica et al., 2019).

Recent variable-order models alter the analysis by making the fractional order time-dependent. One formulation studies

$\epsilon=\nabla u$6

and proves well-posedness by rewriting it as a local-in-time wave equation with a convolution forcing term,

$\epsilon=\nabla u$7

This local reformulation removes a nonphysical initial singularity identified in earlier variable-order work and supports high-order regularity estimates needed for numerical analysis (Jia et al., 8 Nov 2025).

6. Numerical realization and memory compression

The major numerical obstacle is the history term: direct discretization of a Caputo or Riemann–Liouville derivative requires storage and summation over the full past. Two broad strategies dominate the literature. The first is diffusive representation, in which a power-law kernel is expressed as a superposition of decaying exponentials and replaced by memory variables satisfying local ODEs. In the Biot–JKD model the shifted half-order derivative is written as

$\epsilon=\nabla u$8

with

$\epsilon=\nabla u$9

Quadrature then yields a finite-memory Biot-DA model solved by Strang splitting, a fourth-order ADER scheme for the propagative part, and exact integration of the relaxation step. In reported one-dimensional tests, $\nabla \sigma=\rho_0 u_{tt}$0 diffusive variables produced a modeling error of about $\nabla \sigma=\rho_0 u_{tt}$1, a numerical error around $\nabla \sigma=\rho_0 u_{tt}$2, and a global convergence rate about $\nabla \sigma=\rho_0 u_{tt}$3 (Blanc et al., 2012).

The Andrade model uses the same philosophy. Its fractional derivative of order $\nabla \sigma=\rho_0 u_{tt}$4 is represented by diffusive variables $\nabla \sigma=\rho_0 u_{tt}$5 satisfying

$\nabla \sigma=\rho_0 u_{tt}$6

and the wave system is advanced with Strang splitting, exact integration of the diffusive subsystem, and an ADER-4 propagative step. The local PDE system is well posed when quadrature weights are positive, and constrained optimization of the quadrature nodes and weights gives much better agreement with the target Andrade compliance than classical Gauss–Jacobi constructions (Jazia et al., 2013).

The second dominant strategy is sum-of-exponentials compression of Mittag–Leffler or power-law kernels. For constant-$\nabla \sigma=\rho_0 u_{tt}$7 viscoelasticity, a nearly optimal SOE approximation

$\nabla \sigma=\rho_0 u_{tt}$8

is constructed using generalized Gaussian quadrature. The resulting memory variables satisfy local ODEs, and in three-dimensional homogeneous viscoacoustic simulations the method achieved relative error below $\nabla \sigma=\rho_0 u_{tt}$9 with $\kappa(\omega)=\epsilon(\omega)/\sigma(\omega)$0 for $\kappa(\omega)=\epsilon(\omega)/\sigma(\omega)$1; in the full viscoelastic case, values such as $\kappa(\omega)=\epsilon(\omega)/\sigma(\omega)$2 and $\kappa(\omega)=\epsilon(\omega)/\sigma(\omega)$3 yielded about $\kappa(\omega)=\epsilon(\omega)/\sigma(\omega)$4 relative error while storing 62 extra memory wavefields (Guo et al., 2023).

Finite element realizations of SOE acceleration follow the same principle. One mixed finite element/Newmark scheme for a fractional viscoelastic stress law rewrites the constitutive equation into a Mittag–Leffler convolution and then approximates

$\kappa(\omega)=\epsilon(\omega)/\sigma(\omega)$5

reducing memory complexity from $\kappa(\omega)=\epsilon(\omega)/\sigma(\omega)$6 to $\kappa(\omega)=\epsilon(\omega)/\sigma(\omega)$7 and computation complexity from $\kappa(\omega)=\epsilon(\omega)/\sigma(\omega)$8 to $\kappa(\omega)=\epsilon(\omega)/\sigma(\omega)$9, with $k^2(\omega)=\omega^2\rho_0\kappa(\omega)$0. A related fast backward Euler finite element method for a generalized fractional Maxwell–Zener model derives the fully discrete error bound $k^2(\omega)=\omega^2\rho_0\kappa(\omega)$1 and reports the same asymptotic cost reduction. Variable-order formulations exploit Toeplitz translational invariance instead, using a divide-and-conquer algorithm to reduce computational complexity from $k^2(\omega)=\omega^2\rho_0\kappa(\omega)$2 to $k^2(\omega)=\omega^2\rho_0\kappa(\omega)$3 (Yuan et al., 2024, Yuan et al., 16 Jul 2025, Jia et al., 8 Nov 2025).

7. Applications, extensions, and disputed points

Fractional viscoelastic wave propagation is used across acoustics, geophysics, and imaging because it combines broadband attenuation with compact constitutive parameterizations. The fractional Zener formulation has been emphasized for medical ultrasound, photoacoustic tomography, MR elastography, and related imaging modalities in which soft tissue exhibits power-law attenuation and dispersion. In geophysics, constant-$k^2(\omega)=\omega^2\rho_0\kappa(\omega)$4 time-fractional models are used for 3D seismic simulation and inversion because they capture both amplitude loss and phase delay while avoiding large parameter sets typical of discrete multi-relaxation models (Nasholm et al., 2012, Guo et al., 2023).

Marine sediments provide a notable physically grounded application. Grain-shearing models for saturated, unconsolidated sediments show that compressional waves obey a Kelvin–Voigt fractional-derivative equation and shear waves obey a fractional diffusion-wave equation. The latter predicts coexistence of diffusion and wave propagation in the shear channel, a result that the fractional formulation captures exactly. This is one of the clearest cases in which a fractional constitutive description is explicitly derived from a grain-scale mechanism (Pandey et al., 2015).

Extensions continue in several directions. Space-fractional strain measures introduce long-range spatial interactions and can drive the model to a non-propagating disturbance as $k^2(\omega)=\omega^2\rho_0\kappa(\omega)$5 (Atanackovic et al., 2014). Variable-order formulations let the memory strength evolve with material state or time, producing what the literature describes as diffusive wave propagation in viscoelastic media with evolving physical property (Jia et al., 8 Nov 2025). Another line replaces the singular power-law kernel by the Atangana–Baleanu–Caputo derivative with a non-singular Mittag–Leffler kernel,

$k^2(\omega)=\omega^2\rho_0\kappa(\omega)$6

and reports non-exponential, Mittag–Leffler-like energy decay in one-dimensional seismic-wave simulations (Demir et al., 15 Dec 2025).

Several issues remain active points of differentiation rather than settled doctrine. One is whether a particular fractional law preserves finite propagation speed; the literature gives affirmative and negative examples depending on constitutive class, glass modulus, and kernel structure. Another is whether a proposed fractional stress definition is physically admissible: some finite-strain formulations are rejected on frame-indifference or thermodynamic grounds even though they reduce to plausible small-strain forms. A third is interpretive: some models derive the fractional order from microphysics, relaxation spectra, or $k^2(\omega)=\omega^2\rho_0\kappa(\omega)$7-theory, whereas others leave it as an effective parameter. Taken together, these results indicate that “fractional viscoelastic wave propagation model” is not a single equation but a class of hereditary wave theories unified by nonlocal constitutive memory and differentiated by how that memory is justified, constrained, analyzed, and discretized (Oparnica et al., 2019, Berjamin et al., 2023)