Maxwell rheology and the memory stress
The mantle is modelled as an incompressible Maxwell body: it responds elastically on short timescales and viscously on long ones, with a single Maxwell time \(\tau_M = \eta/\mu\) per layer. VILMA integrates this constitutive law explicitly in the time domain — there is no normal-mode decomposition and no Laplace transform. This is the structural choice that makes laterally varying (3D) viscosity essentially free.
The Maxwell body
For a deviatoric stress \(\boldsymbol{\tau}\) and deviatoric strain \(\boldsymbol{\varepsilon}\), the Maxwell rheology is
\[ \dot{\boldsymbol{\tau}} + \frac{\mu}{\eta}\,\boldsymbol{\tau} \;=\; 2\mu\,\dot{\boldsymbol{\varepsilon}}, \]
with \(\mu(r)\) the shear modulus and \(\eta(r)\) the viscosity. Two limits are exact and are used as validation anchors:
- Elastic / lithosphere: \(\eta\to\infty \Rightarrow\) Hooke’s law \(\boldsymbol{\tau} = 2\mu\boldsymbol{\varepsilon}\). The lithosphere is taken as \(\eta = \infty\) exactly; a large-but-finite value would seed a spurious slow mode that contaminates multi-cycle runs.
- Fluid / inviscid core: \(\mu\to 0 \Rightarrow\) the layer carries no deviatoric stress (free-slip).
Explicit memory-stress scheme
Following Martinec (2000) (§3, eqs. 23–25), the total stress at time step \(i\) is the instantaneous elastic stress — computed by the same elastic operator every step — plus a memory stress \(\boldsymbol{\tau}^{V}\) that encodes the relaxation history:
\[ \boldsymbol{\tau}^{V,i} \;=\; (1 - M)\,\boldsymbol{\tau}^{V,i-1} \;-\; 2\mu\,M\,\boldsymbol{\varepsilon}^{\,i}, \qquad M \equiv \frac{\mu\,\Delta t}{\eta}. \]
This is the \(\omega = 1\) (fully explicit) update. The memory stress enters the right-hand side of the field equations as a dissipative forcing
\[ F_{\text{diss}} \;=\; -\!\int \boldsymbol{\tau}^{V,i} : \delta\boldsymbol{\varepsilon}\ \mathrm{d}V \qquad\text{(eq. 35).} \]
The crucial consequence: the left-hand-side operator never changes in time. It is assembled, equilibrated, and LU-factored once, then reused at every step with a fresh right-hand side. Each step costs a back-substitution, not a re-factorization (see Solver). The time step enters only through \(M\), so changing \(\Delta t\) is a rescale, never a re-factorization.
Memory time schemes
The memory update is selected by scheme. All schemes share the fixed operator; they differ only in how \(\boldsymbol{\tau}^{V}\) is advanced over a step.
scheme |
Update | Order, stability | Step control |
|---|---|---|---|
"fe" (default) |
\(\boldsymbol{\tau}^{V,i} = (1-M)\,\boldsymbol{\tau}^{V,i-1} - 2\mu M\,\boldsymbol{\varepsilon}^{i}\) | 1st order, conditionally stable | equal sub-steps with \(M \le\) cfl |
"trap" |
\(\boldsymbol{\tau}^{V,i} = \dfrac{(1-M/2)\,\boldsymbol{\tau}^{V,i-1} - \mu M\,(\boldsymbol{\varepsilon}^{i-1}+\boldsymbol{\varepsilon}^{i})}{1+M/2}\) | 2nd order, A-stable (Crank–Nicolson) | adaptive, step-doubling error estimate |
"fe" is the \(\omega = 1\) scheme of Martinec (2000) and matches how VILMA-v1 marches its memory. "trap" is implicit in the end-of-step strain \(\boldsymbol{\varepsilon}^{i}\), which itself depends on \(\boldsymbol{\tau}^{V,i}\). The commit therefore iterates the step endpoint to a consistent strain–memory pair (up to max_couple_iter passes). The exponential ("etd1") and backward-Euler ("be") updates are kept in the single-degree stepper as controls for testing; the field driver advances them as "fe" and "trap" respectively.
With the explicit update, the memory stress at step \(i\) depends only on the previous step’s stress and strain — both known. Nothing in the operator depends on viscosity, so a spatially varying \(\eta(r,\theta,\phi)\) changes only the scalar \(M\) that multiplies known fields pointwise. The angular harmonics stay decoupled in the solve. A normal-mode method, by contrast, must re-solve a degree-by-degree eigenvalue problem whenever the radial structure changes and cannot represent lateral viscosity at all without abandoning the mode basis.
From layers to the same kernel in 1D and 3D
In the 1D (radially symmetric) case the memory stress evolves directly on the tensor spherical-harmonic coefficients (Martinec 2000, §9): there is no spatial grid. Per radial element it is stored for the four spheroidal tensor components \(\lambda\in\{1,2,5,6\}\) (eq. 109), and the strain coefficients follow from the nodal radial functions \(U,V\) via eq. 87–88. The dissipative right-hand side is a two-point radial Gauss quadrature of the spectral double-dot \(\sum_\lambda \lVert Z^\lambda\rVert^2\,\tau^{V,\lambda}\,\delta\varepsilon^\lambda\) with the component norms \(\{1,\ J/2,\ 2J^2,\ 2J(J-2)\}\), \(J = \ell(\ell+1)\) (eqs. 110, B13).
With laterally varying viscosity the pointwise product \(M(\theta,\varphi)\,\tau\) mixes spheroidal and toroidal stress. The toroidal components \(\lambda\in\{3,4\}\) (norms \(J/2\), \(J(J-2)/2\)) can then be carried as well (l_toroidal), and the toroidal displacement \(W\) they force feeds back into the spheroidal field. See Toroidal flow.
In the 3D (laterally varying viscosity) case the operator, the strain coefficients and the dissipative right-hand side are unchanged; only the memory advance differs. \(M = \mu\Delta t/\eta(\theta,\varphi)\) becomes a field, so for each radial element the memory and strain tensors are synthesized on the Gauss–Legendre grid from their six dyadic components (Martinec 2000, B10–B11), advanced pointwise, and analysed back to spectral space (pseudo-spectral). The per-element Maxwell routines (strain_coeffs, ve_strain_constants, dissipative_rhs, advance_memory) are shared verbatim between the 1D stepper and the field driver. There is no duplicated algorithm and no separate “3D solver”.
Three details matter in practice:
- Only genuinely 3-D elements pay for the grid round trip. An element is treated as 3-D when its lateral spread of \(\log_{10}\eta\) exceeds
visc3d_tol(default \(10^{-3}\) dex). Every other Maxwell element collapses to its lateral-mean rate and advances on the cheap per-degree spectral path, as in VILMA-v1’s 1-D/3-D layer split. A laterally uniform field therefore flags no element as 3-D and costs exactly the 1D run. - The update is analysed as an increment. The grid step is written \(\boldsymbol{\tau}^{+} = \boldsymbol{\tau} - M(\boldsymbol{\tau} + 2\mu\boldsymbol{\varepsilon})\), and only the increment is transformed back. Analysing \(\boldsymbol{\tau}^{+}\) itself would put the transform’s round-off (\(\sim 10^{-13}\) relative to \(|\boldsymbol{\tau}|\)) on a change of size \(M|\boldsymbol{\tau}|\). For the stiffest Maxwell elements (\(M\sim10^{-10}\) at \(\eta\sim10^{30}\) Pa s) that error is of order one per step.
- Elastic and fluid layers stay memory-free. A 3-D field modulates only the Maxwell layers; the elastic lithosphere of the earth model remains exactly elastic and laterally uniform, whatever the file prescribes there.
The lateral field is read as absolute \(\log_{10}\eta(\text{lon},\text{lat},r)\) (visc_3d_file, default the Bagge et al. (2021) field), bounded by visc_log10_min / visc_log10_max, and optionally perturbed by f_visc_sd standard deviations for uncertainty ensembles.
Stability
The explicit scheme is conditionally stable: \(|1-M| < 1\) bounds the step by the shortest Maxwell time present,
\[ \Delta t \;<\; \frac{2\,\eta_{\min}}{\mu}. \]
The "fe" stepper keeps a factor-2 margin: it splits each coupling interval into equal sub-steps with \(M \le\) cfl (default 1). A guard rolls a sub-step back and halves it if the memory norm grows abnormally; it normally never fires. A viscosity floor (visc_log10_min, default \(10^{19.5}\) Pa s) caps the fastest relaxation of a 3-D field and so sets a usable step. "trap" is A-stable and chooses its step by accuracy alone. Elastic layers (\(M\to 0\)) freeze and carry no memory; fluid layers (\(\mu = 0\)) carry no memory either.
vilma_viscoelastic holds both the 1D per-degree stepper (ve_degree) and the shared per-element kernel. The memory state is the only prognostic field carried across a restart (vilma_io). Time integration is validated by test_relax (a held load relaxes elastic \(\to\) fluid, \(t_{\text{relax}}\propto\eta\), step-converged) and is benchmarked in Sea-level equation.
A modal reduction was tried
For an incompressible layered Maxwell Earth the per-degree response is a finite sum of exponential relaxation modes, with only a handful of modes dominating the surface response at each degree. This suggests a natural reduction: replace the full radial memory tensor by \(K\) scalar amplitudes \(\xi_k(\ell,m)\) per coefficient, each advancing by its exact exponential. Radially varying viscosity is then exact in the limit \(K \to \infty\), and lateral viscosity is folded in as a rate modulation of each mode by the local viscosity in the depth band that mode samples.
We implemented and validated such a reduced solver, but decided against retaining it. The measured speed-up over the full explicit scheme was modest, and the approximation error grew rapidly once lateral viscosity contrasts or strongly time-varying loads were introduced. Since the explicit memory-stress scheme is already cheap — one back-substitution per step, with no angular coupling — the added complexity of maintaining a second response operator was not justified. The full solver is retained as the single response path, and the reduced solver has been removed from the code; it survives in the git history up to commit cbf77bd.