Toroidal flow under lateral viscosity
A surface load on a radially symmetric Earth drives only spheroidal deformation: radial and consoidal displacement, \(U\) and \(V\). Once viscosity varies laterally this is no longer true. The Maxwell relaxation then also drives toroidal flow \(W\): horizontal, divergence-free motion with no radial component. That flow feeds back into the spheroidal field. VILMA2 carries this toroidal degree of freedom in its 3D path, as VILMA1 does.
The toroidal treatment is implemented, unit-tested (see Validation) and on by default (l_toroidal = .true.), because it changes local relative sea level substantially (below). l_toroidal = .false. is the faster option for quick runs that do not depend on local sea level near the ice sheets. With the toroidal field on, every 3-D element synthesizes, advances and analyses six tensor channels instead of four, a second per-degree drift solve is added for \(W\), and the memory state grows by about 50 %. In the lateral-viscosity benchmark (test C-i, \(\ell_{\max} = 128\), 64 threads) a run takes 1.8 times as long with the toroidal field on (18.6 h against 10.5 h to 100 kyr; see lateral viscosity). A deglaciation at \(\ell_{\max} = 170\) with 3-D viscosity takes 1.47 times as long (2.67 against 1.82 min per kyr, 64 threads, model f0c4678). With l_toroidal = .false. the model is the spheroidal-only 3D model. Radially symmetric and laterally uniform runs are unaffected either way.
Effect on a realistic deglaciation. In the last deglaciation (GLAC-1D, Bagge et al. 2021 viscosity, \(\ell_\max = 170\), model f0c4678) the toroidal field changes relative sea level, measured against the LGM reference, by 0.6 m rms over the globe at 21 ka, 2.1 m at 10 ka and 3.0 m at present day, and by up to 92 m locally (10 ka, around the margins of the Cordilleran, Laurentide and Fennoscandian ice sheets). The present-day uplift rate changes by 0.10 mm/yr rms and up to 2.6 mm/yr; the ocean-mean sea level by 0.06 m. For comparison, the whole lateral-viscosity signal (3-D minus radial mean) is 17 m rms at 21 ka. With the toroidal field on, VILMA2 also agrees about twice as well with VILMA1, which carries it: the rms difference at 21 ka falls from 4.1 to 2.2 m. Global quantities are insensitive to it; local sea level near the former ice sheets is not. The numbers are in validation/tables/toroidal.md.
Why lateral viscosity drives toroidal flow
The displacement is expanded in vector spherical harmonics (see Spectral reduction),
\[ \mathbf{u} = \sum_{jm}\Bigl[\,U_{jm}(r)\,\mathbf{S}^{(-1)}_{jm} + V_{jm}(r)\,\mathbf{S}^{(1)}_{jm} + W_{jm}(r)\,\mathbf{S}^{(0)}_{jm}\,\Bigr], \qquad \mathbf{S}^{(0)}_{jm} = \mathbf{e}_r\times\nabla_1 Y_{jm}, \]
and the strain and memory stress are expanded in the six tensor spherical harmonics \(Z^{(\lambda)}_{jm}\), \(\lambda = 1,\dots,6\) (Martinec 2000, App. B). The spheroidal field lives in \(\lambda\in\{1,2,5,6\}\) and the toroidal field in \(\lambda\in\{3,4\}\).
On a radially symmetric Earth the two families are decoupled:
- the left-hand-side operator has no term coupling \(W\) to \(U\), \(V\), \(F\) or \(\Pi\) (no pressure, self-gravity or mixed shear contribution);
- the surface forcing has no \(\delta W\) term (Martinec 2000, eq. 84), so a load never forces \(W\) directly;
- with \(M = \mu\Delta t/\eta(r)\) the memory update is diagonal in \(\lambda\), so a spheroidal memory stays spheroidal.
\(W\) is therefore identically zero, and the solver needs only the four spheroidal channels. With a laterally varying \(\eta(\theta,\varphi)\) the third point fails. The memory update \(\boldsymbol{\tau}^{+} = (1-M)\,\boldsymbol{\tau} - 2\mu M\,\boldsymbol{\varepsilon}\) multiplies the stress by a field \(M(\theta,\varphi)\) pointwise. That product mixes all six tensor components. As Martinec (2000) notes after eq. (110), lateral viscosity removes both the decoupling and the degeneracy of the spheroidal and toroidal displacements. The loop closes in three steps:
- The spheroidal memory, multiplied by \(M(\theta,\varphi)\), acquires \(\lambda = 3,4\) components.
- These force \(W\) through the dissipative right-hand side.
- The toroidal strain of \(W\), multiplied again by \(M(\theta,\varphi)\), feeds back into \(\lambda \in \{1,2,5,6\}\) and hence into uplift and geoid.
Perturbation reasoning predicts \(W\) at first order in the lateral viscosity contrast and its feedback on uplift at second order, which the tests confirm.
Toroidal flow vanishes for axisymmetric configurations (load and viscosity symmetric about a common axis), not for every configuration with a shared mirror plane. Under a reflection the toroidal potential is odd. A \(Y_{20}\) load over a \(\cos 2\varphi\) viscosity pattern is mirror-symmetric, yet drives \(W \propto \sin 2\varphi\).
Toroidal strain and memory
For \(\mathbf{u} = W(r)\,\mathbf{e}_r\times\nabla_1 Y_{jm}\) the strain is
\[ \boldsymbol{\varepsilon} = \Bigl(W' - \frac{W}{r}\Bigr) Z^{(3)}_{jm} + \frac{W}{r}\,Z^{(4)}_{jm}, \]
with the dyadic forms (Martinec 2000, B10–B11, in terms of the angular functions \(E,F,G,H\))
\[ Z^{(3)} = -F\,\mathbf{e}_{r\theta} + E\,\mathbf{e}_{r\varphi}, \qquad Z^{(4)} = G\,\mathbf{e}_{\theta\varphi} - H\,(\mathbf{e}_{\theta\theta} - \mathbf{e}_{\varphi\varphi}), \]
the toroidal companions of \(Z^{(2)}\) and \(Z^{(6)}\). Their norms (B13) are
\[ \lVert Z^{(3)}\rVert^2 = \tfrac{J}{2}, \qquad \lVert Z^{(4)}\rVert^2 = \tfrac{1}{2}J(J-2), \qquad J = j(j+1). \]
On the P1 radial elements, with \(W = W^k\psi_k + W^{k+1}\psi_{k+1}\) and \(\varepsilon = a/h + b\,\psi_k/r + c\,\psi_{k+1}/r\) (Martinec 2000, eqs. 87–88), the two strain rows the spheroidal kernel omits are
| \(\lambda\) | \(a\) | \(b\) | \(c\) |
|---|---|---|---|
| 3 | \(W^{k+1} - W^k\) | \(-W^k\) | \(-W^{k+1}\) |
| 4 | \(0\) | \(W^k\) | \(W^{k+1}\) |
The strain map is block-diagonal: \(\lambda = 3,4\) depend only on \(W\) and \(\lambda \in \{1,2,5,6\}\) only on \(U,V\). The toroidal memory is therefore carried as two extra channels of the same memory arrays (channels 5 and 6), and the spheroidal and toroidal kernels are independent. At degree 1, \(Z^{(4)}_{1m}\equiv 0\) (\(J = 2\)), so the \(\lambda = 4\) slot carries no memory there, as for \(\lambda = 6\).
The toroidal operator
The toroidal part of the shear energy (Martinec 2000, eq. 80) is
\[ 2\!\int\!\mu\,\Bigl[\lVert Z^{(3)}\rVert^2\,\varepsilon^{3}\,\delta\varepsilon^{3} + \lVert Z^{(4)}\rVert^2\,\varepsilon^{4}\,\delta\varepsilon^{4}\Bigr] r^2\,\mathrm{d}r, \]
which per element \(k\) discretizes to
\[ \mu_k\Bigl\{ J\bigl[I^1_{\alpha\beta} - I^3_{\beta\alpha} - I^3_{\alpha\beta} + I^6_{\alpha\beta}\bigr] + J(J-2)\,I^6_{\alpha\beta}\Bigr\}\,W^\alpha\,\delta W^\beta, \]
using the same Appendix-C integrals as the spheroidal block. No density, gravity or pressure enters. Rather than adding \(W\) as a fifth interleaved field to the spheroidal system, VILMA2 gives it its own per-degree operator. That operator is symmetric and tridiagonal, and it is factored once per degree and reused for all orders, loads and steps, like the spheroidal one (see Solver). The spheroidal operator is untouched.
Two conditions make it non-singular:
- Fluid regions. Where \(\mu = 0\) (the inviscid core) there is no shear energy and \(W\) has no stiffness. These nodes are pinned, \(W = 0\). This is a Dirichlet condition on a degree of freedom the fluid does not have.
- Degree 1. At \(j = 1\) a rigid rotation, \(W \propto r\), has no strain and is a null mode. It is removed by one KKT border row per solid shell, \(\sum_k W^k\!\int\!\psi_k\,r^3\,\mathrm{d}r = 0\), which states \(\int\mathbf{x}\times\mathbf{u}\,\mathrm{d}V = 0\) (no net rotation) for that shell. A structure with a solid inner core has two such shells. This is a gauge: a rigid rotation moves no output.
Coupling into the time step
Because the surface load never forces \(W\), the toroidal field has no elastic gain. It is pure relaxation drift. Within the field driver a step proceeds as follows:
begin_stepsolves, for every \((j,m)\), the spheroidal drift and the toroidal drift \(W\) from the dissipative forcing of the current memory, the latter with the toroidal operator.- The SLE fixed point calls
apply()as usual. \(W\) does not enter uplift or geoid directly, so the affine structure is unchanged. - The memory advance follows. In each genuinely 3-D element all six dyadic components of memory and strain, spheroidal and toroidal, are synthesized on the Gauss grid and multiplied pointwise by \(M(\theta,\varphi)\). The increment is analysed back into all six channels. This is where the two families exchange stress. Laterally uniform elements advance their toroidal channels on the per-degree spectral path, where they simply relax.
The toroidal channels are switched on only when some element is genuinely 3-D (its lateral \(\log_{10}\eta\) spread exceeds visc_3d_tol) and l_toroidal = .true.. The memory arrays then widen from four to six channels (about 50 % more memory state). A radially symmetric or laterally uniform run keeps four channels and pays nothing. Once on, the channels stay on for the rest of the run: a later laterally uniform field no longer drives \(W\), but existing toroidal memory must still relax.
Restart files record the channel count. A spheroidal (four-channel) restart can be read into a toroidal run, with the toroidal memory starting at zero. Reading a toroidal restart into a spheroidal run is an error, since it would discard memory.
Configuration and output
| Parameter | Group | Default | Effect |
|---|---|---|---|
l_toroidal |
&vilma |
.true. |
carry \(W\) and \(\lambda = 3,4\) once a 3-D element exists |
l_visc_3d |
&vilma |
.false. |
load a lateral viscosity field (required for any toroidal flow) |
visc_3d_tol |
&vilma |
\(10^{-3}\) dex | lateral spread above which an element is 3-D |
file_hor |
&ctl |
"" |
write the surface horizontal displacement |
The horizontal-displacement file holds east/north components on the Gauss grid of file_out: the total (u_east, u_north, spheroidal plus toroidal) and the toroidal part alone (u_east_tor, u_north_tor). The toroidal part is exactly zero whenever no toroidal field is carried.
Validation
No community benchmark exercises toroidal flow: all of them are radially symmetric or axisymmetric. The implementation is instead checked by internal consistency, symmetry and scaling tests:
| Test | Check | Result |
|---|---|---|
test_assembly |
\(W\) stiffness equals the energy rebuilt from the strain rows and norms | \(2\times10^{-16}\) |
test_tensor_sh |
six channels reproduce an arbitrary symmetric tensor field; channel cross-talk at nlat \(= 2\ell_{\max}+2\) |
\(9\times10^{-15}\) (four channels miss 35 %); \(\le 8\times10^{-13}\) |
test_sht |
toroidal synthesis is \(\mathbf{e}_r\times\nabla_1 T\) (SHTns’ native field is its negative) | round-off |
test_toroidal |
\(W = 0\) for an axisymmetric configuration | asserted: \(\lvert W\rvert/\lvert U\rvert \le 10^{-12}\) |
test_toroidal |
\(W\) driven by a \(\cos 2\varphi\) viscosity pattern; reflection selection rule | asserted: \(\lvert W\rvert > 10^{-6}\lvert U\rvert\); rule to \(10^{-10}\) (measured \(4\times10^{-16}\)) |
test_toroidal |
scaling with contrast \(\delta\): \(W \propto \delta\), uplift feedback \(\propto \delta^2\) | asserted: slopes within 0.1 of 1 and 0.2 of 2 (measured 0.987 and 1.986) |
test_toroidal |
with the toroidal field off, no \(W\) is reported | exactly zero |
test_rotinv |
coaxial cap and low-viscosity zone: \(W = 0\) on the pole; off the pole \(W/U\) at discretization level, falling with resolution | \(2.6\times10^{-4}\to1.5\times10^{-6}\) (\(\ell_{\max}\) 16 \(\to\) 32) |
test_restart |
restart with toroidal channels; four-channel file into a toroidal run | bit-for-bit continuation; spheroidal memory exact, toroidal zero |
test_toroidal is an idealized case. It holds a degree-2, order-0 load on the M3-L70-V01 Earth at \(\ell_{\max} = 12\) for 40 forward-Euler steps of 40 yr, with a viscosity contrast of 0.2 dex, and drives the response operator alone. It shows that the toroidal field is driven, has the right symmetry, scales as perturbation theory predicts and vanishes where it must. The measured values in the table are those recorded in doc/design-toroidal.md.
The full model with the toroidal field on is checked against VILMA1, which carries it, in the last deglaciation (validation/, run V7): with it on, the rms difference in relative sea level between the two codes at 21 ka falls from 4.1 to 2.2 m (area-weighted), and the present-day uplift-rate difference from 0.13 to 0.08 mm/yr. The restart tests run it through the sea-level equation and assert reproducibility. No community benchmark exercises it, and its effect on horizontal displacement has not been compared with another code.
Radially symmetric and laterally uniform runs carry no toroidal field, so their results do not depend on l_toroidal.
vilma_radial_fe (toroidal_operator: assembly, fluid-node pinning, degree-1 rotation borders), vilma_viscoelastic (strain_coeffs_tor, ve_strain_constants_tor, dissipative_rhs_tor, advance_memory_tor), vilma_tensor_sh (\(Z^{(3)}\), \(Z^{(4)}\) as channels 5 and 6), vilma_sht (toroidal vector transforms), vilma_response (the \(W\) drift, enable_toroidal, the six-channel 3-D advance), and vilma_io (restart channel count, file_hor). The design notes and derivation are in the repository’s doc/design-toroidal.md.