# Topographic Drag at the Core-Mantle Interface

R. Monville, D. Cébron, D. Jault

► To cite this version:

R. Monville, D. Cébron, D. Jault. Topographic Drag at the Core-Mantle Interface. Journal of Geophysical Research: Solid Earth, 2025, 130 (4), 10.1029/2024JB029770. http://dx.doi.org/10.1029/2024JB029770

HAL Id: hal-05033958

https://hal.science/hal-05033958v1

Submitted on 14 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# JGR Solid Earth

## RESEARCH ARTICLE

10.1029/2024JB029770

**Key Points:**

* Topographic drag at the core-mantle boundary (CMB) is calculated using a semi-analytical model including rotation, buoyancy and magnetic effects

* The fluid waves responsible for the stress are identified, revealing the key role of Rossby waves

* The main part of the pressure torque at the CMB is due to the largest topography wavelengths

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

R. Monville,
remy.monville@univ-grenoble-alpes.fr

**Citation:**

Monville, R., Cébron, D., & Jault, D. (2025). Topographic drag at the core-mantle interface. *Journal of Geophysical Research: Solid Earth*, 130, e2024JB029770. https://doi.org/10.1029/2024JB029770

Received 24 JUN 2024

Accepted 18 MAR 2025

**Author Contributions:**

**Conceptualization:** R. Monville,

D. Cébron, D. Jault

**Formal analysis:** R. Monville, D. Cébron,

D. Jault

**Funding acquisition:** D. Cébron

**Investigation:** R. Monville, D. Cébron,

D. Jault

**Methodology:** R. Monville

**Project administration:** D. Cébron

**Software:** R. Monville

**Supervision:** D. Cébron, D. Jault

**Validation:** R. Monville, D. Cébron,

D. Jault

**Visualization:** R. Monville

**Writing – original draft:** R. Monville,

D. Cébron, D. Jault

**Writing – review & editing:** R. Monville,

D. Cébron, D. Jault

## Topographic Drag at the Core-Mantle Interface

R. Monville¹ ID, D. Cébron¹ ID, and D. Jault¹ ID

¹University Grenoble Alpes, University Savoie Mont Blanc, CNRS, IRD, University Gustave Eiffel, Grenoble, France

<figure><img src="image_2.png" /></figure>

**Abstract** The length of day variations with periods from five to one hundred years are mainly due to core-mantle interactions. Assuming a differential velocity between the core and the mantle, we investigate the pressure coupling on a core-mantle boundary (CMB) interface with topography. Including rotation, buoyancy, and magnetic effects in local models of the CMB, we provide a taxonomy of the waves radiated by the core flow along the topography. We obtain the local stress with a perturbation approach and a semi-analytical spectral model built upon these waves. We incorporate planetary curvature effects by considering a “non-traditional” β-plane approximation suited for deep fluid layers and long topography wavelengths. We calculate weakly non-linear flows and characterize the wave drag mechanism. Unlike previous works, our analysis is not restricted to strong stratification or short wavelengths. It reveals the significant impact of the Rossby waves on stress. We also show that these waves are drastically modified when considering two-dimensional topographies instead of simple ridges. For a buoyancy frequency *N* at least comparable to the rotation frequency, the main factors defining the stress are *N* and $\sqrt{U_0}$ for the small velocity amplitudes $U_0$ relevant for the Earth's core. We document the departures from this scaling law as the velocity is increased. The main part of the CMB pressure torque is due to the topography with the largest horizontal length scale. We calculate the minimum stratification for the topographic torque to produce discernible changes in the length-of-day.

**Plain Language Summary** The Earth's rotation period varies over time. Part of these variations come from the flows within its metallic liquid core. Indeed, these flows exert a drag force on the solid mantle, modifying the Earth's rotation. Like winds over mountains, the core motions generate a pressure force on the topography of the core-mantle boundary (CMB). By contrast with the atmosphere, electromagnetic effects have however to be accounted for. Relying on physical models of this drag, Earth's rotation data can provide insights into the Earth's deep interior. Previous drag models only provided estimates in certain extreme limits and for simple shapes of the core-mantle boundary (one-dimensional ridges). To unlock these limitations, our new approach combines analytical and numerical developments. The waves radiated by the flow along the topography are calculated using a local model that includes rotation, magnetic and stratification effects. We find that the largest scales of the topography contribute the most to the drag. We show that the drag is drastically changed when going beyond the ridge geometry. We identify the fluid waves responsible for the drag, revealing the key role of Rossby waves. We calculate the minimum stratification for the topographic torque to produce discernible changes in the length-of-day.

## 1. Introduction

The rotation of a planet is affected by the dynamics in its liquid core. Accurately measuring these variations is thus a way of probing its deep interior. Such rotation is accurately tracked for the Earth (Thompson et al., 2017). The core coupling with the surrounding mantle originates from viscous, pressure, gravity or electromagnetic effects (Roberts & Aurnou, 2012). Current estimates of these couplings do not account for some of the observations. For instance, the Earth's retrograde nutations record implies a strong dissipative component for the coupling at the Core-Mantle Boundary (CMB). This remains difficult to explain with our current models. Visco-magnetic (Deleplace & Cardin, 2006) and turbulent (Buffett & Christensen, 2007) effects do not seem sufficient to explain this coupling. A topographic drag has thus been also considered (Buffett, 2010). It requires high electrical conductivity of the lowermost mantle, which is possible, and strong stratification at the top of the core.

Variations in the length of the day (LOD) also illustrate the difficulties encountered by our current models of core-mantle couplings. These variations are mainly due to the different exchanges of angular momentum between the solid Earth and the fluid layers. On timescales shorter than 1 year, variations of the LOD are mainly due to the coupling with the atmosphere (Gross et al., 2004; Rosat & Gillet, 2023). The modeled contribution of the atmosphere, with the addition of less significant ocean effects, provides a satisfactory explanation for these short

© 2025. The Author(s).

This is an open access article under the
terms of the Creative Commons

Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Typical topography height as a function of wavelength, from the RMS of the observed power spectrum at the Earth's surface (solid green, Rexer & Hirt, 2015), and from a simulation (dashed green and red, Coltice et al., 2019) of the Earth's mantle convection with plate tectonics (model 5 in Arnould et al., 2020, averaged on 300–450 Myr, courtesy of T. Frasson), using the vertical stress (equilibrium with gravity). The SBF2012-T model for Earth's core-mantle boundary (CMB) (Koelemeijer, 2021), based on seismic body waves, is shown in blue, as well as the CMB flattening departure (from hydrostatics) given by studies of Earth's nutations (Dehant et al., 2022). The viscous Ekman layer thickness is 0.1–1 m at the CMB (light orange zone). The model of CMB topography used in this work is shown with the solid orange line, with a scaling in $l^2$ at small wavelength $l$ (Puica et al., 2023).</figcaption></figure>

periods. For longer periods, atmospheric or oceanic motions explain only a small part of the signal. Most of the LOD variations has been attributed to exchanges of angular momentum with the core (Finlay et al., 2023; Holme & De Viron, 2013; Jault et al., 1988). Variations in the mass distribution of the ice sheets play a role at these time scales, Kiani Shahvandi et al. (2024) quantified their contributions and showed that they are small compared to the one of the core.

The transfer of angular momentum between the core and the mantle associated with LOD has been attributed to gravitational, electromagnetic or topographic effects. The gravitational coupling has been proposed to explain the variations in LOD (Buffett, 1996; Davies et al., 2014). This gravitational torque results from a coupling between the inner core and the mantle, with a geoid deformed by mass heterogeneities. This theory is supported by the seismological inference of large-scale deep regions at the bottom of the mantle together with observations of Earth's body tides (Lau et al., 2017). The electromagnetic torque is often considered to explain some of the variations in LOD (Holme, 1998). Although the torque amplitude may be of the right order of magnitude, the LOD time series measurements cannot be recovered with the electromagnetic torque alone (Gillet et al., 2015; Schwaiger et al., 2024).

Topographic coupling has also been considered. This coupling arises when a boundary deflects a flow. In the core, it appears because of the CMB deviations from an exact sphere. This concerns not only the ellipsoidal shape of the CMB but also medium and short wavelength bumps. This has been proposed as a mechanism by numerous studies (e.g., Acheson, 1975; Braginsky, 1998;

Buffett, 2010; Glane & Buffett, 2018; Hide, 1969; Jault & Mouël, 1989; Jault, 2020; Kuang & Chao, 2001; Moffatt, 1977). The amplitude of this coupling is related to the stratification, the intensity of the magnetic field, and its orientation.

Even if it remains largely unknown, the CMB topography is expected to have a height quite similar to that of the Earth's surface, with a broad wavelength spectrum. Figure 1 illustrates the expected typical topography, including seismological constraints (for a detailed review on this topic, see Koelemeijer, 2021). In addition to seismological advances (e.g., Koroni et al., 2023; Muir et al., 2022), some insights are also provided by numerical simulation of mantle convection (Coltice et al., 2019) or by crystallisation/dissolution models (Narteau et al., 2001). Such simulations provide information about the small wavelengths (< 100 km), currently undetectable by seismology. Like at the Earth's surface (Forte & Rowley, 2022), topography is both related to lateral heterogeneity (isostasy) and to dynamical effects (through mantle viscosity and convection, Heyn et al., 2020). This topography is key to a possible pressure torque between the core and the mantle.

The torque responsible for the variation of LOD has been expected to arise primarily from the axial rotation component of the core surface flow (Hide, 1969; Jault, 2020), although the axisymmetric part of the flow is only a minor component of the total flow (Finlay et al., 2023). The axial component of the core velocity field with the largest amplitude varies over periods of the order of centuries (Pais & Jault, 2008). This is longer than the timescale of the velocity component excited by the topography (≈ 1–30 yrs for the topography considered in this article). This has led most authors, including ourselves, to consider the relative motion between the core and the mantle as stationary. The stress models on the CMB depend on the typical velocity $U_0$ of this axial component. Moffatt (1977) and Jault (2020) assumed $U_0 = 10^{-4} \text{ ms}^{-1}$, which is approximately the average zonal flow at the core surface (see Gillet et al., 2015, for estimations of core surface velocities).

Calculating the magneto-hydrodynamic (MHD) core flows in the presence of topographic effects remains a challenge, in particular within a global (near spherical) geometry. For instance, Gerick et al. (2020) investigate the core-mantle coupling in ellipsoidal geometries, but only for a mantle with infinite electrical conductivity. To avoid the difficulties associated with non-spherical computational domains, one can rather exploit the small CMB topography slope and calculate the core flows with perturbed boundary conditions (see Kuang & Chao, 2001, for such simulations in spherical geometries). To be numerically stable and converged, these simulations require
however a fluid viscosity that is much larger than the core one. Moreover, this perturbation approach is only valid when the topography height is the smallest length scale of the problem. This drastically hampers this approach by imposing the topography height to be smaller than the thickness of the viscous layer (Triana et al., 2021). Alternative models have thus been developed, relying on local models of the CMB to perform analytical perturbative developments in the inviscid regime. No viscous layer is then present, and topographic effects can be characterized in the limit appropriate for the Earth's core. However, such a local model is only relevant when the solution is localized near the CMB. Solutions confined next to the boundary can, for instance, be obtained with a horizontal magnetic field (Moffatt, 1977), or with a strong density stratification at the top of the core, as studied by Braginsky (1998), Glane and Buffett (2018) and Jault (2020).

It has long been proposed that the upper part of the core has a stratified layer (Braginsky, 1998), but its thickness and the strength of this stratification are debated. The most common explanation is the accumulation of light elements at the CMB. Alternatively, Landeau et al. (2016) found that a thick stratified layer at the top of the core may result from merging between the proto-Earth and projectile cores following a giant impact, such as the moon forming event. Bouffard et al. (2020) concluded that such a primordial layer would survive convective erosion. Insights on this subject come from seismic data (Helffrich & Kaneshima, 2010; Irving et al., 2018; Kaneshima, 2018; Van Tent et al., 2020), mineral physics, thermal evolution model, and also magnetic field observations and dynamo simulations (Gastine et al., 2020; Olson et al., 2018). Buffett et al. (2016) predict a 140 km thick layer with a buoyancy frequency (Brunt-Väisälä) close to the rotation frequency from the fluctuation of the dipole magnetic field. This stratified layer may also result from a slow cooling of the CMB. A thermally stratified layer is hard to characterize with seismic studies, since seismic velocities in liquid iron have a negligibly small temperature dependence at the Earth's core condition (Ichikawa et al., 2014). A better understanding of the topographic torque may constrain this stratification.

Building on recent works, we develop a local model of the CMB using a perturbation approach. We aim to unlock previous limitations. First, we consider more realistic 2D topographies. Second, our model is not restricted to previously studied asymptotic limits, and we investigate where the earlier models lose their validity. Third, effects from the global geometry are included through β-planes and variations of the model with latitude. Fourth, our model can go beyond the first-order solution. Our higher-order solutions provide insights into the weakly non-linear regime and the associated topographic feedback on the flow.

The article is structured as follows. We first introduce the model and the numerical method. We then study the propagation of MHD waves forced by topography and their contribution to the topographic stress. We perform a detailed comparison of the present model with previous studies. Before applying our results to the Earth's core, we finally discuss weakly non-linear effects and the remaining limitations of our model. Possible ways out are proposed as perspectives.

## 2. Model

### 2.1. Geometry and Governing Equations

