Reference (5)

Liquid-Phase

When the melt fraction of a magma layer exceeds a critical value $F_m \gtrsim F_{m,c}$, its response to tidal forcing becomes fluid-like rather than solid-like. In this regime, tidal deformation is governed by the linearized Laplace tidal equations (LTEs), which describe the dynamics of a thin, global fluid shell of uniform thickness $H$ on a planet of radius $R_p$.

Neglecting nonlinear advection and mean flows, momentum and mass conservation for the magma ocean read

\[\partial_t \vec{u} + \sigma_R \vec{u} + \vec{f} \times \vec{u} + g \nabla \zeta = g \nabla \zeta_{\mathrm{eq}}, \qquad \partial_t \zeta + \nabla \cdot (H \vec{u}) = 0,\]

where $\vec{u}$ is the horizontal velocity, $\zeta$ the tidally induced surface displacement, $\zeta_{\mathrm{eq}} = U_T/g$ the equilibrium tide associated with the tidal potential $U_T$, $g$ the surface gravity, and $\sigma_R$ a Rayleigh drag frequency parameterizing dissipation in the magma ocean.

We focus on the strongly damped, highly viscous limit ($\sigma_R \gg 2\Omega$), appropriate for magma oceans, in which dissipative processes dominate over rotational effects. In this high-Ekman-number regime, the Coriolis term is negligible and the tidal response becomes overdamped. Dissipation is encapsulated by $\sigma_R$, which should be interpreted as an effective damping rate accounting for viscous resistance, boundary-layer friction, form drag, and porous-flow dissipation near the rheological transition.

Transforming the governing equations to the frequency domain and expanding tidal quantities in spherical harmonics yields a closed-form expression for the complex, frequency-dependent tidal Love number $k_\ell(\sigma)$.

\[k_{\ell}(\sigma) = -\frac{\varrho_\ell \bar{\sigma}_\ell^2}{\sigma\tilde{\sigma} - \bar{\sigma}_\ell^2},\]

where $\varrho_\ell = 3/(2n+1)(\rho_f/\rho_b)$ is the degree-$\ell$ density ratio between the fluid and the bulk part of the planet. The characteristic frequency is given by

\[\bar{\sigma}_\ell = \sqrt{\frac{\mu_\ell g H_\mathrm{magma}}{R^2}},\]

where $H_\mathrm{magma}$ the ocean depth, $R$ the radius of the planet, and with

\[\mu_\ell = \ell(\ell+1),\]

the eigenvalue of the nth degree spherical harmonic, and

\[\tilde{\sigma} = \sigma - i \sigma_R.\]

Radially resolved dissipation

Given that there is no consensus on where in a fluid layer most tidal energy is dissipated, we provide several parameterized radial dissipation profiles. These profiles represent different physical assumptions about how turbulent mixing, wave breaking, or viscous damping distribute tidal energy with depth. The available profiles should be interpreted as idealized parameterizations of where energy is dissipated within the fluid layer.

Uniform dissipation

The uniform profile assumes dissipation is constant throughout the layer.

This corresponds to the simplest assumption that turbulence or small-scale mixing distributes energy evenly with depth. It is useful as a baseline case when no information about the dissipation structure is available.

Exponential dissipation

The exponential profile concentrates dissipation near the lower boundary and decreases with height:

\[D(z) \propto e^{-z/H_R}\]

This can approximate scenarios where dissipation is strongest near the interface with the underlying solid layer, for example when tidal flows interact with boundary roughness or generate shear-driven turbulence.

Linear dissipation

The linear profile decreases linearly with height and reaches zero at a characteristic scale (H_R).

This provides a simple way to model dissipation that is concentrated toward the bottom of the fluid layer but not as strongly localized as an exponential profile.

Quadratic dissipation

The quadratic profile falls off more steeply with height:

\[D(z) \propto (1 - z/H_R)^2\]

This concentrates dissipation even more strongly toward the base of the layer. It can approximate cases where turbulence is primarily generated near the boundary and decays rapidly away from it.

Dynamic dissipation

The dynamic profile introduces a depth-dependent mixing length

\[\ell_{\text{mix}} = \min(z, H_R)\]

so that the dissipation scale adjusts with distance from the boundary.

This mimics mixing-length arguments commonly used in geophysical and astrophysical fluid dynamics, where turbulence intensity depends on the available eddy size.

Dynamic-interpolated dissipation (dynamic_interp)

The dynamic_interp profile (the default sigma_R_prf) replaces the single depth-dependent shape function above with a composite of three physically distinct contributions, each assigned its own share of the total dissipated power so that the three shares sum exactly back to the total:

\[P(z) = P_{\mathrm{sbd}}(z) + P_{\mathrm{drag}}(z) + P_{\mathrm{fric}}(z), \qquad \int P\,dV = E_{\mathrm{sbd}} + E_{\mathrm{drag}} + E_{\mathrm{fric}} = E_{\mathrm{total}}.\]

• Shear/bulk/Darcy component ($E_{\mathrm{sbd}}$): the shear-, bulk-, and Darcy-heating already computed at the interface (e.g. from a neighbouring mush segment) is carried into the fluid layer with a Gaussian-shoulder radial decay from the top interface, $\propto \exp[-(z/H_{\mathrm{decay}})^{1.5}]$. The decay length starts at the configured scale height $H_R$, but is reduced (found by bisection) if needed so that this component never exceeds 25% of the frequency-independent baseline energy $E_{\infty}$ (the energy dissipated in the $\sigma_R \to \sigma_{R,\infty}$ bulk-fluid limit).
• Drag component ($E_{\mathrm{drag}} = \max(E_{\infty} - E_{\mathrm{sbd}}, 0)$): the remaining bulk Rayleigh-drag dissipation, distributed with a smooth ($\tanh$) sigmoid that switches on around the depth $z_{\mathrm{visc}}$ where the viscosity profile first drops to the pure-liquid viscosity visc_l, i.e. at the base of the mush-to-liquid transition.
• Friction component ($E_{\mathrm{fric}} = \max(E_{\mathrm{total}} - E_{\infty}, 0)$): the excess dissipation associated with the interfacial (frequency-dependent) drag term, represented as a Gaussian pulse centred below the interface (at $z = 2H_R$, width $H_R/4$) and forced to vanish at $z=0$.

This construction is designed to interpolate smoothly between the mush-transition heating handed off from the solid side and the bulk fluid response, without requiring the user to hand-pick a single idealized shape.

Choosing a dissipation profile

Since the physical location of tidal energy dissipation is uncertain in many systems, these profiles should be viewed as parameterized hypotheses. Comparing results across multiple profiles is often more informative than adopting a single choice.

In practice:

Function Documentation

run_fluid1d(omega, rho, radius, visc, ρ_ratio, P_b, R; n=2, sigma_R=1e-3, sigma_inf=1e-7, sigma_R_prf="uniform", H_R=1e3, efficiency=0.3, visc_l=1e2)

compute_fluid_lovenumbers(omega, R, H_magma, g, ρ_ratio, n, σ_R)