Solitary Wave Benchmarks in Magma Dynamics

We present a model problem for benchmarking codes that investigate magma migration in the Earth’s interior. This system retains the essential features of more sophisticated models, yet has the advantage of possessing solitary wave solutions. The existence of such exact solutions to the nonlinear problem make it an excellent benchmark problem for combinations of solver algorithms. In this work, we explore a novel algorithm for computing high quality approximations of the solitary waves and use them to benchmark a semi-Lagrangian Crank-Nicholson scheme for a finite element discretization of the time dependent problem.

💡 Research Summary

The paper introduces a benchmark problem for numerical codes that simulate magma migration in the Earth’s interior, focusing on a reduced nonlinear PDE system that retains the essential physics of more complex magma dynamics models while admitting exact solitary‑wave solutions. The governing equations (2) describe the evolution of porosity φ and compaction pressure P in a viscously deformable porous matrix. By eliminating the pressure variable, the authors obtain a single scalar equation (3) for φ with a far‑field condition φ→1 as |x|→∞. This equation possesses localized traveling‑wave solutions (solitary waves) that maintain a fixed shape while moving at a constant speed c. Such exact solutions are ideal for benchmarking because they provide a known reference state against which numerical distortion can be measured.

To compute these solitary waves with high accuracy, the authors develop a novel algorithm based on sinc‑collocation, a spectral method that expands the unknown function in shifted and scaled sinc basis functions. They first symmetrically extend the radial solution φ_c(r) to an even function on the whole real line (equation 10) so that the problem is posed on (−∞,∞). The unknown u=φ−1 is approximated by a finite sum C_M(u,h)(x)=∑_{k=−M}^{M} u_k sinc((x−kh)/h). The spacing h is chosen as a function of M (equation 15) to match the exponential decay of both the solution and its Fourier transform, guaranteeing near‑optimal convergence. Differentiation and integration operators are represented by explicit matrices D^{(l)} (equations 18–21), while the singular 1/x ∂_x operator is handled by a specially constructed matrix ˜D^{(1)} (equation 24). Substituting these discretizations into the even‑extended version of (3) yields a nonlinear algebraic system F(u)=0 (equation 25).

Because the system is highly nonlinear, a good initial guess is crucial. The authors employ a two‑stage continuation strategy. First, for the one‑dimensional case (d=1) they start from a small‑amplitude analytical approximation derived from the Korteweg‑de Vries limit (equations 26–28). They then continue in the wave speed parameter c, incrementally solving the discretized system for a sequence of c values up to the desired speed. Second, for higher dimensions (d>1) they treat the spatial dimension d as a continuation parameter, gradually increasing d from 1 to the target dimension while using the previously computed solution as the initial guess. This double continuation (in c and in d) enables robust convergence even for strongly nonlinear parameter regimes.

Numerical experiments demonstrate the method’s rapid convergence. By increasing the number of sinc nodes M, the amplitude φ_c(0) converges exponentially, with errors dropping below 10⁻⁸ for modest M (tables 1–3). The authors compute families of solitary waves for various (c, n, m) and dimensions, showing that wave amplitude generally grows with dimension, and that the dispersion relation between c and amplitude is markedly more nonlinear for (n=2, m=1) than for (n=3, m=0). They also observe a failure to converge for (m=1, d=3), suggesting limitations of the continuation approach in certain highly nonlinear regimes.

Having obtained high‑quality stationary solitary‑wave profiles, the authors use them as initial conditions for a time‑dependent benchmark. They implement a semi‑Lagrangian Crank‑Nicholson scheme combined with finite‑element spatial discretization to solve the full coupled hyperbolic‑elliptic system (34) for porosity and compaction pressure. In a moving reference frame, they simulate the off‑center collision of two 2‑D solitary waves with speeds c=5 and c=7 (Figure 3). The waves retain their shapes and propagate at the expected speeds, confirming the temporal scheme’s accuracy and stability when initialized with the precise solitary‑wave solutions.

In summary, the paper makes three principal contributions: (1) a sinc‑collocation algorithm that efficiently computes solitary‑wave solutions of the magma migration model on an unbounded domain with exponential convergence; (2) a robust double‑continuation strategy that supplies reliable initial guesses for both the wave speed and spatial dimension; and (3) a demonstration that these high‑fidelity stationary solutions serve as excellent benchmarks for time‑dependent solvers, exemplified by a semi‑Lagrangian Crank‑Nicholson finite‑element implementation. The work provides a valuable reference point for the geophysical modeling community and suggests future directions such as improving continuation techniques for challenging parameter regimes (e.g., m=1, d≥3) and extending the approach to more complex, fully coupled magma dynamics models.