Westervelt Equation Overview

• Westervelt equation is a quasilinear PDE modeling finite-amplitude sound propagation in viscous media with an amplitude-dependent inertia coefficient.
• It features strong damping and versatile formulations—including pressure, velocity potential, and periodic models—enabling robust energy analysis and simulation.
• Its study underpins advanced numerical methods, inverse problem strategies, and ultrasound imaging techniques in nonlinear acoustics.

The Westervelt equation is a quasilinear model of finite-amplitude sound propagation in viscous or thermoviscous media. In a standard pressure formulation with acoustic pressure fluctuation $u(t,x)$, sound speed $c>0$, diffusivity $b>0$, and nonlinearity parameter $k$, it is written as

$u_{tt}-c^2\Delta u-b\,\Delta u_t=k\,(u^2)_{tt},$

or equivalently,

$$(1-2ku)\,u_{tt}-c^2\Delta u-b\,\Delta u_t=2k\,(u_t)^2.$$

This equivalent form exhibits the quasilinear dependence of the principal time-derivative coefficient on the acoustic amplitude, a feature that governs both analysis and computation. Across recent literature, the equation appears in pressure, velocity-potential, periodic, absorbing-boundary, thermoacoustic, and time-fractional variants, and it serves as a central model for nonlinear acoustics, ultrasound simulation, and parameter identification (Meyer et al., 2015).

1. Classical formulations and model variants

In nonlinear acoustics, the pressure formulation is often expressed as

$$p_{tt}-c^2\Delta p-b\,\Delta p_t=\partial_t^2\!\left(\frac{\beta}{\rho c^2}p^2\right),$$

or, after setting $k=\beta/(\rho c^2)$,

$u_{tt}-c^2\Delta u-b\,\Delta u_t=k\,(u^2)_{tt}.$

Using $(u^2)_{tt}=2u\,u_{tt}+2(u_t)^2$, one obtains the quasilinear form $c>0$0, which makes the effective inertia explicit (Wilke, 2022).

A complementary formulation uses the acoustic velocity potential $c>0$1. In dimensionless variables, one form studied in the literature is

$c>0$2

with acoustic velocity $c>0$3 and pressure variation $c>0$4. This formulation is particularly convenient for energy-preserving discretization, because it admits a canonical gradient-system structure and an associated continuous energy identity (Egger et al., 2019).

The equation is closely related to the Kuznetsov equation but is not identical to it. The Westervelt model arises as a reduced nonlinear acoustics equation under appropriate simplifications, including the neglect of convective terms and retention of strong damping through a $c>0$5 contribution (Wilke, 2022). In periodic settings with Robin boundary conditions, the equation is also written as

$c>0$6

with time-periodic conditions $c>0$7, $c>0$8, and impedance-type boundary terms $c>0$9, thereby linking nonlinear propagation to multiharmonic analysis (Rainer et al., 2024).

Although the classical Westervelt equation is rooted in fluid acoustics, a recent derivation from deformation theory extends the model to nonlinear lossy solids and tissues. In that framework, the pressure equation takes the generalized form

$b>0$0

with coefficients determined by elastic moduli, geometric nonlinearity, and dissipative parameters. This places the classical Westervelt equation as a leading-order member of a broader hierarchy of nonlinear pressure-wave models (Caruso et al., 17 Mar 2025).

2. Quasilinearity, damping, and energy structure

The defining structural feature of the Westervelt equation is the amplitude-dependent coefficient in front of $b>0$1. In the pressure form

$b>0$2

or, in a normalized variant,

$b>0$3

the equation loses parabolic-type coercivity when $b>0$4 or $b>0$5 approaches zero. For this reason, smallness conditions that keep the coefficient strictly positive are not merely technical conveniences; they are structural conditions preventing degeneracy (Simonett et al., 2016).

The strong damping term $b>0$6 or $b>0$7 gives the equation a mixed hyperbolic-parabolic character. This is the mechanism behind maximal-regularity results, temporal smoothing, and exponential decay in several analytical works. It also explains why the Westervelt equation behaves differently from nondiffusive weakly nonlinear wave models (Meyer et al., 2015).

Energy identities make this structure transparent. In the velocity-potential formulation, one energy is

$b>0$8

and under homogeneous Neumann boundary conditions it satisfies