We aim to build a local model approximating the region around the nearly spherical CMB. The top of the outer core is modeled as a rapidly rotating fluid, stratified in density, and permeated by a magnetic field. The fluid is assumed to be inviscid, of uniform electrical conductivity $\sigma$ and magnetic diffusivity $\eta$. Above the CMB, the mantle is assumed to be electrically insulating. Both the fluid and the solid rotate at the rotation vector $\Omega_0$, which is assumed to be constant. The fluid-solid boundary deviates slightly from a sphere of radius $R$. The spherical coordinates $(r, \theta, \phi)$ are then defined with respect to the rotation axis directed by $\Omega_0$. The spherical CMB is approximated by its tangent plane at the colatitude $\theta = \theta_0$ and longitude $\phi = \phi_0$. This allows the use of a local model and a Cartesian frame centered at $(\theta_0, \phi_0)$. The spherical basis $[1_r, 1_\theta, 1_\phi]$ is thus transformed into a Cartesian basis with $1_x = 1_\phi, 1_y = -1_\theta$ and $1_z = 1_r$ (see Figure 2). We investigate the fluid half-space $z < 0$, the solid half-space $z > 0$, and assume periodicity in the $x$ and $y$ directions (see Figure 3). The boundary topography $h(x,y)$ is described with Fourier series

$$h(x, y) = \epsilon_l \sum_j \Re[A_j e^{i(k_H \cdot r)}] \quad (1)$$
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Perturbing the global geometry (left) by small wavelength topography, a local Cartesian box model can be used (right), imposing a steady basic flow $U_0$ and a magnetic field $B_0$ in a frame rotating at $\Omega_0$, with the gravity vector $g_0 = -g_0I_z$. We denote the colatitude $\theta$, the radius $r$ and the longitude $\phi$.</figcaption></figure>

in dimensionless form. We note $\mathbf{r}$ the position vector. As the unit of length, we choose $\tilde{k}^{-1}$ related to the largest topography wavelength $2\pi\tilde{k}^{-1}$ in the direction $\mathbf{1}_x$ of the base flow (such that $x, y$ and $z$ are dimensionless). This defines the non-dimensional height $\epsilon_t = h_0\tilde{k}$ of the topography, where $h_0$ is the (dimensional) maximum topography height. In this work, we assume $\epsilon_t \ll 1$ in order to handle the topography slope as a perturbation. For each topography component, $A_j$ is the amplitude and $k_{Hj} = k_{xj}\mathbf{1}_x + k_{yj}\mathbf{1}_y$ is the horizontal wave vector. We consider the departures from a steady uniform flow $U_0$, of magnitude $U_0$. We consider that the global magnetic field $B_0$, of magnitude $B_0$ at the pole ($\theta = 0$), is steady and uniform in the local frame. The reference frame rotates at the rotation vector $\Omega_0 = \Omega_0\Omega$, with $\Omega_0$ the magnitude of $\Omega_0$. In the following, we will use local approximations of the global rotation vector (e.g., in the $\beta$-plane approximation, see Section 2.2). In the local model, $\Omega$ is then not necessarily a unit vector nor uniform, while $\Omega_0$ is constant.

The gravity acceleration $g_0 = -g_0\mathbf{1}_z$ is assumed to be steady and uniform (with the gravity magnitude $g_0$), and buoyancy effects are included under the Boussinesq approximation. We use $U_0$ as the velocity unit, then the time unit is $(U_0\tilde{k})^{-1}$. Like Jault (2020), we choose $U_0\sqrt{\rho_r\mu_0}$ as the magnetic field unit, using the fluid reference density $\rho_r$ as the density unit and taking the fluid magnetic permeability equal to the permeability $\mu_0$ of the vacuum. Using a linear equation of state for the changes in density, the dimensionless buoyancy force $g_0\rho(U_0^2\rho_r\tilde{k})^{-1}\mathbf{1}_z$, with $\rho$ the fluid density, is written as $[a_0(z) + a]\mathbf{1}_z$. The buoyancy anomaly $a$ is then an unknown variable of the problem, and the dimensionless basic state $a_0(z)$ is characterized by the local (dimensional) Brunt-Väisälä frequency $N$ through $N^2 = (U_0\tilde{k})^2\partial_z a_0(z)$, where $\partial_z$ denotes the non-dimensional $z$ derivative. Like Braginsky (1998), Glane and Buffett (2018) and Jault (2020), we consider steady perturbations of velocity $\mathbf{u}$, magnetic field $\mathbf{b}$, pressure $p$ and buoyancy $a$, governed by

$$ (\mathbf{u} \cdot \nabla)\mathbf{u} + 2Ro^{-1}\Omega \times \mathbf{u} = -\nabla p + a\mathbf{1}_z + (\nabla \times \mathbf{b}) \times \mathbf{b} \quad (2a) $$

$$ (\mathbf{u} \cdot \nabla)a = -u_z Fr^{-2}, \quad (2b) $$

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Typical topography shape handled by our code (height exaggerated). Simple 1D ridge shape of equation $h = \epsilon_t \cos x = \epsilon_t \Re(\exp(ix))$ (left) has been considered in several previous studies (Buffett, 2010; Glane & Buffett, 2018; Jault, 2020). The 2D egg-box topography $h = \epsilon_t \cos x \cos y = \epsilon_t \Re[\exp(i(x - y)) + \exp(i(x + y))]/2$ has been considered by Braginsky (1998).</figcaption></figure>
<table><caption>**Table 1**
Dimensionless Numbers and Their Typical Values for the Earth's Core-Mantle Boundary for Different Topography Wavelengths</caption><thead><tr><th>Parameters</th><th>Expression</th><th>k&#773;<sup>-1</sup> = 5 km</th><th>k&#773;<sup>-1</sup> = 100 km</th></tr></thead><tbody><tr><td>Rossby</td><td>Ro = U<sub>0</sub>k&#773;/&Omega;<sub>0</sub></td><td>2.74 &times; 10<sup>-4</sup></td><td>1.37 &times; 10<sup>-5</sup></td></tr><tr><td>Froude</td><td>Fr = U<sub>0</sub>k&#773;/N</td><td>2.74 &times; 10<sup>-4</sup></td><td>1.37 &times; 10<sup>-5</sup></td></tr><tr><td>Magnetic Reynolds</td><td>Rm = U<sub>0</sub>/(&eta;k&#773;)</td><td>6.25 &times; 10<sup>-1</sup></td><td>12.5</td></tr><tr><td>Alfv&#773;6<sup>a</sup></td><td>Al = U<sub>0</sub>&radic;(p<sub>r</sub>&mu;<sub>0</sub>/B<sub>0</sub>)</td><td>2.24 &times; 10<sup>-2</sup></td><td>2.24 &times; 10<sup>-2</sup></td></tr><tr><td>Global wavelength</td><td>&chi; = (k&#773;R)<sup>-1</sup></td><td>1.43 &times; 10<sup>-3</sup></td><td>2.86 &times; 10<sup>-2</sup></td></tr><tr><td>Lehnert</td><td>Le = Ro/Al</td><td>1.22 &times; 10<sup>-2</sup></td><td>6.12 &times; 10<sup>-4</sup></td></tr><tr><td>Elsasser</td><td>&Lambda; = RoRm/Al<sup>2</sup></td><td>3.41 &times; 10<sup>-1</sup></td><td>3.41 &times; 10<sup>-1</sup></td></tr><tr><td>Interaction parameter</td><td>Rm/Al<sup>2</sup></td><td>1.24 &times; 10<sup>3</sup></td><td>2.49 &times; 10<sup>4</sup></td></tr></tbody></table>

Note. Parameters: N/Ω₀ = 1, U₀ = 10⁻⁴ m/s, and B₀ = 5 × 10⁻⁴ T. <sup>a</sup>Al is not visible in Equation 2 as it only appears through the amplitude of the dimensionless basic magnetic field (see Equation 10).

$$ \nabla \times (\mathbf{u} \times \mathbf{b}) = -Rm^{-1} \nabla^2 \mathbf{b}, \qquad (2c) $$

with the solenoidal constraints $\nabla \cdot \mathbf{u} = \nabla \cdot \mathbf{b} = 0$, and where diffusive effects are neglected in Equations 2a and 2b. The dimensionless parameters are defined in Table 1, and p is the reduced pressure (which includes centrifugal and hydrostatic effects). Although the velocity and density unknowns are only defined in the fluid domain, the magnetic field perturbation should also be obtained in the solid. Within the insulating solid, the magnetic field **b** can be obtained from a scalar potential $\psi$ using the equations

$$ \mathbf{b} = -\nabla\psi, \quad \nabla^2\psi = 0. \qquad (3) $$

Boundary conditions complement the dynamical equations. We consider semi-infinite geometries in the z direction and impose the perturbations to vanish at infinity. The domain is periodic in the x and y directions. Since the fluid is inviscid, the flow boundary condition reduces to the non-penetration condition $\mathbf{u} \cdot \mathbf{n}|_{z=h} = 0$, with the boundary normal unit vector $\mathbf{n} = |\mathbf{1}_z - \nabla h|^{-1}(\mathbf{1}_z - \nabla h)$ directed from the fluid to the solid. Since the solid is electrically insulating and the fluid is not a perfect conductor, the magnetic field is continuous across the boundary $z = h$. Currents and buoyancy anomalies do not require specific boundary conditions, as they are fully determined by **b** and **u**.

The fluid stress on the insulating solid is due to the pressure, and the topographic stress is then given by

$$ \langle \mathbf{F} \rangle = \frac{1}{S} \int_S \mathbf{p} \mathbf{n} \, d\Sigma, \qquad (4) $$

with the elementary surface $d\Sigma = |\mathbf{1}_z - \nabla h| dS$ of the topography, and with the elementary area $dS = dx dy$ of the surface S, projected on the $z = 0$ plane, for one period of the topography. The brackets $\langle \rangle$ indicate the horizontal spatial average in x and y. In our framework, the stress of the fluid on the solid only exists in the presence of dissipative processes (here magnetic only), which occur either near the boundary or further away when waves are radiated. The diffusionless limit is singular: in the limit of vanishing dissipation, the calculated stress tends to a value that could be non-zero. This is a documented issue in the atmosphere and ocean community (e.g., Bell, 1975a), where they use the Sommerfeld (or radiation) condition (Sommerfeld, 1912) at infinity to prevent the energy propagating from the bumps from returning to the boundary. This issue is related to the geometry of our problem (semi-infinite and periodic).
## 2.2. Approximations of the Coriolis Force in the Local Frame

To include the effect of the Coriolis force in our model we need to relate the rotation of the sphere to the one of the local model. The rotation vector $\Omega = (0, \Omega_y, \Omega_z) = (0, \sin \theta, \cos \theta)$ can be simplified by considering only its vertical component since it is the most influential for horizontal flows. The Cartesian (local) coordinate y is related to the spherical (global) colatitude $\theta$ through $y = \sin(\theta_0 - \theta)/\chi$, with the parameter $\chi = (\tilde{k}R)^{-1}$ defined in Table 1. The exact spherical relation $\Omega_z = \cos \theta$ can then be expanded in Taylor series around the coordinate y = 0. As we do not use the spherical colatitude $\theta$ in the following, $\theta_0$ will be noted $\theta$ hereinafter. The series expansion of $\Omega_z$ leads to

$$ \Omega_z = \cos \theta + \chi y \sin \theta + O(y^2). \qquad (5) $$

At zeroth order, $\Omega_z = \cos \theta$, which is the f-plane approximation. Variations of the vertical component need to be included for large lengthscales. We will see below, with Equations 15–17, that it is needed for large enough wavelengths $\chi > O(Ro)$. Then, $\Omega$ is no longer uniform in the local model. When we include the first order term $\chi y \sin \theta$ of Equation 5 we obtain the $\beta$-plane approximation relevant for thin fluid layers (Rossby, 1939). For details of the $\beta$-plane approximation, see Pedlosky (2013).

In deep layers, the vertical component of the flow may be important, and it is relevant also to consider the tangential component of the rotation vector,

$$ \Omega = (\cos \theta + \chi y \sin \theta)\mathbf{1}_z + \sin \theta \mathbf{1}_y, \qquad (6) $$

which we will refer to as a tilted traditional $\beta$-plane approximation in the following.

The non-traditional $\beta$-plane includes, in addition, a linear variation $-2\chi y \cos \theta \mathbf{1}_y$ for the tangential component (Grimshaw, 1975). But it violates $\nabla \cdot \Omega = 0$, necessary for the material conservation of the potential vorticity (Dellar, 2011). To solve this problem, an additional term $2\chi z \cos \theta \mathbf{1}_z$ is then required, related to the vertical variation of the radial component of $\Omega$. This leads to

$$ \Omega = [0, \sin \theta - 2\chi y \cos \theta, (1 + 2\chi z)\cos \theta + \chi y \sin \theta]^T, \qquad (7) $$

corresponding to Equations 8.2 and 8.3 of Dellar (2011).

These various approximations of $\Omega$ introduce linear dependencies in the Cartesian coordinates within the equations. To obtain dynamical equations that are independent of the coordinates and do not involve the pressure, we consider the vorticity equation by curling Equation 2a. We thus obtain $\beta$-effect terms that are independent of the coordinates. Terms that still depend linearly on y and z are neglected since the $\beta$-plane approximation is expanded around $x = y = z = 0$. Combined with, there are 9 scalar equations for 7 scalar unknowns ($u_x, u_y, u_z, b_x, b_y, b_z$ and a). In the ridge topography case and a flow perpendicular to the ridge, we remove the equations in the direction parallel to the ridges. Otherwise, one horizontal equation is automatically satisfied if the other is, and any of them can then be removed. The pressure field is obtained a posteriori with a component of Equation 2a.

