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. VILMA2 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 admits laterally varying (3D) viscosity without changing the operator. The extra cost of a 3D run lies in the memory advance, and in the shorter time step that low viscosities require.
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 VILMA1 marches its memory. "trap" is implicit in the end-of-step strain \(\boldsymbol{\varepsilon}^{i}\), which itself depends on \(\boldsymbol{\tau}^{V,i}\). The step endpoint is therefore iterated to a consistent strain–memory pair. In the coupled model this iteration runs inside the sea-level solve, so the ocean load converges together with the memory. It stops when the surface drift changes by less than \(10^{-6}\) (relative) between passes, or after sle_max_mem_iter passes (default 20).
scheme takes "fe" or "trap"; any other value stops the run. The exponential (ETD1) and backward-Euler updates exist only in the single-degree stepper (ve_degree), as controls for the convergence tests (test_etd1, test_couple_order); they are not settings of the model.
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
visc_3d_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 VILMA1’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 of the earth model. A layer declared elastic (
rheology = 0) stays exactly elastic and laterally uniform, whatever the file prescribes there. The same holds for a fluid layer. The lid rule below is a separate, optional rule that acts on Maxwell elements.
The 3-D viscosity input chain
With l_visc_3d = .true. the lateral field is built once at initialization, in a fixed order. Each stage acts on the result of the previous one.
- Read.
visc_3d_fileholds absolute \(\log_{10}\eta\) [Pa s]. Withvisc_3d_structure = 0it is a field on (lon, lat, r), the radius in meters, with the variable namesname_visc,name_visc_lon,name_visc_latandname_visc_r. Withvisc_3d_structure\(\ge 1\) the file is in the FesmData layout,log10_visc(lon, lat, depth, structure), with depth in meters below a surface radius of 6371 km. The value is then the index of the structure to read, and thename_visc*names are ignored. The default is the Bagge et al. (2021) field of the FesmData productGLOBAL-T170_SolidEarth-Bagge2021.nc, structure 6 (v_1.0_s16). The field is interpolated in \(\log_{10}\eta\) onto the Gauss grid and the radial nodes: bilinear in longitude (periodic) and latitude, linear in radius. Outside the latitude or radius range of the file the nearest value is used. An axis may be stored in either direction. A non-monotonic axis stops the run. - Perturb (
f_visc_sd, default 0). The field is shifted byf_visc_sdstandard deviations, for uncertainty ensembles. The standard deviation is read from the variablename_visc_sdif one is named (lon-lat-r layout only). Otherwise it is relative,f_visc_rel\(\times\,\lvert\log_{10}\eta\rvert\) (default 0.1). - Shift (
visc_shift_file, default"", none). The file holds a shiftdlog10_visc(lon, lat, r)in dex, with the fixed nameslon,latandr(radius in meters). It is read and interpolated as the field itself and added to it. This changes viscosity by region and depth, for example per ensemble member. Outside the radial range of the file the shift is zero, so a file that covers only the upper mantle shifts only the upper mantle. A NaN, Inf or fill value in the viscosity, its standard deviation or the shift stops the run, naming the file and variable.visc_shift_fileneedsl_visc_3d = .true.; the run stops otherwise. - Clamp to [
visc_log10_min,visc_log10_max] (defaults 19.5 and 30.0). The floor bounds the fastest relaxation and so the time step (see Stability). - Lid rule (
visc_3d_lid_depth, default 0, off). The nodal field is first averaged to the radial elements (the mean of \(\log_{10}\eta\) at the two nodes). In a Maxwell element that lies wholly abovevisc_3d_lid_depth[m], every grid point whose \(\log_{10}\eta\) exceedsvisc_3d_lid_log10max(default 22.0) is then set tovisc_log10_max. Such points relax over millions of years, so a held-load spin-up never equilibrates them. The rule splits the lid into points that relax within thousands of years and points that do not relax on the run’s time scale. The capped points still carry memory, so they are close to elastic, not exactly elastic. The rule acts on elements, not nodes, because a node at the lid base is shared with the element below.test_visc_lidasserts that only lid elements above the threshold change. - 1-D/3-D split (
visc_3d_tol, default \(10^{-3}\) dex). Elements whose lateral spread exceeds the tolerance run on the 3-D path, as described above. A larger tolerance moves weakly varying elements to the 1-D path. This lowers the cost of the memory advance and approximates their lateral structure by its mean.
The rotation channels relax on the radial profile of the field that results from this chain.
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. vilma_earth_structure (load_visc_3d) and vilma_response build the 3-D field. The memory stress is the largest part of the state that a restart carries, but not all of it: see Memory storage and restart. 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.