$b>0$9

Thus, exact conservation occurs when $k$0, while $k$1 yields viscous dissipation (Egger et al., 2019).

For absorbing-boundary pressure formulations, a natural energy is

$k$2

and the boundary law

$k$3

produces boundary dissipation in addition to the interior term $k$4 (Simonett et al., 2016).

In the nondiffusive regime, the weakly nonlinear geometric-optics limit is different. For

$k$5

the leading profile equation is of Burgers type,

$k$6

so wave steepening, rather than parabolic smoothing, governs the principal asymptotics (Eptaminitakis et al., 2022).

3. Initial-boundary value problems and long-time dynamics

For bounded $k$7-domains with homogeneous Dirichlet boundary conditions, the Westervelt equation admits a global small-data theory in an optimal maximal-$k$8-regularity framework. If the initial data are sufficiently small and satisfy the appropriate compatibility conditions, there exists a unique global solution

$k$9

and the solution converges to zero exponentially fast as $u_{tt}-c^2\Delta u-b\,\Delta u_t=k\,(u^2)_{tt},$0 (Meyer et al., 2015).

A more refined analysis with nonlinear absorbing boundary conditions of order zero establishes global well-posedness for small initial data, instantaneous temporal regularization, and convergence to equilibria. In that setting the full problem is

$u_{tt}-c^2\Delta u-b\,\Delta u_t=k\,(u^2)_{tt},$1

and each equilibrium belongs to

$u_{tt}-c^2\Delta u-b\,\Delta u_t=k\,(u^2)_{tt},$2

Solutions starting sufficiently close to an equilibrium converge exponentially to a possibly different equilibrium, reflecting the one-dimensional center manifold generated by constants (Simonett et al., 2016).

Instantaneous smoothing is a recurring feature. In the absorbing-boundary theory, for every $u_{tt}-c^2\Delta u-b\,\Delta u_t=k\,(u^2)_{tt},$3 and $u_{tt}-c^2\Delta u-b\,\Delta u_t=k\,(u^2)_{tt},$4,

$u_{tt}-c^2\Delta u-b\,\Delta u_t=k\,(u^2)_{tt},$5

hence $u_{tt}-c^2\Delta u-b\,\Delta u_t=k\,(u^2)_{tt},$6 (Simonett et al., 2016). A parallel parameter-trick argument in the non-isothermal Westervelt–Pennes system yields analogous instantaneous temporal regularization for both acoustic and thermal variables (Wilke, 2022).

These analytical results also clarify a recurrent misconception. The equation is sometimes described as a weakly nonlinear wave equation with lower-order corrections, but the main highest-order coefficient depends on the solution itself, and the damping term is analytically decisive. The small-data restriction therefore preserves the equation’s admissible regime rather than merely simplifying the proofs (Meyer et al., 2015).

4. Boundary conditions, periodic forcing, and coupled thermoacoustic models

Artificial boundary modeling is a major part of Westervelt-equation analysis. A family of nonlinear absorbing boundary conditions was constructed by pseudo-differential calculus in one and two space dimensions. The zero-order condition

$u_{tt}-c^2\Delta u-b\,\Delta u_t=k\,(u^2)_{tt},$7

acts as a nonlinear amplitude-dependent impedance. First-order conditions introduce tangential derivatives and memory-type terms, and numerical experiments showed that the pseudo-differential conditions were more accurate than classical Engquist–Majda absorbing boundary conditions in the tested regimes (Kaltenbacher et al., 2014).

Periodic excitation leads naturally to multiharmonic formulations. For a bounded domain with Robin boundary conditions,

$u_{tt}-c^2\Delta u-b\,\Delta u_t=k\,(u^2)_{tt},$8

and time-periodic forcing, the nonlinear periodic Westervelt equation generates higher harmonics through quadratic frequency autoconvolution. The corresponding truncated multiharmonic system consists of coupled Helmholtz-type equations for harmonic amplitudes $u_{tt}-c^2\Delta u-b\,\Delta u_t=k\,(u^2)_{tt},$9, and an iterative harmonic-by-harmonic solution scheme converges for sufficiently small excitation (Rainer et al., 2024).

The equation also appears in thermoacoustic couplings. In the non-isothermal Westervelt–Pennes model,