## 2.3. Taylor Expansion of the Fluid Stress on the Solid Boundary

Writing the pressure $p$ as an expansion in powers of $\epsilon_t$, as $p = \sum_{m \ge 1} \epsilon_t^m p_m(x, y, z)$, and the topography as $h(x, y) = \epsilon_t \hat{h}(x, y)$, the series expansion of $pn$ at the boundary is then

$$ pn|_{z=h} = \left[ \begin{array}{l} -\sum_{m \ge 0} \epsilon_t^m \partial_x \hat{h}(x, y) \sum_{k=0}^{m-1} \hat{h}^k(x, y) k!^{-1} \partial_z^k p_{m-k-1}(x, y, z) \bigg|_{z=0} \\ -\sum_{m \ge 0} \epsilon_t^m \partial_y \hat{h}(x, y) \sum_{k=0}^{m-1} \hat{h}^k(x, y) k!^{-1} \partial_z^k p_{m-k-1}(x, y, z) \bigg|_{z=0} \\ \qquad \sum_{m \ge 0} \epsilon_t^m \sum_{k=0}^{m} \hat{h}^k(x, y) k!^{-1} \partial_z^k p_{m-k}(x, y, z) \bigg|_{z=0} \end{array} \right] \qquad (8) $$
<figure><img src="image_6.png" /><figcaption>Figure 4. Flowchart description of the mathematical method at order 1. This corresponds to the generalized version of the method used in previous studies (e.g., Glane & Buffett, 2018; Jault, 2020).</figcaption></figure>

Ignoring the basic state $p_0$, the leading orders of this series are

$$pn|_{z=h} = \begin{bmatrix} -\epsilon_t^2 p_1(x, y, 0) \partial_x \hat{h}(x, y) + O(\epsilon_t^3) \\ -\epsilon_t^2 p_1(x, y, 0) \partial_y \hat{h}(x, y) + O(\epsilon_t^3) \\ \epsilon_t p_1(x, y, 0) + O(\epsilon_t^2) \end{bmatrix}, \quad (9)$$

where only the uniform terms contribute to the mean stress (Equation 4). This shows that the leading order of the horizontal boundary stress (Equation 4) induced by the topography perturbation is of order $\epsilon_t^2$. However, only the first-order solution of the pressure ($p_1$) is necessary to obtain the stress of order $\epsilon_t^2$. This has been used in previous studies to obtain the horizontal stress with only the linear solution (Glane & Buffett, 2018; Jault, 2020). With Equation 8 we see that this can be extended to higher orders (e.g. order 2 solution provides the tangential pressure stress at order 3).

## 3. Numerical Method

### 3.1. Semi-Analytical Approach

We aim to calculate the perturbation caused by topography on an imposed basic state that satisfies the governing equations and the boundary conditions with a flat boundary ($\epsilon_t = 0$). The unknowns are the sum of this basic state and a perturbation, such that

$$u = u_0 + \sum_{m=1}^{\infty} \epsilon_t^m u_m, \quad b = A l^{-1} b_0 + \sum_{m=1}^{\infty} \epsilon_t^m b_m, \quad a = \sum_{m=1}^{\infty} \epsilon_t^m a_m, \quad (10)$$

with the dimensionless background magnetic field $b_0 = (b_{0x}, b_{0y}, b_{0z})$, and flow $u_0 = 1_x$. The background magnetic field $b_0$ varies with latitude, but remains uniform within the local box. When not specified, we consider a dipolar magnetic field $b_0 = \sin(\theta)/21y - \cos\theta1z$. The series expansions in $\epsilon_t$ allow us to linearize the equations. Combining the expansions (10) with the dynamical equations (which include solenoidal conditions on the velocity and magnetic fields), we get a hierarchy of linear forced problems that can be solved sequentially.

We write the unknown perturbations ($u_m, b_m, a_m$) as a sum of harmonic plane waves of the form $\exp(ik \cdot r) = \exp[i(k_x x + k_y y + k_z z)]$. Similar expressions will be used to calculate a posteriori the pressure and the surface stress. As illustrated in Figure 4, injecting this ansatz into the boundary conditions and dynamical equations imposes a set of {$k_{Hj}$} = {[k_x, k_y]}$_j$, where $k_H = [k_x, k_y]$ is the horizontal part of $k$. Note that this set of {$k_{Hj}$} reduces to the topography wavevectors only at order one (when all the other physical quantities are uniform, as is the case here). For example, with a ridge topography $h = \epsilon_t \sin(k_x x)$, {$k_{Hj}$} = {[k_x, 0], [-k_x, 0]} at order 1, and {$k_{Hj}$} = {[2k_x, 0], [-2k_x, 0], [0, 0]} at order 2. For orders $m > 1$, previous orders also force a set of {$k_l$} (full wavevector including $k_z$) via the linearized “non-linear” terms of the equations (see Figure 5).

Using this ansatz in the governing equations, we obtain, at any order $m$, a linear system $\mathcal{L}(X_m) = \mathcal{N}(\sum_{l \le m-1} X_l)$ with $X_m$ the unknowns, $\mathcal{L}$ and $\mathcal{N}$ the linear operators for the current or previous orders, respectively. This system can be divided into a homogeneous system $\mathcal{L}(X_h) = 0$ and a particular one $\mathcal{L}(X_p) = \mathcal{N}(\sum_{l \le m-1} X_l)$ for orders $m \ge 2$, with $X_m = X_h + X_p$. Indeed, the forced solution $X_p$ of the dynamical equations does not a priori verify the boundary conditions, requiring the use of $X_h$ to obtain a complete solution of both the dynamical equations and the boundary conditions.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Flowchart description of the mathematical method at any order. Solving of the homogeneous part of the solution is described in detail in Figure 4.</figcaption></figure>

The homogeneous part $X_h = \epsilon_t^m \sum_j H_j e^{i(k_{Hj} \cdot r + k_{zj}z)}$ of the solution is obtained by considering each horizontal harmonic j separately. We obtain a set of linear problems of the form $A_{(k_z)}H = 0$. To ensure the existence of nonzero solutions, the determinant of A should be 0. We express this determinant as a polynomial in $k_z$ using the Berkowitz (1984) algorithm, which provides a fraction-free polynomial determinant. We numerically find each of the roots using the *polyroots* function of the arbitrary precision package “mpmath” (The mpmath development team, 2023), which is based on the *Durand-Kerner* algorithm (Kerner, 1966). For each j, we obtain a set of $k_{zjq}$.

For each $k_{zjq}$, we calculate the associated nonzero solution(s) $K_{jq}$ that are the so-called null spaces or kernels of $A_{(k_{zjq})}$, with $H_j = \sum_q C_{jq} K_{jq}$. Null spaces are found using the singular value decomposition of the matrix, which is based on the Golub and Reinsch (1971) algorithm. The eigenvalues $k_{zjq}$ can be degenerate. In such a case, they yield several solutions $K_{jq}$. In the following, we consider that the list $k_{zjq}$ contains all $k_z$, even identical ones, and their associated solution $K_{jq}$. The homogeneous solution vector is thus written as

$$X_h = \epsilon_t^m \sum_j \left[ \exp[i(k_{xj}x + k_{yj}y)] \sum_q (C_{jq}K_{jq} \exp(ik_{zjq}z)) \right]. \quad (11)$$

The particular solutions of the system are written as $X_p = \epsilon_t^m \sum_l P_l e^{i(k_l \cdot r)}$. The corresponding linear system of equations can be solved directly for each l. Finally, adding the particular and homogeneous solutions, the constants $C_{jq}$ are obtained from the boundary conditions.

The calculation involves numerous $k_z$ wave numbers, which results in very large expressions. To tackle this difficulty and the extreme values of the dimensionless parameters, we combine symbolic and arbitrarily precise
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (left) Number of required harmonics versus perturbation order *m*, for a basic flow $\mathbf{u}_0 = \mathbf{1}_x$. The topography is either a ridge ($h = \cos(k_x x)$, solid lines) or a 2D egg-box ($h = \cos(k_x x) \cos(k_y y)$, dashed). The basic magnetic field is either zero (green) or not (orange, $\mathbf{b}_0 = \mathbf{1}_z$). Horizontal dotted line: typical values for previous works (Glane & Buffett, 2018; Jault, 2020) based on linear perturbations. Horizontal dashed line: Typical calculation that can be achieved in one day on a laptop. (right) Normalized mean residual for non-penetration boundary condition versus $\epsilon_t$ for three different orders of topography perturbation *m*. Colors correspond to $N/\Omega = 20$ (blue, $Fr = 2.16 \times 10^{-5}$) and $N/\Omega = 1$ (red, $Fr = 4.32 \times 10^{-4}$). Circles: typical cases of Glane and Buffett (2018). Vertical dashed lines: critical value $\epsilon_t^c$ below which the method is convergent. Gray shading: high residual or divergent series. Parameters: $Ro = 4.32 \times 10^{-4}$, $Rm = 9.92$, $Al = 8.62 \times 10^{-2}$, $\chi = \theta = 0$, $\Omega = \mathbf{b}_0 = \mathbf{1}_z$.</figcaption></figure>

calculations, performing linear algebra automatically. In our code, written in Python, symbolic calculations are achieved using the “Sympy” package (Meurer et al., 2017). Extreme parameters can lead to badly conditioned matrices. Handling these matrices requires high-precision calculations, requiring the use of the package “mpmath”. The code is then fast enough to allow systematic exploration of parameter ranges. An order 4 calculation typically requires about 10 hours on one computing core for the simplest case of one harmonic topography (e.g., $h = \epsilon_t \sin x$).

## 3.2. Convergence Tests

Higher-order solutions are desired for their better accuracy, possibly allowing larger topography heights, but also because they can quantify the validity limits of lower-order solutions. However, the associated numerical cost can quickly become prohibitive. Figure 6a shows that the number of harmonics involved in the calculation grows exponentially with the order. Moreover, this figure also shows that the number of harmonics can already be quite large at the lowest order $m = 1$. Having optimized our implementation, our code is typically capable of calculating hydromagnetic solutions at order 3 with 2D topographies (requiring ~1000 harmonics).

We can now use these higher-order solutions to quantify convergence with the order. Figure 6b shows the residual error on the non-penetration boundary conditions for two different values of $N/\Omega_0$. We observe that the method only converges below a certain critical value $\epsilon_t^c$ of $\epsilon_t$ (shown in the figure by vertical lines), where $\epsilon_t^c$ is nearly independent of the order $m$. But the error can already be quite large for $\epsilon < \epsilon_t^c$, the case $N = 20\Omega$ of Glane and Buffett (2018) (blue point) requires an order $m = 3$ to reduce the error around 10%. These convergence issues are related to the smallest length scale of the problem (vertical dashed lines), either originating from the geometry or from the values of $k_{zjq}$. In the case considered in Figure 6, this limit is related to stratification through $\max(k_{zjq})$, which is proportional to $Fr^{-1}$ for strong stratification (Jault, 2020). When $\epsilon < \epsilon_t^c$, Figure 6 also confirms that higher orders allow the use of higher topographies for the same accuracy (e.g., one order of magnitude larger between orders $m = 1$ and $m = 3$). As allowed by our approach, considering such larger topography heights may be important to account for geophysical observations (Jault, 2020). To illustrate the accuracy gain provided by higher orders, Figure 7 compares the flow at order 1, 2 and 4 for a typical hydromagnetic calculation. The boundary condition is satisfied identically at each order, but the flow does not fully satisfy the boundary conditions of a finite-height topography. With higher orders, the flow is markedly modified near the boundary to better ensure the non-penetration boundary condition $\mathbf{u} \cdot \mathbf{n} = 0$.

Within our perturbative framework, diffusion can be restored in the momentum and buoyancy equations. This requires additional small parameters, notably for the velocity amplitude. Weakly non-linear effects in boundary
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Flow streamlines at order 1, 2 and 4 (top to bottom). Residual on the non-penetration condition $u \cdot n = 0$ is shown projected on a line a the top of each figure. Parameters: Table 1 for $\tilde{k}^{-1} = 100$ km, $\chi = \theta = 0$, $\epsilon_t = 6 \times 10^{-3}$ with 1D (ridge) topography $z = \epsilon_t \cos(x)$, and $\Omega = b_0 = 1_z$.</figcaption></figure>

layers can then be studied, with or without topography. However, this imposes the topography height to be smaller than any boundary layers. This is a severe constraint for the Earth's core, where the viscous layer is expected to be 0.1–1 m thick. This limitation is removed by considering an inviscid fluid.

# 4. Results

## 4.1. Horizontal Stress

We first map the stress as a function of two dimensionless parameters that measure the flow amplitude and the electrical conductivity (Section 4.1.1). Having found sharp transitions, we investigate what waves are present in each of these regions and which of these waves carry the stress (Section 4.1.2). The robustness of the results is investigated in Section 4.1.3, varying the topography shape and the approximation for the Coriolis term.

### 4.1.1. Survey of the Parameter Space

