Rotational feedback (polar motion)
A redistribution of surface mass (ice, ocean water) and the resulting solid-Earth deformation change the planet’s inertia tensor, which perturbs the rotation vector — true polar wander (TPW). The shifted rotation in turn alters the centrifugal potential, itself a degree-2 surface load that drives further deformation and sea-level change. VILMA implements this two-way rotational feedback in the time domain, consistent with the explicit memory-stress machinery.
Rotation is implemented (vilma_rotation), coupled self-consistently into the sea-level equation, persisted across restarts, and tested against the Spada et al. (2011) polar-motion benchmark (test_rotation, test_rotation_sle). It is on by default (rotation = .true.) for real-Earth runs; the non-rotating community benchmarks set rotation = .false. explicitly.
Formulation
The rotational theory follows the time-domain approach of Martinec and Hagedoorn (2014), which composes with a viscoelastic, 3D-viscosity solver. On the GIA timescale the Chandler wobble is neglected, leaving a quasi-static Liouville equation for the polar-motion vector \(\mathbf{m}\),
\[ \Bigl[\,1 - \tfrac{k^{T}(t)}{k_{s}}\,\Bigr] * \mathbf{m}(t) = \Psi_{L}(t), \qquad \Psi_{L}(t) = \frac{I(t)}{C-A}, \qquad I(t) = \bigl[\delta(t) + k^{L}(t)\bigr] * I_{\rm rigid}(t), \]
where \(*\) is time convolution, \(I_{\rm rigid}\) is the rigid-body inertia perturbation from the surface load, \(k^{L}\) and \(k^{T}\) are the loading and tidal degree-2 Love numbers, \(k_s\) is the secular (fluid-limit) number, and \(C-A\) is the moment difference.
The implementation carries two degree-2 viscoelastic channels on the same per-degree operator — a loading channel forced by the rigid inertia (from the load’s (2,1) spherical-harmonic coefficient, i.e. Gauss quadrature of the actual load, so no axisymmetry is assumed) and a tidal channel forced by the centrifugal potential. Each carries its own Maxwell memory. Both channels are radially symmetric and relax on the same radial viscosity as the load response: a viscosity file’s radial profile, taking the lateral geometric mean where the file varies laterally, and the layer table only when no file is read. The polar motion \(\mathbf{m}\) is solved algebraically each step and is affine in the load, so it sits inside the SLE ↔︎ rotation fixed point safely; the channel memories are committed once with the converged \(\mathbf{m}\).
Coupling through the sea-level equation
The centrifugal potential of the perturbed rotation,
\[ \Lambda = -\Omega^{2} a^{2}\,\sin\theta\cos\theta\,(m_1\cos\phi + m_2\sin\phi), \]
is a degree-2, order-1 field. It is the first-order change of \(\tfrac12\Omega^2 a^2\sin^2\theta'\) about the displaced pole \(\hat z + m_1\hat x + m_2\hat y\): a point the pole moves toward comes closer to the spin axis, so sea level falls there. It enters the SLE as the rotational sea-surface contribution \(s_{\rm rot} = (1+k^{T})\Lambda/g - h^{T}\Lambda/g\) (geoid minus uplift). The rotation steps with the solid Earth: inside every SLE solve, each iteration’s load sets the polar motion at the end of that solver sub-step, and its \(s_{\rm rot}\) enters the same iteration’s sea surface, so the rotation ↔︎ SLE fixed point converges together with the sea level. The channels are committed where the load response’s memory advances. A sub-step longer than the channels’ explicit stability limit is sub-cycled with the load held. The rotational potential is thus treated as a load like any other — self-consistent, not a bolt-on correction.
Design constraints honoured
- Degree-2, tidal Love numbers. The rotational deformation uses tidal (not loading) Love numbers; the centrifugal forcing is purely degree 2.
- The secular number \(k_s\). By default \(k_s\) is the model’s own fluid limit \(k^T_f\): the tidal response with every Maxwell layer fluidized, which reproduces Spada test 3/2. Deriving \(k_s\) from the assumed lithosphere thickness leads to the well-known “lithosphere-thickness paradox” in rotational stability. The observed-flattening value \(k_s = 3G(C-A)/(a^5\Omega^2)\) (Mitrovica et al. 2005), recommended for deep time, is computed as well.
rotation_k_s(> 0) prescribes \(k_s\) androtation_c_minus_asets \(C-A\). - Couple through the SLE. The rotational potential enters the sea-level equation self-consistently (above).
Validation
The Liouville polar-motion solve is checked against Spada et al. (2011) test 3/2 (Spada et al. 2011) (model fluid \(k_s \approx 0.93\)–\(0.94\)), and the rotation ↔︎ SLE per-step fixed point is exercised by test_rotation_sle. See the validation ladder.
vilma_rotation holds the two degree-2 channels, the Liouville solve, and the s_rot SLE contribution; vilma_sle iterates the rotation ↔︎ SLE fixed point and vilma_timestep steps and rolls it back with the load response; vilma_io persists the polar motion and channel memory on restart. On by default (rotation = .true.); set rotation = .false. for non-rotating benchmarks.