$$(1-2ku)\,u_{tt}-c^2\Delta u-b\,\Delta u_t=2k\,(u_t)^2.$$0

the acoustic coefficients depend on temperature, and the pressure couples to a bioheat equation through a source term $$(1-2ku)\,u_{tt}-c^2\Delta u-b\,\Delta u_t=2k\,(u_t)^2.$$1. Local and global well-posedness, instantaneous regularization, and exponential decay near equilibrium were established in an $$(1-2ku)\,u_{tt}-c^2\Delta u-b\,\Delta u_t=2k\,(u_t)^2.$$2-$$(1-2ku)\,u_{tt}-c^2\Delta u-b\,\Delta u_t=2k\,(u_t)^2.$$3 maximal-regularity framework (Wilke, 2022).

Replacing Fourier heat conduction with the hyperbolic Cattaneo law yields the Westervelt–Pennes–Cattaneo model. Local well-posedness and nondegeneracy under small pressure data were proved together with a singular limit $$(1-2ku)\,u_{tt}-c^2\Delta u-b\,\Delta u_t=2k\,(u_t)^2.$$4, in which the Cattaneo system converges to the Fourier-coupled Westervelt–Pennes model (Benabbas et al., 2023). A later global theory showed existence and exponential decay under small lower-order Sobolev norms even when higher-order norms are arbitrarily large, provided the acoustic diffusivity $$(1-2ku)\,u_{tt}-c^2\Delta u-b\,\Delta u_t=2k\,(u_t)^2.$$5 is present (Benabbas et al., 2024).

Time-fractional attenuation provides another extension. In the fractionally damped model

$$(1-2ku)\,u_{tt}-c^2\Delta u-b\,\Delta u_t=2k\,(u_t)^2.$$6

inhomogeneous Neumann boundary control and distributed control were analyzed, and existence of globally optimal controls together with adjoint-based first-order optimality conditions were derived (Nikolić et al., 19 Nov 2025).

5. Numerical discretization and computational methods

Several numerical frameworks are now established for the Westervelt equation. For homogeneous Dirichlet problems, a semi-discrete $$(1-2ku)\,u_{tt}-c^2\Delta u-b\,\Delta u_t=2k\,(u_t)^2.$$7 finite element approximation yields optimal a priori spatial error bounds in $$(1-2ku)\,u_{tt}-c^2\Delta u-b\,\Delta u_t=2k\,(u_t)^2.$$8-based norms under sufficiently small data and mesh size. The analysis combines a linear variable-coefficient wave estimate with a Banach fixed-point argument and uses inverse estimates together with Scott–Zhang interpolation to prevent degeneration of the discrete coefficient $$(1-2ku)\,u_{tt}-c^2\Delta u-b\,\Delta u_t=2k\,(u_t)^2.$$9 (Nikolić et al., 2019).

Structure-preserving discretization has also been developed in the velocity-potential formulation. A canonical gradient-system reformulation supports space-time Galerkin–Petrov approximations of formally arbitrary order, and with a consistent quadrature-based inner product the scheme preserves a discrete energy identity exactly. In the linear case $$p_{tt}-c^2\Delta p-b\,\Delta p_t=\partial_t^2\!\left(\frac{\beta}{\rho c^2}p^2\right),$$0, the method coincides with Lobatto–IIIA; for $$p_{tt}-c^2\Delta p-b\,\Delta p_t=\partial_t^2\!\left(\frac{\beta}{\rho c^2}p^2\right),$$1, it preserves the nonlinear energy identity whereas Lobatto–IIIA, Gauß–Runge–Kutta, and Newmark do not (Egger et al., 2019).

Time integration by operator splitting exploits the hyperbolic-parabolic decomposition induced by the Westervelt nonlinearity and damping. For the first-order Lie–Trotter splitting applied to the velocity-potential form, a global error estimate confirms stability and first-order convergence under sufficient regularity, and numerical experiments show that the favored decomposition reduces computational cost by separating a nonlinear diffusion subproblem from a simpler wave-like subproblem (Kaltenbacher et al., 2013).

For curved geometries and manifolds, discrete exterior calculus provides a coordinate-free spatial discretization. The DEC formulation approximates the Laplace–Beltrami operator on simplicial complexes and leads to explicit, implicit, and semi-implicit schemes for