We have systematically surveyed the hydromagnetic topographic stress for different configurations. The leading order of the stress is $O(\epsilon_t^2)$. Theoretically, the electrical conductivity of the fluid is measured either by the interaction parameter $Rm/Al^2$ in the absence of rotation or by the Elsasser number $\Lambda = \sigma B_0^2 / (\rho_r \Omega_0)$ when the flow is rapidly rotating. We investigate how the stress component $\langle F_x \rangle$ along the flow varies with $R_0$ and the interaction parameter (Figure 8, left). We also examine the angle between the base flow and the boundary stress (Figure 8, right). The figures are shown for small constant values of the Lehnert number $L_e$ and the Froude number $F_r$ (see Table 1 for their definitions). Rotation and magnetic effects dominate in the figures' left and top parts, respectively. Examining first the amplitude of $\langle F_x \rangle$ in the figure on the left, we find various regimes and sharp transitions between these regimes. For a local Rossby number $R_0 \cos \theta < 1$ (with $\theta = \pi/4$ in Figure 8), magnetic effects are important for the stress when the local Elsasser number $\Lambda \cos \theta$, which compares the magnetic and Coriolis forces, is larger than 1 (above the dashed-dotted line). In the opposite case $R_0 \cos \theta > 1$, the parameter space is divided into two domains, on either side of $Rm/Al^2 = 1$ (black dashed line). We characterize these regions by considering the asymptotic limits $Rm/Al^2 \gg 1$ (negligible diffusion in the induction equation, i.e. ideal MHD) and $Rm/Al^2 \ll 1$ (insignificant role of the magnetic field i.e. hydrodynamic case). The lowest stress values are in the bottom middle region (D). In this area, there are large variations in stress as a function of $Rm$, and the stress vanishes in the limit of vanishing conductivity. In contrast, the stress increases with magnetic diffusion in region A. The stress is independent of the electrical conductivity in most domains (regions B, C, E). Following Braginsky (1999), we interpret this result as the signature of wave drag. Indeed, the stress is generated only by wave drag and/or dissipative processes (viscosity, ohmic dissipation...). Therefore, the stress is
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Normalized stress (left) and angle between the basic flow and the boundary stress (right) as a function of the interaction parameter $Rm/Al^2$ and $Ro$ at $\theta = \pi/4$ for the 2D egg-box topography of Figure 3, a non-traditional $\beta$-plane approximation (Equation 7) and a dipolar magnetic field. The black dashed-dotted line is $\Lambda = 1$ ($\Lambda > 1$ above this line) and the horizontal black dashed one is $Rm/Al^2 = 1$. The different visible regimes are denoted with a capital letter. Parameters: $\chi = 5.74 \times 10^{-3}$ (which is $\bar{k}^{-1} = 20$ km for the Earth's core radius), $Fr = 6.86 \times 10^{-5}$, $Le = 3.06 \times 10^{-3}$.</figcaption></figure>

generated by wave drag if the fluid is inviscid and non-conducting ($R_m \to 0$). In domain B, the average stress increases as $Ro^{-2}AlFr^{-1}$. In domain C, the stress increases as $Ro^{-1/2}Fr^{-1}$, which agrees with the results obtained by Braginsky (1998) and Jault (2020) (they found a dimensional stress scaling as $\sqrt{U_0}$ that corresponds to a non-dimensional stress scaling as $U_0^{-3/2}$). The stress in region E depends on the stratification (not shown in the figure), increasing as $Fr^{-1}$ as for stratified hydrodynamic flows (Athanassiadou, 2003). Finally, the stress projected onto the base flow is always positive. For small $Ro$, this stress is also parallel to the basic flow (see the right figure). We now investigate which types of waves are responsible for the stress and the transitions from one regime to another.

### 4.1.2. Which Waves Generate the Topographic Drag?

In Figure 9 (top parts), we reproduce the average scaled stress $\langle F_x \rangle / e_t^2$ (blue lines) for the two extreme values $Rm/Al^2 = 10^{10}$ (top row) and $Rm/Al^2 = 10^{-10}$ (bottom row) reached in Figure 8. We retrieve the abrupt stress changes already observed in Figure 8. To interpret these transitions, we also show (bottom part of each figure) the real (left panels) and imaginary (right panel) parts of the vertical wavenumbers $k_z$ of the modes entering the solution. At the order considered here, $\langle F_x \rangle$ can be obtained by summing the stress associated with each $k_z$. We show the relative contribution of each $k_z$ to the stress in color. In these figures, we also plot the theoretical dispersion relationships for a variety of waves, as given in the hydrodynamic limit (interaction parameter $Rm/Al^2 = 0$), and the diffusionless limit ($Rm^{-1} = 0$) by Finlay (2008) and Salhi et al. (2017). The wave frequency, $U_0\bar{k}$, determined by the motion of the fluid relative to the topography, corresponds to the dimensionless frequency 1 in the fluid frame.

We first investigate the hydrodynamic limit (top figures). Inertial waves are waves that propagate due to the Coriolis restoring force. When inertial waves are modified by buoyancy force they become inertia-gravity waves, also named inertial-gravity waves. These waves exist only in a specific frequency range (Bell, 1975b). This domain of existence (red vertical dotted lines), calculated here in the $f$-plane approximation, is

$$1 < Ro_l = \frac{Ro}{2 \cos \theta}, \qquad (12)$$

for $Fr \ll 1$, as used in Figure 9 (see the derivation of Equation 12 in Appendix B1). These waves (dashed blue lines) reduce to internal gravity waves (dashed-dotted red lines) with dispersion relationship given by
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Normalized stress (blue solid line at the top of each figure) as a function of Ro and associated absolute value of the real (left) or imaginary (right) part of vertical wavenumbers $k_z$. Vertical wave numbers are color-coded (red to yellow) to show the relative contribution of each $k_z$ to the total stress. Top figures (a, b) are calculated in the hydrodynamic limit with $RmAl^{-2} = 10^{-10}$. The bottom ones (c,d) are in the ideal MHD limit $RmAl^{-2} = 10^{10}$. Other parameters are the same as Figure 8 ($\chi = 5.74 \times 10^{-3}$, $Fr = 6.86 \times 10^{-5}$, $Le = 3.06 \times 10^{-3}$, $\theta = \pi/4$, $h = \epsilon_t \cos x \cos y$, non-traditional $\beta$-plane approximation (Equation 7) and dipolar magnetic field). Other lines show diffusionless MHD wave dispersion relationships (see Appendix A or Salhi et al., 2017). The blue vertical lines correspond to the theoretical limits for internal Rossby waves (Equation 15, dotted) and MAC-Rossby waves (Equation 17, dashed, with $\tau = 9$). The red vertical dotted line is Equation 12 and the dashed one is Equation 16.</figcaption></figure>

Equation A2 for $Ro \gg 1$ (see the dispersion relationship of inertia gravity waves in Appendix B2). Propagation of internal waves, also known as internal gravity waves, is made possible by the gravitational restoring force. When generated by flow over topography, they are named Lee waves or mountain waves. Outside their domain of existence (left of the red vertical dotted line, $Ro \lesssim 2^{1/2}$), they become evanescent (large $|\tilde{\nu}(k_z)|$ in Figure b) and the average stress drops suddenly to diffusion-dependent values. The stress varies proportionally to $1/|\tilde{\nu}(k_z)|$ in these zones (where $k_z$ is the wave number that contributes the most to the stress). In the f-plane approximation, no other stress transition occurs. But including the $\beta$-effect, as in Figures 8 and 9, internal Rossby waves (solid black lines in Figure 9a) differ from inertia-gravity waves at low $Ro$, increasing the average stress drastically. These waves are Rossby waves (waves generated by the variation of angle between the boundary and the rotation vector
with latitude) modified by the buoyancy force. These internal Rossby waves only exist between two values of Ro, or below a certain value of Fr, which are all given by the roots of

$$ (k_H^2 - A^2 Fr^2) Ro^2 + 2 \sin \theta (B Fr^2 - 1) \chi Ro + C Fr^2 \chi^2 = 0, \quad (13) $$

where A, B and C, given in the appendix B3, are only function of θ. The bound in Fr is then

$$ Fr^2 < \frac{(R_{0k})^2 - 2\chi R_0 \sin \theta}{(R_0 A)^2 - 2\chi R_0 B \sin \theta - \chi^2 C}, \quad (14) $$

which reduces to $Fr < k_H A^{-1}$ for $\chi \ll 1$, and further simplifies into $Fr < 1$ for ridge geometries ($k_H = 1$, $k_y = 0$). The two bounds in Ro can also be obtained exactly, giving

$$ 8\chi Fr^2 \cos \theta \cot \theta < Ro < 2\chi k_H^{-2} \sin \theta \quad (15) $$

for $Fr \ll 1$. The lower bound, related to the $\chi^2$ term $CFr^2\chi^2$ in Equation 13, does not exist when considering Equation 6, and requires the more sophisticated approximation (Equation 7). This transition occurs at $Ro \approx 10^{-10}$ in Figure 8 or 9, and is thus not visible. The upper bound (dotted blue line) is at $Ro \approx 4 \times 10^{-3}$.

In the MHD case (bottom figures), internal waves, which are almost diffusionless, also exist in a specific range of frequency. This domain of existence depends on Ro, but also on Al. It is given by (vertical red-dashed lines, $Al \approx 2.915 \times 10^4$)

$$ Al > 2k_H Fr^{-1} \cos \theta, \quad (16) $$

for a dipolar magnetic field $b_0$, and in the relevant limits $Rm \gg 1$, $Fr \ll 1$ and $Ro \gg 1$ considered in Figures 9c and 9d. The stress drops suddenly at smaller $Al$ (left of the vertical red-dashed lines). Similarly to the hydrodynamic case, this happens at a transition from large to smaller $Ro$. However, the stress remains here independent of diffusion. This originates from MAC (Magneto-Archimedes-Coriolis, dashed green) waves, that arise from magnetic and rotation effects. Alfvén waves (dashed purple) contribute only marginally to the stress. In the f-plane approximation, no other stress transition occurs. Including $\beta$-effects, MAC Rossby waves (solid black lines in Figures c and d) differ from MAC waves in a finite range $Ro$. For a dipolar magnetic field $b_0$, this domain is bounded by the following approximate formula (Appendix B)

$$ k_H^2 Ro = \chi \sin \theta \left[ 1 \pm \sqrt{1 - (Frk_H)^2 \left[ \omega_1 (k_y/\chi)^2 + \omega_2 \cot^2 \theta \right]} \right], \quad (17) $$

for $Fr, Ro \ll 1$, where $(\omega_1, \omega_2)$ are (1,0), (9,0) and (9,16) when considering the three Coriolis force approximations given by Equations 5–7, respectively. Equation 17 illustrates the importance of 2D topographies ($k_y \neq 0$) and density stratification ($Fr \neq 0$) when considering $\beta$-plane effects ($\chi \ll 1$ but $\chi \neq 0$). In Figures 9c and 9d, this corresponds to the zone between the two vertical dashed blue lines with $2.6 \times 10^{-6} \lesssim Ro \lesssim 4 \times 10^{-3}$. A branch of MAC Rossby waves differs from MAC waves by its imaginary part (solid black line in Figure d). These waves participate in the boundary conditions. As a result, the MAC waves that carry the stress have a lower amplitude, reducing the total stress drastically. Stratification turns out to be key in all stress mechanisms investigated here.

Large stress values are obtained at the frontiers of Rossby and Rossby-MAC wave domains. Calculating the eigenmodes, we obtain free steady (MAC-) Rossby modes that can be excited by the topography-forced waves at these frontiers. Looking at the associated wavenumber, at resonances, the wavevector is found to be perpendicular to the rotation vector in the hydrodynamic case, and perpendicular to the magnetic field in the ideal MHD limit (which we exploited to find Equation 17).
<figure><img src="image_12.png" /></figure>

**Figure 10.** Normalized stress (colors) as a function of interaction parameter *Rm*/*Al*² and *Ro* at θ = π/4 for (left) a ridge *h* = *ε*ᵢ cos *x*, Ω = cos θ*1*, (f-plane) and *b*₀ = −cos θ*1*. The right side is with a 2D topography *h* = *ε*ᵢ cos *x* cos *y*, a non-traditional *β*-plane approximation (Equation 7) and a dipolar magnetic field. The stress of Figure (a, b) is estimated with *χ* = 0.0014, that is, *k*⁻¹ = 5 km for the core. For (c, d) the stress is shown for *χ* = 0.0287, that is, *k*⁻¹ = 100 km for the core. The star represents the Earth's core (purple). The bounds for the internal waves are the red dotted line (Equation 12) and the dashed line (Equation 16). The limits for Rossby waves are the blue dotted line (Equation 15) and the dashed blue ones (Equation 17). The dashed-dotted black line is Λ cos *θ* = 1. The type of wave that contributes the most to the total stress is written within each zone. Parameters: (top) *Le* = 1.22 × 10⁻² (bottom) *Le* = 6.12 × 10⁻⁴, (see Table 1 for other parameters at the two wavelengths).

4.1.3. Influence of Topography Shapes and β-Effects

Using Figure 10, we investigate the combined role of topography wavelength and β-effect. The color map shows
the normalized average stress for two approximations of the Coriolis force (f-plane on the left and non-traditional
β effect on the right) and two different topography wavelengths (k⁻¹ = 5 km, top and k⁻¹ = 100 km, bottom).
Note that Figure 8 presents an intermediate case between Figures 10b and 10d. The region A, where Rossby-MAC
waves propagate, expands with the topography wavelength, as shown by Equation 17. In the limit of large
topography wavelength χ ≫ Fr, this domain covers the region delimited by Equation 15 and Λ > 1. For the 1D
<table><caption>**Table 2**
Approximations Used in Previous Studies, All Relaxed by Our Code</caption><thead><tr><th></th><th>Braginsky (1998)</th><th>Glane and Buffett (2018)</th><th>Jault (2020)</th></tr></thead><tbody><tr><td>Advection</td><td>×</td><td>✓</td><td>✓</td></tr><tr><td>2D bumps</td><td>×</td><td>×</td><td>x<sup>a</sup></td></tr><tr><td>β-plane</td><td>✓<sup>b</sup></td><td>×</td><td>✓<sup>b</sup></td></tr><tr><td>Any tilt of <b>B</b> and Ω</td><td>×</td><td>×</td><td>✓<sup>c</sup></td></tr><tr><td>Asymptotic limits</td><td>Λ &lt;&lt; 1</td><td></td><td>Fr &lt;&lt; 1, Ro &lt;&lt; 1<br>RmFrAl<sup>-1</sup> &lt;&lt; 1<sup>d</sup></td></tr></tbody></table>

ridge topography with β-effects, studied by Jault (2020), the region A extends to infinitely small Ro, and the region B at small Ro is then inexistent (not shown). This can be inferred from Equation 17.

We observe dramatic differences in the stress value between the f-plane and the β-plane approximations, especially for large wavelength topographies (bottom figures). At small Rossby numbers, estimates of the stress calculated with the two approximations remain different for small wavelengths of the topography too.

## 4.2. Comparison With Previous Studies

The topographic stress has been studied before, using various approximation and asymptotic methods (see Table 2). These previous studies have allowed us to validate our code. We have first reproduced the hydrodynamic stress on a 2D-topography of Equation 1 due to inertia gravity waves

$$ \langle F_x \rangle = \epsilon_t^2 \sum_j \left| \text{Re} \frac{k_{xj}}{k_{Hj}} A_j^2 \sqrt{(Fr^{-2} - 1)(1 - R_o^{-2})} \right| \quad (18) $$

Equation 18 is equivalent to Equation 20 of Legg (2021), see also Bell (1975b) and Nikurashin and Ferrari (2010). Then, we have successfully validated our hydromagnetic stresses against previous results, obtained for 1D topographies.

Figure 11 (left) compares stress values from our work and other studies (Braginsky, 1998; Glane & Buffett, 2018; Jault, 2020). We take the β-plane approximation and consider a 1D ridge. Results are shown for a large

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Average stress as a function of Rm for Fr = 1.37 × 10⁻⁵ (left) and as a function of Fr for Rm = 12.5 (right). Ridge topography (k_y = 0). We compare our model (solid black) with those of Jault (2020) (dotted blue) and Braginsky (1998) (dashed green), which both include β-plane effects. Vertical lines in the left figure are Λ = 1 (green dashed) and Rm = Al/Fr (blue dotted). The scaling at RmFr/Al ≫ 1 is also indicated (purple). Vertical dashed lines in the right figure show Fr = Le/(2 cos θ) in purple (Equation 20) and Fr = 1 in blue (Equation 14). The waves generating the stress in each domain are given. Parameters: b₀ = 1, θ = π/4, with a traditional β-plane approximation. Others parameters are given by Table 1 for k̃⁻¹ = 100 km.</figcaption></figure>

<footnote><sup>a</sup>Preliminary results (see Appendix C of Jault (2020)). <sup>b</sup>Ω is present only via its vertical component. <sup>c</sup>Only two orientations of <b>B</b> (vertical and horizontal). <sup>d</sup>Limit of the analytical expression after his Equation 48.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Topographic mean stress. (left) As a function of $\chi$, for different approximations of the Coriolis force: f-plane (dotted blue), traditional $\beta$-plane (dashed green), tilted and uniform (dashed-dotted red), tilted with $\beta$-effects on z-component (solid black). Note that values obtained with non-traditional $\beta$-plane (Equation 7) are superimposed to this latter case (Equation 6). The green and black vertical dashed lines are Equation 17 for ($\omega_1, \omega_2$) = (1,0) and (9,0) respectively. The dotted vertical black line is $\chi = \Lambda$. (right) As a function of $Fr$, considering only the $\epsilon_i^2$ (teal) or the $\epsilon_i^4$ (red) contribution to the stress. Scalings for small $Fr$ are noted in the figure. Parameters: Table 1 for $\bar{k}^{-1} = 100$ km, $h = \epsilon_i \cos x \cos y$, $\theta = \pi/4$.</figcaption></figure>

topography wavelength of $2\pi\bar{k}^{-1} \sim 600$ km. All the models agree for moderate $Rm$ values and Earth-like values of the other parameters. Previous investigations focused on Elsasser numbers $\Lambda$ differing at most by two orders of magnitude from 1. Jault (2020) extended the study of Braginsky (1998), which was limited to $\Lambda < 1$ (see his Figure 3). He thought his results to be valid when $RmFr^2/Al^2 = \Lambda^{-1}Fr^2/Ro \ll 1$ and obtained analytical expressions for the stress (his equation B14) assuming that the stress is mostly carried by only one mode, the mode of largest $|k_z|$. It turns out that the wavenumber modulus $|k_z|$ of a second mode becomes as large as the wavenumber modulus of the first one when $RmFr/Al$ reaches a value of order unity. Then both modes are responsible for the stress and the calculation of Jault (2020) is not valid anymore: it requires indeed $RmFr/Al \ll 1$ (see the vertical blue dotted line in Figure 11). For completeness, we derive the asymptotic expression for the stress in the opposite limit $RmFr/Al \gg 1$. The stress arising from the two modes with the largest moduli cancels out at the zeroth order in the small parameter $Al/RmFr$ and we find

$$ \langle F_x \rangle = \frac{\sqrt{Al}}{8\sqrt{2} Fr^{5/2} Rm} \epsilon_i^2, \qquad (19) $$

as detailed in the Supplemental Material.

We also examine the stress dependency on $Fr$ (Figure 11, right figure). Both Braginsky (1998) and Jault (2020) found a stress varying as $Fr^{-1}$. In the considered limit $Ro \ll 1$, all calculations agree for

$$ Fr < Le(2 \cos \theta)^{-1}, \qquad (20) $$

where the stress is due to MAC waves. For $(2 \cos \theta)^{-1}Le < Fr < 1$, the stress is due to internal Rossby waves (see Appendix C). It however keeps its scaling as $Fr^{-1}$. We observe an increase in the stress at the transition $Fr = Le/(2 \cos \theta)$. The stress drops abruptly for $Fr > 1$ and becomes independent of $Fr$. The residual stress is then associated with evanescent (heavily damped) waves.

We can then compare the traditional $\beta$-plane approximation with more sophisticated approximations of the Coriolis force. In Figure 12 (left) we show the normalized stress as a function of the amplitude $\chi$ of the $\beta$ effect, for Earth's core parameters. Stress is given for four different approximations of the Coriolis force, from the simplest f-plane to the non-traditional $\beta$-plane. Tilting the rotation vector has a major impact on the stress when $\chi$ is smaller than its value from Equation 17, see differences between blue/green and red/black curves. We find that the results obtained with the non-traditional $\beta$-plane (Equation 7, not shown) are very close to those obtained with traditional tilted $\beta$-plane (Equation 6). The additional terms introduced by the non-traditional $\beta$-plane, which involve z variations of $\Omega$, play a significant role when the perturbations propagate far from the boundary. This is, for instance, the case for Quasi-Geostrophic flows. This figure displays the case where the flow is dominated by
<figure><img src="image_15.png" /><figcaption>**Figure 13.** (left) Mean flow of the order $m = 2$ perturbation as a function of depth. The blue line is $u_x$ and the orange one is $u_y$. (right) Typical decay length for the mean flow of the order $m = 2$ (dashed-dotted black) compared to the smallest wavelength (dashed black) as a function of $Le$. The red solid line is the order 2 mean flow at the boundary. The vertical black dashed line is the Earth's $Le$ used in the left figure. The green vertical dashed line is $Le = Fr$. Other parameters: Table 1 for $\bar{k}^{-1} = 100$ km, $\theta = 0$, with an egg-box topography.</figcaption></figure>

stratification effects ($N/\Omega_0 = 1$) while rotation effects dominate if $N/\Omega_0 \ll 1$. Because we have found the pressure stress to be small for Quasi-Geostrophic flows, the new terms introduced in the non-traditional approximation have no significant impact on our study.

In Figure 12 (right), we separate the contribution of higher orders up to $O(\epsilon_t^4)$ to the stress. The $O(\epsilon_t^3)$ contribution vanishes because of the monochromatic topography (like any odd degree contribution). Non zero average $O(\epsilon_t^3)$ stress requires a topography with three (non necessarily distinct) wavenumbers $k_1, k_2$ and $k_3$ such that $k_1 \pm k_2 \pm k_3 = 0$. This happens when the topography has a continuous spectrum. Order 4 contributions mostly enhance the stress, but also reduce it when the stress at order 2 presents a local maximum. High orders seem important to better describe the singularities. Order 2 and 4 contributions to the stress scale respectively as $Fr^{-1}$ and $Fr^{-3}$ for $Fr \ll 1$, and become independent of $Fr$ for large $Fr$.

### 4.3. Weakly Non-Linear Feedback on the Basic Flow

The non-linear interaction between the waves driven by the imposed differential velocity $u_0 = 1_x$ yields a horizontally averaged flow $u_2$ that opposes $u_0$ at the solid-mantle interface (Figure 13). This agrees with the orientation of the stress on the solid, which is almost parallel to $u_0$ (see Figure 8 right). Interestingly, the flow $u_2$ is not confined to a thin sheet close to the boundary. It keeps instead a significant amplitude in a layer of thickness comparable to the topography wavelength $k_x^{-1} = 1$. However, there is a weaker variation of $u_2$ on a much smaller vertical scale next to the boundary. As $\epsilon_t$ is increased to $O(10^{-1})$, the flows $u_2$ and $u_0$ have similar strengths, invalidating the starting hypothesis of the calculation. This puts a limit on $\epsilon_t$ that nicely agrees with the one shown in Figure 6 (right, see the dashed red line, for $m = 2$ and $N/\Omega = 1$ as in this subsection).

Figure 13 (right) shows that our conclusions do not depend strongly on $Le$. The typical vertical length scale of $u_2$ remains of the order $k_x^{-1} = 1$. Similarly, the relative amplitudes of $u_2$ and $u_0$ at the boundary do not vary much with $Le$, except for a transition near $Le = Fr$.

Beyond these first results on the topographic feedback on the flows, we would like to complete this study with a full description of the shear layer attached to the boundary. However, the perturbation method does not allow the description of the full feedback on the basic state, requiring new approaches (see the discussion in Section 5).

### 4.4. Mean Electric Currents and Flows Across the Local Box

Our periodic model is not isolated. Consequently, the total forces acting on the fluid are not the opposite of the forces acting on the solid. Instead, they are exactly zero since we assume that there is no fluid acceleration in our model. The same kind of reasoning holds when considering only the magnetic forces. The total magnetic force on the solid is zero since we treat the solid domain as an insulator. However, the flux of electrical currents across the fluid domain yields a non-zero magnetic force on the fluid. Similarly, the Coriolis force integrated into the fluid
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Magnitude of Lorentz (blue/purple), Coriolis (yellow/green) and $(u \cdot \nabla)u$ (orange/red) forces as a function of $Ro$, calculated between the bumps (see Equation 21, solid lines) and in the fluid interior (dashed). The black line is the boundary pressure stress. The solid vertical line is $\Lambda = 1$, the dashed blue line is $Le = 1$. The parameters are those of Table 1 for $\tilde{k}^{-1} = 100$ km. The calculations are for a ridge topography with $\Omega = b_0 = 1_z$.</figcaption></figure>

domain would be zero if this domain was closed. Instead, the flux of mass yields a non-zero integral of the Coriolis force in the fluid domain.

The linear solution induces a second-order volumetric flow rate and an electric current between the bumps (in the volume bounded by $z = h$ and $z = -e_r$). Looking at the plane $y = 0$, at leading order, these fluxes can be written as

$$V = \int_S u_y(x, z) \, dS = \int_S h(x)u_y(x, 0) \, dx, \quad J = \int_S j_y(x, z) \, dS = \int_S h(x)j_y(x, 0) \, dx, \quad (21)$$

where $S$ is a plane of constant $y$, and $j = \nabla \times b$ the electrical current density. These fluxes yield mean Coriolis $2Ro^{-1}\Omega_z V$ and Lorentz forces $b_0 J$ (Jault, 2020). Since we calculate second-order solutions, we can obtain the forces that correspond to the volumetric flow rate and the current in the fluid interior (between $z = 0$ and $z = -\infty$) and compare them with those calculated between the bumps (noted respectively interior and boundary in Figure 14). First, the Lorentz, Coriolis, and inertial $(u \cdot \nabla)u$ forces in the interior cancel out. Second, the inertial contribution at the surface exactly balances its contribution in the interior as a consequence of the non-penetration boundary condition and incompressibility. For an Elsasser number $\Lambda < 1$, the Coriolis force is stronger than the Lorentz force and vice versa for $\Lambda > 1$. In both cases, the force calculated between the bumps predominates.

This discussion means that it does not suffice to stitch together local models to obtain a global model. The total fluxes of momentum and electrical current vary from one local box to the next (Figures 13 and 14). This problem is highlighted if we consider the azimuthal average of the fluxes. The global torque about the rotation axis is obtained from the azimuthal average drag integrated over the core surface. At each latitude, the azimuthally averaged drag is a function of the flux of mass and electrical current in the meridional direction. At the poles, these fluxes vanish and we cannot reconcile the zero flux at the poles with non-zero fluxes at other latitudes. We can instead introduce the local calculation of the topographic stress on the mantle in a global model developed independently.

## 5. Discussion and Conclusion

We now discuss the geophysical consequences of our study. Our systematic study allows us to identify the various regimes, which is useful as the Earth lies near the transition between different domains (Figure 10). Figure 15 (left) shows the dimensional value of the stress as a function of the basic flow amplitude $U_0$. For small Earth-like values of $U_0$, the dimensional stress scales as $U_0^{1/2}$ as previously found by Braginsky (1998), Glane and Buffett (2018) and Jault (2020). This regime ends with a sharp drop in the shear value at $U_0 \sim 10^{-1}$ m.s$^{-1}$ for $\tilde{k}^{-1} = 100$ km and $U_0 \sim 5 \times 10^{-3}$ m.s$^{-1}$ for $\tilde{k}^{-1} = 5$ km. In Appendix D, we detail our reasoning to characterize the waves that carry the stress, for $\tilde{k}^{-1} = 5$ km. For small values of $U_0$, MAC waves are responsible for
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Average stress (left) as a function of basic flow velocity for a unit height, $\tilde{k}^{-1} = 100$ km (orange) and $\tilde{k}^{-1} = 5$ km (green), or (right) as a function of topography wavelength $l = 2\pi\tilde{k}^{-1} = 2\pi\chi R$ for the Earth's core parameters, a unit height (blue), or our core-mantle boundary height model (red) shown in Figure 1. An egg-box topography with a non-traditional $\beta$-plane approximation and a dipolar magnetic field is considered in both cases. Green and orange vertical dashed lines are given by the limits of Inertia-gravity waves for corresponding values of $\tilde{k}$ (see Equation B5). Equation 17 gives the blue vertical dashed line. Red horizontal lines are the length of the day derived mean stress values $2.4 \times 10^{-3}$ Nm$^{-2}$ and $2.7 \times 10^{-2}$ Nm$^{-2}$, from Roberts and Aurnou (2012) and Glane and Buffett (2018), respectively. Parameters: $\rho_r = 10^4$ kg.m$^{-3}$, $B = 5 \times 10^{-4}$ T, $\eta = 0.8$ m$^2$.s$^{-1}$, $R = 3486$ km, $\Omega_0 = N = 7.29 \times 10^{-5}$ rad.s$^{-1}$, $\theta = \pi/4$ (and $U_0 = 10^{-4}$ m.s$^{-1}$ for right figure).</figcaption></figure>

most of the stress, although other waves, especially Rossby-MAC waves, participate in the boundary conditions and are a key component of the solution. The drop in stress value occurs when evanescent gravity waves replace MAC waves as the main carrier of the stress. For large velocity amplitudes, we do not find a simple scaling of $\langle \tilde{F}_x \rangle$ with $U_0$. A scaling in $U_0^{-1/2}$ is nevertheless found for very large $U_0$.

For a unit height of the topography, Figure 15 (right) shows that the topographic mean stress varies over several decades in a non-trivial way. For all wavelengths and $U_0 = 10^{-4}$ m.s$^{-1}$ the stress is carried by MAC waves. The sharp variation at 30 km corresponds to the emergence of Rossby-MAC waves, as discussed in the paragraph above. This transition is given by Equation 17 (blue dashed vertical line). This (scaled) stress keeps the same order of magnitude below 50 km, before decreasing as $l^{-1/2}$ for larger wavelengths. Our results can be combined with estimates of Earth's core topography height from Figure 1 to provide the average stress (red curve). The increase in topography height with wavelength $l$, in $l^2$ for $l < 10^3$ km, more than compensates for the decrease in scaled stress. In conclusion, the most important topography components in the CMB for pressure coupling are those with the largest wavelengths. This shows the importance of properly accounting for the global geometry through the $\beta$-plane.

Figure 15 (right) shows that the topographic drag can explain the mean stress derived from LOD observations, for wavelengths larger than $10^3$ km when using a buoyancy frequency $N/\Omega_0 = 1$ and the CMB height model of Figure 1. Conversely, one could rather use this mean stress to constrain the stratification. In Figure 16, we vary the buoyancy frequency and the topography wavelength, showing in colors the topography height required to explain the LOD observation. It reveals that a non-zero stratification is necessary to retrieve the desired mean stress, requiring buoyancy frequencies $N/\Omega_0 \ge 0.1-1$. Note that for most of our applications perturbation is restricted to a thin layer (i.e., large $k_z$). Moreover, the stress is carried out mostly by the largest $k_z$. This allows us to model our stratified layer using an infinite thickness layer. This approximation should be questioned in the case of stratification restricted to a thin layer or small $k_z$.

For typical CMB values, the topographic stress obeys a simple function of $N$, $\Omega_0$, $l$ and $h_0$. In a region encompassing the red curve in Figure 16, we find (using $l = 2\pi\tilde{k}^{-1}$)

$$ \langle \tilde{F}_x \rangle = \rho_r F(\theta) \frac{\sqrt{\rho_r \mu}}{B_0} N \Omega_0 \sqrt{\eta U_0 / l h_0^2}, \quad (22) $$
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Topography height $h_0$ (color) necessary to reproduce the mean stress derived from the length of the day variations (typical value of $2.4 \times 10^{-3}$ Nm$^{-2}$, Glane & Buffett, 2018), varying the topography wavelength $l$ and the Brunt-Väisälä frequency $N$. The solid red line shows the height $h_0$ of our core-mantle boundary height model shown in Figure 1. The dashed blue line shows Equation 17. Parameters: $\theta = \pi/4$, $\rho_r = 10^4$ kg.m$^{-3}$, $B = 5 \times 10^{-4}$ T, $\eta = 0.8$ m$^2$.s$^{-1}$, $R = 3486$ km, $\Omega_0 = 7.29 \times 10^{-5}$ rad.s$^{-1}$, $U_0 = 10^{-4}$ m.s$^{-1}$, and $\Omega$ from Equation 7.</figcaption></figure>

with the function $F = O(1)$ of latitude $\theta$ displayed in the Supplemental Material. This function is nearly constant for $\theta \in [0, \pi/3]$, where it reaches its maximum $\max(F) \sim 0.7$. Interestingly, Equation 22 can be retrieved as an extension of Jault (2020) to the egg-box geometry (see details in the Supplemental Material). It is specifically obtained in the regime $\Lambda \lesssim 1$, where the study of Braginsky (1998) is also valid. This gives analytically the prefactor as $\sqrt{2\pi}/4 \approx 0.63$, in perfect agreement with our numerical results. As noted by Jault (2020), the stress is nearly independent of $\Lambda$ in the geophysical regime $\Lambda \sim 1$. Therefore, with Earth-like parameters, the main factors defining the stress are $N$ and $U_0^{1/2}$.

The electromagnetic stress on an electrically conducting mantle and the topographic stress oppose the fluid motion. The electromagnetic coupling scales linearly with the differential velocity at the CMB (Roberts, 1972; Schwaiger et al., 2024). Therefore, the two couplings do not scale similarly with $U_0$. The topographic coupling's $U_0^{1/2}$ scaling implies that it always predominates in the limit of small differential velocities. The critical value of the velocity below, which the topographic stress dominates, depends on the parameters, notably on the horizontal length scale of the topography. At small wavelengths, the electromagnetic stress $\langle \tilde{F}_x \rangle_{EM} = B_{0z}^2 U_0 G$ predominates even for quite small values of velocity (with the mantle conductance $G$ and

the vertical component of the magnetic field $B_{0z}$). In contrast, the topographic stress largely dominates the electromagnetic one for large topography wavelengths ($> 1000$ km), using parameters of Figure 15 with our model of topography, and a mantle conductance of $5 \times 10^7$ S (Kuvshinov et al., 2021). The different scalings of the electromagnetic and topographic couplings mean that they are not interchangeable.

On intradecadal timescales, the Earth's core angular momentum is carried by torsional Alfvén waves (Finlay et al., 2023). This gives a special role to the equator where the waves are reflected, as discussed for electromagnetic coupling by Schaeffer and Jault (2016). Topography in the equatorial region can also play a role in the reflection of these waves. In addition, Buffett and Matsui (2019) studied MAC waves trapped at the equator by a latitudinal variation of the imposed magnetic field and in the presence of stratification. The topographic stress associated with these waves can be calculated in our framework by including a linear variation of the magnetic field.

Following Buffett (2010), we could adapt the approach presented here to investigate the time-dependent coupling associated with the diurnal differential rotation between the core and the mantle. This would allow us to seek again an explanation for the dissipation of the retrograde annual nutation by including $\beta$-plane, inertia, and Lorentz force effects in the Earth's core model. These effects might reduce the large stratification required at the CMB to explain the nutation data according to Buffett (2010). While raising additional difficulties, this is an ongoing work that has already been implemented in our code by using a double perturbation series.

A global estimate of the topographic torque in spherical geometry is necessary to compare with the geophysical observations. We suggest basing it on an integration of the local stress over the CMB, with a latitudinal dependence of the flow, magnetic field, and angle between the rotation vector and the normal to the sphere, as discussed in this work. A direct coupling of global and local models may require too much computational resources. An alternative approach consists in pre-calculating local topographic stresses for various rotation vectors, basic bulk flows, and magnetic fields. At each time step, the global model would rely on this pre-calculated parameterization of the topographic stress to drive the flow at the boundary and evolve the interior fields with time. We would thus replace the stress-free boundary condition, often used in geodynamo models, with an imposed stress boundary condition. In contrast with the local models, the global model would include the inertia term $\partial_t u$.
## Appendix A: Dispersion Relationships of MHD Waves

Considering an ansatz $\exp(i(\omega t + k_x x + k_y y + k_z z))$ for space and time dependencies, we can write the diffusionless dispersion relationship (without $\beta$-plane) of Alfvén, internal and inertial waves as, respectively (Salhi et al., 2017).

$$ \omega_M = (b_0 \cdot k) Al^{-1}, \qquad (A1) $$

$$ \omega_A = k_H k^{-1} Fr^{-1}, \qquad (A2) $$

$$ \omega_C = 2(\Omega \cdot k) Ro^{-1} k^{-1}, \qquad (A3) $$

with the wave vector $k$, and the magnitude $k_H = \sqrt{k_x^2 + k_y^2}$ of the horizontal wave number $k_H$. We can then write the dispersion relationship of inertia-gravity waves as

$$ \omega = \sqrt{\omega_C^2 + \omega_A^2}, \qquad (A4) $$

of Magneto Coriolis Waves as

$$ (\omega^2 - \omega_C \omega - \omega_M^2)(\omega^2 + \omega_C \omega - \omega_M^2) = 0, \qquad (A5) $$

and of MAC waves as

$$ (\omega^2 - \omega_M^2)(\omega^2 - \omega_M^2 - \omega_A^2) - \omega^2 \omega_C^2 = 0. \qquad (A6) $$

## Appendix B: Domain Boundaries of Figures 10 and 11

Figures 10 and 11 exhibit limits delimiting various parameter domains.

### B1. Range of $Ro$ for Inertia Gravity Waves When $Rm/Al^2 \ll 1$

Imposing a zero determinant for the system of dynamical equations (without a magnetic field and an f-plane approximation), we obtain

$$ -Rok_x^4 + 4k_x^2 Ro^{-1} \cos^2 \theta k_z^2 - k_x^6 - k_x^4 k_y^2 + Fr^{-2} (k_x^4 + k_x^2 k_y^2) = 0, \qquad (B1) $$

which gives the dispersion relationship for inertia-gravity waves. This gives the following solution for $k_z$,

$$ k_z = \pm R_o Fr^{-1} k_H \sqrt{\frac{1 - (Fr_k)^2}{(R_o k_x)^2 - 4 \cos^2 \theta}}, \qquad (B2) $$

whose real part can be searched to be non-zero. Assuming $(Frk_x)^2 < 1$, this leads to

$$ R_o k_x > 2 \cos \theta, \qquad (B3) $$

which is independent of $Fr$. Note that $\beta$-effect modifies Equation B3 as

$$ Ro > \sqrt{\chi^2 \sin^2 \theta + 4 \cos^2 \theta + \chi \sin \theta}, \qquad (B4) $$

which gives significantly different values, only close to the equator. With a tilted rotation vector (and no beta plane), $k_z$ changes sharply in value at
$$Fr^2 Ro^2 k_H^2 k_x^4 - 4Fr^2 k_x^4 k_y^2 - 4Fr^2 k_x^4 \cos^2(\theta) - Ro^2 k_H^2 k_x^2 + 4k_H^2 \cos^2(\theta) = 0. \quad (B5)$$

In Figure 15 (left), β-plane effects are negligible for the limits of Inertia-gravity waves.

## B2. Range of Al for Inertia Gravity Waves When $Rm^{-1} = 0$

To find Equation 16, that is

$$Al = 2k_H Fr^{-1} \cos \theta, \quad (B6)$$

which delimits the transition between MAC and internal waves, we consider the governing dynamical equations in the limit $Rm^{-1} = 0$, without β-plane effects. Imposing a zero determinant, and removing the smallest terms (regarding their values in our study), we obtain

$$\frac{k_x^2 \cos^4 \theta}{Al^4} k_z^6 - \frac{2k_x^4 \cos^2 \theta}{Al^2} k_z^4 + \frac{Al^2 Fr^2 k_x^6 + k_x^4 \cos^2 \theta + k_x^2 k_y^2 \cos^2 \theta}{Al^2 Fr^2} k_z^2 - \frac{k_x^6 + k_x^4 k_y^2}{Fr^2} = 0 \quad (B7)$$

whose complex solutions for $k_z$ are given by

$$k_z = \pm \sqrt{\frac{Al}{2Fr} \frac{\sqrt{AlFrk_x^2 \pm \sqrt{\Delta}}}{|\cos \theta|}}, \quad (B8)$$

with $\Delta = Al^2 Fr^2 k_x^4 - 4k_x^2 \cos^2 \theta - 4k_y^2 \cos^2 \theta$. Equation B6 is then obtained from Equation $\Delta = 0$ by using $k_x = 1$.

## B3. Rossby Waves Limits

Equations 14, 15 and 17 are, respectively, derived in the hydrodynamic and ideal MHD limits. Prescribing the $k_z$ values, as observed in our results: aligned with the rotation axis in the hydrodynamic case ($k_z = -k_y \tan \theta$), or with the magnetic field in the ideal MHD limit (here for a dipolar field, $k_z = (k_y/2) \tan \theta$). In the hydrodynamic limit, we obtain the following matrices (for which we need to obtain a zero determinant)

$$\begin{bmatrix} k_x (Rok_y \tan \theta + 2is) & -4\chi_c + 2ik_y s & k_x^2 Ro - 2ik_y s \tan \theta - ik_x R_o Fr^{-2} \\ k_x (Rok_y + 2ic) & -Rok_x^2 + 2\chi_s + 2ik_y c & 4\chi_c - 2ik_y s & 0 \\ ik_x & ik_y & -ik_y \tan \theta & 0 \\ 0 & 0 & -1 & ik_x \end{bmatrix} \quad (B9)$$

with $s = \sin \theta$ and $c = \cos \theta$. In Equation 13, given by the determinant of matrix B9, we have.

$$A^2 = (k_H^2 - k_y^2)^2 + [(k_H k_y)^2 - k_y^4] (\cos \theta)^{-2}, \quad (B10)$$

$$B = k_H^2 + 2k_y^2 + k_y^2 (\cos \theta)^{-2}, \quad (B11)$$

$$C = 16 \cos^2 \theta \quad (B12)$$

Equation 15 is obtained for $Fr \ll 1$. In the ideal MHD limit, we obtain the matrix
$$ \mathcal{M} = \begin{bmatrix} \mathcal{A} & 0 \\ 0 & \mathcal{B} \end{bmatrix} \qquad (B13) $$

where $\mathcal{A}$ is given by

$$ \begin{bmatrix} k_x(-Rok_y \tan \theta + 4is)/2 & -4\chi - ik_y s & Rok_x^2 + ik_y s \tan \theta - ik_x R_o F r^{-2} \\ k_x(Rok_y + 2i) & -Rok_x^2 + 2\chi s + 2ik_y & 4\chi - 2ik_y s & 0 \\ ik_x & ik_y & (i/2)k_y \tan \theta & 0 \\ 0 & 0 & -1 & ik_x \end{bmatrix}, \qquad (B14) $$

noting $s = \sin \theta$, $c = \cos \theta$, and with

$$ \mathcal{B} = \begin{bmatrix} ik_x & ik_y & (i/2)k_y \tan \theta \\ 0 & ik_x & 0 \\ 0 & 0 & ik_x \end{bmatrix}. \qquad (B15) $$

Equation 17 can be obtained by canceling the determinant of $\mathcal{M}$ in the limit $Fr, Ro \ll 1$. Note that the matrix $\mathcal{M}$ given by Equation B13 is given for non-traditional $\beta$-plane and is naturally modified for other approximations.

## Appendix C: Regime $Le < Fr$

Equation 20 is found by calculating the dispersion relationship of MAC-Rossby waves in the parameter regime of Figure 11. We obtain the polynomial equation

$$ \begin{aligned} & -\frac{k_z^6}{Al^4} + k_z^4 \left( \frac{2}{Al^2} - \frac{2\chi \sin \theta}{Al^2 R_o} - \frac{1}{Al^4} \right) \\ & + k_z^2 \left( -1 + \frac{2\chi \sin \theta}{R_o} + \frac{4\cos^2 \theta}{R_o^2} + \frac{2}{Al^2} - \frac{2\chi \sin \theta}{Al^2 R_o} - \frac{1}{Al^2 Fr^2} \right) \\ & - 1 + \frac{2\chi \sin \theta}{R_o} + \frac{1}{Fr^2} - \frac{2\chi \sin \theta}{Fr^2 R_o} = 0. \end{aligned} \qquad (C1) $$

The transition between MAC and MAC-Rossby occurs when the coefficient of $k_z^2$ is canceled. We obtain

$$ Fr^2 = \frac{Ro^2}{Al^2(2Ro\chi \sin \theta + 4\cos^2\theta - Ro^2) + 2Ro(R_o - \chi \sin \theta)} \qquad (C2) $$

which simplifies into $Fr = Le/(2 \cos \theta)$ for $Ro \ll 1$. For Froude numbers smaller than the transition (Equation C2), the imaginary part of $k_z$ vanishes.

Figure C1 shows the equivalent of Figure 8 for $Le < Fr$. In this opposite limit, we can see that internal Rossby waves (non-magnetic) are even present at large values of the interaction parameter $Rm/Al^2$. This explains the stress increase for $Fr > Le$ in Figure 11.
<figure><img src="image_19.png" /></figure>

**Figure C1.** Normalized stress (colors) as a function of the interaction parameter $Rm/Al^2$ and $Ro$, at $\theta = \pi/4$ for the 2D topography of Figure 3, a non-traditional $\beta$-plane approximation (Equation 7), a dipolar magnetic field $b_0 = (\sin \theta)/2l_y - \cos \theta l_z$. Blue dotted lines are Equation 15, and the red one is Equation 12. The black dashed-dotted line is $\Lambda = 1$ and the purple vertical line is $Al = Fr$. Parameters: $\chi = 0.005737$ ($\tilde{k}^{-1} = 20$ km for Earth core radius), $Le = 10^{-4}$, $Fr = 3 \times 10^{-2}$.

## Appendix D: The Waves in Figure 15

To better characterize the waves that carry the stress shown in Figure 15, we adopt an approach similar to that followed in §4.1.2. This figure shows the stress as a function of the velocity $U_0$ for Earth-like values of all other parameters. We compare the stress calculated with $\tilde{k}^{-1} = 5$ km and either $\eta = 0.8 \text{ m}^2 \cdot \text{s}^{-1}$ (geophysical case, top row of Figure D1) or vanishing magnetic diffusivity (highly conducting fluid, bottom row). Switching from Earth-like to high electrical conductivity causes the stress to drop when $U_0 < 10^{-4} \text{ m} \cdot \text{s}^{-1}$. This corresponds to the domain of existence of a branch of Rossby-MAC waves that participate in the boundary condition although associated with negligible stress. In the high conductivity limit and for $U_0 < 10^{-4} \text{ m} \cdot \text{s}^{-1}$, the MAC wave still carries most of the stress but has low amplitude in the presence of Rossby-MAC waves. This explains the drop in the stress observed only in the high conductivity limit. At $U_0 = 5 \times 10^{-3} \text{ m} \cdot \text{s}^{-1}$ and for both $\eta = 0.8 \text{ m}^2 \cdot \text{s}^{-1}$ and vanishing magnetic diffusivity, the change in the slope of $\langle \tilde{F}_x \rangle$ as a function of $U_0$ occurs when inertia-gravity waves replace MAC waves as the main carrier of the stress.
<figure><img src="image_20.png" /><figcaption>**Figure D1.** Normalized stress (blue solid line at the top of each figure) as a function of $U_0$ and associated absolute value of the real (left) or imaginary (right) part of non dimensional vertical wavenumbers $k_z$ for $\bar{k}^{-1} = 5$ km. Vertical wave numbers are color-coded (red to yellow) to show the relative contribution of each $k_z$ to the total stress. Top figures (a, b) are calculated for the parameters of Figure 15, left. The bottom ones (c, d) are in the ideal MHD limit $Rm = 10^{20}$. Other parameters are the same as Figure 15. Other lines show diffusionless MHD wave dispersion relationships (see Appendix A or Salhi et al., 2017).</figcaption></figure>

## Acknowledgments

We thank Thomas Frasson for his help on the CMB height model. We sincerely thank an anonymous referee and Mathieu Dumberry, for their careful reading of the manuscript and the many comments, which greatly helped us to improve the quality of the paper. This work is supported by the ERC under the European Union's Horizon 2020 research and innovation programme via the THEIA project (grant agreement no. 847433). ISTerre is part of Labex OSUG@2020 (ANR10 LABX56). Computations have been performed on the local cluster of ISTerre. DC acknowledges support from GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), and HPC resources (Jean Zay V100) of IDRIS under allocation AD010413621 attributed by GENCI (Grand Equipement National de Calcul Intensif).

## Data Availability Statement

The code used to produce the result of this article is open-source and available at Monville (2024) under the CeCILL Licence. This code is built using the Sympy (Meurer et al., 2017) symbolic computation library and the mpmath arbitrary precision arithmetic library (The mpmath development team, 2023). Figures are plotted using Matplotlib library (Hunter, 2007). Colour maps from Crameri (2018) prevent visual distortion of the data and exclusion of readers with colour-vision deficiencies.

## References

* Acheson, D. J. (1975). On hydromagnetic oscillations within the Earth and core-mantle coupling. *Geophysical Journal International*, **43**(1), 253-268. https://doi.org/10.1111/j.1365-246X.1975.tb00634.x

* Arnould, M., Coltice, N., Flament, N., & Mallard, C. (2020). Plate tectonics and mantle controls on plume dynamics. *Earth and Planetary Science Letters*, **116439**. https://doi.org/10.1016/j.epsl.2020.116439

* Athanassiadou, M. (2003). Wave and form drag: Their relation in the linear gravity wave regime. *Tellus A: Dynamic Meteorology and Oceanography*, **55**(2), 173-180. https://doi.org/10.3402/tellusa.v55i2.12090

* Bell, T. H. (1975a). Lee waves in stratified flows with simple harmonic time dependence. *Journal of Fluid Mechanics*, **67**(4), 705-722. https://doi.org/10.1017/S0022112075000560

* Bell, T. H. (1975b). Topographically generated internal waves in the open ocean. *Journal of Geophysical Research*, **80**(3), 320-327. https://doi.org/10.1029/JC080i003p00320
Berkowitz, S. J. (1984). On computing the determinant in small parallel time using a small number of processors. *Information Processing Letters*, **18**(3), 147–150. https://doi.org/10.1016/0020-0190(84)90018-8

Bouffard, M., Landeau, M., & Goument, A. (2020). Convective erosion of a primordial stratification atop Earth's core. *Geophysical Research Letters*, e2020GL087109. https://doi.org/10.1029/2020GL087109

Braginsky, S. I. (1998). Magnetic Rossby waves in the stratified ocean of the core, and topographic core-mantle coupling. *Earth Planets and Space*, **50**(8), 641–649. https://doi.org/10.1186/BF03352159

Braginsky, S. I. (1999). Dynamics of the stably stratified ocean at the top of the core. *Physics of the Earth and Planetary Interiors*, **111**(1–2), 21–34. https://doi.org/10.1016/s0031-9201(98)00143-5

Buffett, B. A. (1996). Gravitational oscillations in the length of day. *Geophysical Research Letters*, **23**(17), 2279–2282. https://doi.org/10.1029/96GL02083

Buffett, B. A. (2010). Chemical stratification at the top of Earth's core: Constraints from observations of nutations. *Earth and Planetary Science Letters*, **296**(3–4), 367–372. https://doi.org/10.1016/j.epsl.2010.05.020

Buffett, B. A., & Christensen, U. R. (2007). Magnetic and viscous coupling at the core—Mantle boundary: Inferences from observations of the Earth's nutations. *Geophysical Journal International*, **145**, 152. https://doi.org/10.1111/j.1365-246X.2007.03543.x

Buffett, B. A., Knezek, N., & Holme, R. (2016). Evidence for MAC waves at the top of Earth's core and implications for variations in length of day. *Geophysical Journal International*, **1789–1800**. https://doi.org/10.1093/gji/ggv552

Buffett, B. A., & Matsui, H. (2019). Equatorially trapped waves in Earth's core. *Geophysical Journal International*, **1210–1225**. https://doi.org/10.1093/gji/ggz233

Coltice, N., Husson, L., Faccenna, C., & Arnould, M. (2019). What drives tectonic plates? *Science Advances*, eaax4295. https://doi.org/10.1126/sciadv.aax4295

Crameri, F. (2018). *Scientific colour maps*. Zenodo.

Davies, C. J., Stegman, D. R., & Dumberry, M. (2014). The strength of gravitational core-mantle coupling. *Geophysical Research Letters*, **3786–3792**. https://doi.org/10.1002/2014GL059836

Dehant, V., Campuzano, S. A., De Santis, A., & van Westrenen, W. (2022). Structure, materials and processes in the Earth's core and mantle. *Surveys in Geophysics*, **263–302**. https://doi.org/10.1007/s10712-021-09684-y

Deleplace, B., & Cardin, P. (2006). Viscomagnetic torque at the core mantle boundary. *Geophysical Journal International*, **557–566**. https://doi.org/10.1111/j.1365-246X.2006.03180.x

Dellar, P. J. (2011). Variations on a beta-plane: Derivation of non-traditional beta-plane equations from Hamilton's principle on a sphere. *Journal of Fluid Mechanics*, **674**, 174–195. https://doi.org/10.1017/S0022112010006464

Finlay, C. C. (2008). Course 8 Waves in the presence of magnetic fields, rotation and convection. In *Les houches* (pp. 403–450). Elsevier. https://doi.org/10.1016/s0924-8099(08)80012-1

Finlay, C. C., Gillet, N., Aubert, J., Livermore, P. W., & Jault, D. (2023). Gyres, jets and waves in the Earth's core. *Nature Reviews Earth and Environment*, **4**(6), 1–16. https://doi.org/10.1038/s43017-023-00425-w

Forte, A. M., & Rowley, D. B. (2022). Earth's isostatic and dynamic topography—A critical perspective. *Geochemistry, Geophysics, Geosystems*, **23**(9), e2021GC009740. https://doi.org/10.1029/2021GC009740

Gastine, T., Aubert, J., & Fournier, A. (2020). Dynamo-based limit to the extent of a stable layer atop Earth's core. *Geophysical Journal International*, **1433–1448**. https://doi.org/10.1093/gji/gga250

Gerick, F., Jault, D., Noir, J., & Vidal, J. (2020). Pressure torque of torsional Alfvén modes acting on an ellipsoidal mantle. *Geophysical Journal International*, **338–351**. https://doi.org/10.1093/gji/gga166

Gillet, N., Jault, D., & Finlay, C. C. (2015). Planetary gyre, time-dependent eddies, torsional waves, and equatorial jets at the Earth's core surface. *Journal of Geophysical Research: Solid Earth*, **3991–4013**. https://doi.org/10.1002/2014JB011786

Glane, S., & Buffett, B. A. (2018). Enhanced core-mantle coupling due to stratification at the top of the core. *Frontiers in Earth Science*. https://doi.org/10.3389/feart.2018.00171/full

Golub, G. H., & Reinsch, C. (1971). Singular value decomposition and least squares solutions. In F. L. Bauer, A. S. Householder, F. W. J. Olver, H. Rutishauser, K. Samelson, & E. Stiefel (Eds.), *Handbook for automatic computation: Volume II: Linear algebra* (pp. 134–151). Springer Berlin Heidelberg. https://doi.org/10.1007/978-3-642-86940-2_10

Grimshaw, R. H. J. (1975). A note on the β-plane approximation. *Tellus*, **27**(4), 351–357. https://doi.org/10.1111/j.2153-3490.1975.tb01685.x

Gross, R. S., Fukumori, I., Menemenlis, D., & Gegout, P. (2004). Atmospheric and oceanic excitation of length-of-day variations during 1980–2000. *Journal of Geophysical Research*. https://doi.org/10.1029/2003JB002432

Helffrich, G., & Kaneshima, S. (2010). Outer-core compositional stratification from observed core wave speed profiles. *Nature*, **807–810**. https://doi.org/10.1038/nature09636

Heyn, B. H., Conrad, C. P., & Trönnes, R. G. (2020). Core-mantle boundary topography and its relation to the viscosity structure of the lowermost mantle. *Earth and Planetary Science Letters*, **543**, 116358. https://doi.org/10.1016/j.epsl.2020.116358

Hide, R. (1969). Interaction between the Earth's liquid core and solid mantle. *Nature*, **222**(5198), 1055–1056. https://doi.org/10.1038/2221055a0

Holme, R. (1998). Electromagnetic core—Mantle coupling—I. Explaining decadal changes in the length of day. *Geophysical Journal International*, **132**(1), 167–180. https://doi.org/10.1046/j.1365-246x.1998.00424.x

Holme, R., & De Viron, O. (2013). Characterization and implications of intradecadal variations in length of day. *Nature*, **202–204**, https://doi.org/10.1038/nature12282

Hunter, J. D. (2007). Matplotlib: A 2D graphics environment. *Computing in Science and Engineering*, **90–95**. https://doi.org/10.1109/MCSC.2007.55

Ichikawa, H., Tsuchiya, T., & Tange, Y. (2014). The P-V-T equation of state and thermodynamic properties of liquid iron. *Journal of Geophysical Research: Solid Earth*, **119**(1), 240–252. https://doi.org/10.1002/2013JB010732

Irving, J. C. E., Cottaar, S., & Lekić, V. (2018). Seismically determined elastic parameters for Earth's outer core. *Science Advances*, eaar2538. https://doi.org/10.1126/sciadv.aar2538

Jault, D. (2020). Tangential stress at the core–mantle interface. *Geophysical Journal International*, **221**(2), 951–967. https://doi.org/10.1093/gji/gga048

Jault, D., Gire, C., & Le Mouël, J. L. (1988). Westward drift, core motions and exchanges of angular momentum between core and mantle. *Nature*, **353–356**. https://doi.org/10.1038/333353a0

Jault, D., & Mouël, J.-L. L. (1989). The topographic torque associated with a tangentially geostrophic motion at the core surface and inferences on the flow inside the core. *Geophysical and Astrophysical Fluid Dynamics*, 273–295. https://doi.org/10.1080/03091928908218533

Kaneshima, S. (2018). Array analyses of SmKS waves and the stratification of Earth's outermost core. *Physics of the Earth and Planetary Interiors*, **276**, 234–246. https://doi.org/10.1016/j.pepi.2017.03.006
Kerner, I. O. (1966). Ein gesamtschrittverfahren zur berechnung der nullstellen von polynomen. *Numerische Mathematik*, 8(3), 290–294. https://doi.org/10.1007/BF02162564

Kiani Shahvandi, M., Adhikari, S., Dumberry, M., Mishra, S., & Soja, B. (2024). The increasingly dominant role of climate change on length of day variations. In *Proceedings of the national academy of sciences*. https://doi.org/10.1073/pnas.2406930121.e2406930121.

Koelemeijer, P. (2021). Toward consistent seismological models of the core–mantle boundary landscape. In *Mantle convection and surface expressions* (pp. 229–255). American Geophysical Union (AGU). https://doi.org/10.1002/9781119528609.ch9

Koroni, M., Borgeaud, A., Fichtner, A., & Deschamps, F. (2023). An analysis of core–mantle boundary related seismic waves using full-waveform modelling and adjoint methods. *Geophysical Journal International*, 232(2), 1259–1275. https://doi.org/10.1093/gji/ggac389

Kuang, W., & Chao, B. F. (2001). Topographic core-mantle coupling in geodynamo modeling. *Geophysical Research Letters*, 1871–1874. https://doi.org/10.1029/2000GL012237

Kuvshinov, A., Grayver, A., Toffner-Clausen, L., & Olsen, N. (2021). Probing 3-D electrical conductivity of the mantle using 6 years of Swarm, CryoSat-2 and observatory magnetic data and exploiting matrix Q-responses approach. *Earth Planets and Space*, 67. https://doi.org/10.1186/s40623-020-01341-9

Landeau, M., Olson, P., Deguen, R., & Hirsh, B. H. (2016). Core merging and stratification following giant impact. *Nature Geoscience*, 786–789. https://doi.org/10.1038/ngeo2808

Lau, H. C., Mitrovica, J. X., Davis, J. L., Tromp, J., Yang, H.-Y., & Al-Attar, D. (2017). Tidal tomography constrains Earth's deep-mantle buoyancy. *Nature*, 321–326. https://doi.org/10.1038/nature24452

Legg, S. (2021). Mixing by oceanic lee waves. *Annual Review of Fluid Mechanics*, 53(1), 173–201. https://doi.org/10.1146/annurev-fluid-051220-043904

Meurer, A., Smith, C. P., Paprocki, M., Čertík, O., Kirpichev, S. B., Rocklin, M., et al. (2017). SymPy: Symbolic computing in Python. *PeerJ Computer Science*, e103, e103. https://peerj.com/articles/cs-103.pdf

Moffatt, H. K. (1977). Topographic coupling at the core-mantle interface. *Geophysical and Astrophysical Fluid Dynamics*, 279–288. https://doi.org/10.1080/03091927708242332

Monville, R. (2024). ToCCo [Software]. Zenodo. https://doi.org/10.5281/zenodo.14002373

Muir, J. B., Tanaka, S., & Tkalčić, H. (2022). Long-wavelength topography and multiscale velocity heterogeneity at the core-mantle boundary. *Geophysical Research Letters*, 49(17), e2022GL099943. https://doi.org/10.1029/2022GL099943

Narteau, C., Le Mouël, J. L., Poirier, J. P., Sepúlveda, E., & Shnirman, M. (2001). On a small-scale roughness of the core–mantle boundary. *Earth and Planetary Science Letters*, 49–60. https://doi.org/10.1016/S0012-821X(01)00041-0

Nikurashin, M., & Ferrari, R. (2010). Radiation and dissipation of internal waves generated by geostrophic motions impinging on small-scale topography: Application to the southern ocean. *Journal of Physical Oceanography*, 2025–2042. https://doi.org/10.1175/2010JPO4315.1

Olson, P., Landeau, M., & Reynolds, E. (2018). Outer core stratification from the high latitude structure of the geomagnetic field. *Frontiers in Earth Science*, 140. https://doi.org/10.3389/feart.2018.00140/full

Pais, M. A., & Jault, D. (2008). Quasi-geostrophic flows responsible for the secular variation of the Earth's magnetic field. *Geophysical Journal International*, 421–443. https://doi.org/10.1111/j.1365-246X.2008.03741.x

Pedlosky, J. (2013). *Geophysical fluid dynamics*. Springer Science and Business Media.

Puica, M., Dehant, V., Folgueira, M., Van Hoolst, T., & Rekier, J. (2023). Analytical computation of total topographic torque at the core–mantle boundary and its impact on tidally driven length-of-day variations. *Geophysical Journal International*, 585–596. https://doi.org/10.1093/gji/ggad077

Rexer, M., & Hirt, C. (2015). Ultra-high-Degree surface spherical harmonic analysis using the gauss–legendre and the driscoll/healy quadrature theorem and application to planetary topography models of Earth, Mars and Moon. *Surveys in Geophysics*, 803–830. https://doi.org/10.1007/s10712-015-9345-z

Roberts, P. H. (1972). Electromagnetic core-mantle coupling. *Journal of Geomagnetism and Geoelectricity*, 24(2), 231–259. https://doi.org/10.5636/jgg.24.231

Roberts, P. H., & Aurnou, J. M. (2012). On the theory of core-mantle coupling. *Geophysical and Astrophysical Fluid Dynamics*, 157–230. https://doi.org/10.1080/03091929.2011.589028

Rosat, S., & Gillet, N. (2023). Intradecadal variations in length of day: Coherence with models of the Earth's core dynamics. *Physics of the Earth and Planetary Interiors*, 107053. https://doi.org/10.1029/j.pepi.2023.107053

Rossby, C.-G. (1939). Relation between variations in the intensity of the zonal circulation of the atmosphere and the displacements of the semi-permanent centers of action. *Journal of Marine Research*, 2(1), 38–55. https://elischolar.library.yale.edu/journal_of_marine_research/544/

Salhi, A., Baklouti, F. S., Godeferd, F., Lehner, T., & Cambon, C. (2017). Energy partition, scale by scale, in magnetic Archimedes Coriolis weak wave turbulence. *Physical Review E*, 023112. https://doi.org/10.1103/PhysRevE.95.023112

Schaeffer, N., & Jault, D. (2016). Electrical conductivity of the lowermost mantle explains absorption of core torsional waves at the equator. *Geophysical Research Letters*, 4922–4928. https://doi.org/10.1002/2016GL068301

Schwaiger, T., Gillet, N., Jault, D., Istas, M., & Mandea, M. (2024). Wave-like motions and torques in Earth's core as inferred from geomagnetic data: A synthetic study. *Physics of the Earth and Planetary Interiors*, 346, 107104. https://doi.org/10.1016/j.pepi.2023.107104

Sommerfeld, A. (1912). *Die Greensche Funktion der Schwingungslgleichung* (pp. 309–352). Jahresbericht der Deutschen Mathematiker-Vereinigung.

The mpmath development team. (2023). Mpmath: A Python library for arbitrary-precision floating-point arithmetic (version 1.3.0) [Manual].

Thompson, A. R., Moran, J. M., & Swenson, G. W. (2017). Very-long-baseline interferometry. In A. R. Thompson, J. M. Moran, & G. W. JrSwenson (Eds.), *Interferometry and synthesis in radio astronomy* (pp. 391–483). Springer International Publishing. https://doi.org/10.1007/978-3-319-44431-4_9

Triana, S. A., Trinh, A., Rekier, J., Zhu, P., & Dehant, V. (2021). The viscous and Ohmic damping of the Earth's free core nutation. *Journal of Geophysical Research: Solid Earth*, e2020JB021042. https://doi.org/10.1029/2020JB021042

Van Tent, R., Deuss, A., Kaneshima, S., & Thomas, C. (2020). The signal of outermost-core stratification in body-wave and normal-mode data. *Geophysical Journal International*, 223(2), 1338–1354. https://doi.org/10.1093/gji/ggaa368