$$p_{tt}-c^2\Delta p-b\,\Delta p_t=\partial_t^2\!\left(\frac{\beta}{\rho c^2}p^2\right),$$2

This extends finite-difference-type ultrasound simulation to triangulated and tetrahedralized manifolds (Xie et al., 2010).

For nonlocal attenuation, a semi-discretization in time combines the trapezoidal rule with A-stable convolution quadrature. Because the fractional kernel produces a weak singularity at $$p_{tt}-c^2\Delta p-b\,\Delta p_t=\partial_t^2\!\left(\frac{\beta}{\rho c^2}p^2\right),$$3, an initial correction is added to restore the optimal rate compatible with the weak-start behavior. Without correction the energy-norm error is $$p_{tt}-c^2\Delta p-b\,\Delta p_t=\partial_t^2\!\left(\frac{\beta}{\rho c^2}p^2\right),$$4, whereas the corrected scheme achieves $$p_{tt}-c^2\Delta p-b\,\Delta p_t=\partial_t^2\!\left(\frac{\beta}{\rho c^2}p^2\right),$$5 for kernels behaving like $$p_{tt}-c^2\Delta p-b\,\Delta p_t=\partial_t^2\!\left(\frac{\beta}{\rho c^2}p^2\right),$$6 near zero (Baker et al., 2022).

6. Inverse problems and imaging

The Westervelt equation has become a model for nonlinear ultrasound imaging because harmonic generation carries information about material parameters not present in linear wave data. In a lossless variable-speed model,

$$p_{tt}-c^2\Delta p-b\,\Delta p_t=\partial_t^2\!\left(\frac{\beta}{\rho c^2}p^2\right),$$7

the coefficient of nonlinearity $$p_{tt}-c^2\Delta p-b\,\Delta p_t=\partial_t^2\!\left(\frac{\beta}{\rho c^2}p^2\right),$$8 can be recovered uniquely from the nonlinear Dirichlet-to-Neumann map. The proof uses second-order linearization and Gaussian beam solutions to reduce the problem to a weighted geodesic ray transform (Acosta et al., 2021).

In the nondiffusive weakly nonlinear regime,

$$p_{tt}-c^2\Delta p-b\,\Delta p_t=\partial_t^2\!\left(\frac{\beta}{\rho c^2}p^2\right),$$9

the transmitted high-frequency wave packets acquire a measurable tilt. That tilt is proportional to the X-ray transform of the compactly supported nonlinearity $k=\beta/(\rho c^2)$0, and the second harmonic coefficient is proportional to the same line integral in the small-$k=\beta/(\rho c^2)$1 regime (Eptaminitakis et al., 2022).

Simultaneous identification problems are more demanding. For a Westervelt model with spatially varying sound speed $k=\beta/(\rho c^2)$2 and nonlinearity coefficient $k=\beta/(\rho c^2)$3, the problem can be reformulated in terms of squared slowness $k=\beta/(\rho c^2)$4 and

$k=\beta/(\rho c^2)$5

A frozen Newton method with a range-invariance identity was devised, and linearized uniqueness for simultaneous recovery was proved under suitable two-measurement conditions (Kaltenbacher et al., 2022).

The periodic setting admits even broader identification results. For the Westervelt operator

$k=\beta/(\rho c^2)$6

complex-valued time-periodic solutions excited at high frequency generate first and second harmonic boundary data. Knowledge of those two harmonic Cauchy data sets is sufficient to determine simultaneously the wave-speed coefficient $k=\beta/(\rho c^2)$7, the diffusivity $k=\beta/(\rho c^2)$8, and the nonlinearity $k=\beta/(\rho c^2)$9 (Acosta et al., 12 Sep 2025). In a related periodic framework with Robin-type boundary conditions, Fréchet differentiability of the forward operator and linearized uniqueness for an all-at-once frozen Newton scheme were established for the simultaneous reconstruction of sound speed, diffusivity, and nonlinearity from boundary measurements (Rainer et al., 15 May 2026).

These developments suggest a broad shift in Westervelt-equation research: the model is no longer used only for forward simulation of nonlinear acoustics, but also as a parameter-sensitive imaging equation in which harmonic content, boundary traces, and periodic responses encode recoverable information about the medium.