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

HAL
open science

# A magnetically driven equatorial jet in Europa's ocean

Christophe Gissinger, Ludovic Petitdemange

► To cite this version:

Christophe Gissinger, Ludovic Petitdemange. A magnetically driven equatorial jet in Europa's ocean. Nature Astronomy, 2019, 3 (5), pp.401-407. 10.1038/s41550-019-0713-3. hal-02381718

HAL Id: hal-02381718

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

Submitted on 8 Oct 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A magnetically driven equatorial jet in Europa's ocean

Christophe Gissinger<sup>1,3*</sup> and Ludovic Petitdemange<sup>2,3</sup>

During recent decades, data from space missions have provided strong evidence of deep liquid oceans underneath a thin outer icy crust on several moons of Jupiter<sup>1,2</sup>, particularly Europa<sup>3,4</sup>. But these observations have also raised many unanswered questions regarding the oceanic motions generated under the ice, or the mechanisms leading to the geological features observed on Europa<sup>5,6</sup>. By means of direct numerical simulations of Europa's interior, we show here that Jupiter's magnetic field generates a retrograde oceanic jet at the equator, which may influence the global dynamics of Europa's ocean and contribute to the formation of some of its surface features by applying a unidirectional torque on Europa's ice shell.

Whereas both radiogenic and tidal heating<sup>8,9</sup> produce the energy dissipation necessary for the melting of ice<sup>10,11</sup>, motions in the ocean underneath the Jovian moons are believed to be generated through vigorous thermal convection<sup>12</sup>, hydrothermal plumes<sup>13,14</sup> or double-diffusion convection<sup>15</sup>. Such flows certainly play a dominant role, but may fail at explaining some of the observations if considered alone<sup>7</sup>, strongly suggesting the presence of an additional physical mechanism in these oceans. Because the magnetic dipole axis is tilted by about 10° with the rotation axis of the gaseous giant, Jupiter's moons also experience a time-varying magnetic field with a rotation rate ω, inducing electrical currents in the oceanic salty water<sup>16</sup>.

Here, we argue that as long as the phase lag between the induced field and the Jovian one is non-zero, these induced currents naturally combine with the magnetic field to generate a Lorentz force, leading to a weak magnetohydrodynamic (MHD) process that might play a significant role in the global dynamics of the ocean. We therefore model Europa's interior as a spherical shell (mean radius $R=(R_i + R_E)/2$, thickness $h=R_E-R_i$) of salty water (electrical conductivity $\sigma$, magnetic permeability $\mu_0$, density $\rho$ and kinematic viscosity $\nu$) confined between an inner mantle of silicate rocks (radius $R_i$) and an outer layer (radius $R_E$) of ice crust (see Methods for a definition of the control parameters). We specifically model Europa here, but our results should apply equally to subsurface oceans found in other Jovian moons. To focus on the MHD process, we first present simulations in which thermal buoyancy is neglected. In the last part of this paper and in the Supplementary Information, we show that buoyancy only weakly modifies the magnetically driven jet, but remains crucial to get a full picture of the ocean dynamics.

Figure 1 shows that the rotation of Jupiter's magnetic field induces a planetary-scale recirculation, and generates strong upward and downward turbulent motions at the equator (Fig. 1c). But the most striking feature is the generation of a powerful oceanic jet propagating westward (Fig. 1a,b), and localized in the moon's equatorial region. The Jupiter-Europa system can therefore be regarded as a gigantic induction electromagnetic pump, in which the salty water of the subsurface ocean is electromagnetically pumped at a mean

velocity $\bar{U}$ by the variations of the Jovian magnetic field travelling in the horizontal plane at speed $c=\omega R \sim 230 \text{ m s}^{-1}$. This can be easily understood from the induction equation:

$$ \frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{U} \times \mathbf{B}) + \frac{1}{\mu_0 \sigma} \nabla^2 \mathbf{B} \quad (1) $$

which governs the evolution of the magnetic field $B(t)$ inside Europa's ocean. Using the spherical coordinate system ($e_\phi, e_\rho, e_\phi$), if one assumes that the field can be written ($B_r^\prime e_r + B_\phi^\prime e_\phi$) $e^{i(\phi-\omega t)}$ and that the induced currents $J=\sigma(U(r)-c)B_\phi e_\phi$ are also travelling waves, it follows that the mean Lorentz force acting on the ocean can be written (see Methods):

$$ \bar{\mathbf{F}} = \frac{\sigma B_0^2 R^2 (c - \bar{U})}{h^2 (8 + 2\mu_0^2 \sigma^2 R^2 (c - \bar{U})^2)} e_\phi \quad (2) $$

This expression of the Lorentz force is well known in the context of electromagnetic pumps<sup>17,18</sup>, and describes how a small driving of the ocean is produced as long as $\sigma$ is finite. An important feature is the phase lag $\phi_i$ between the Jovian field $B_0$ and the induced one, which controls the torque applied on the ocean. Our simulations span a large range of values of $\phi_i$, but cases corresponding to Jovian moons exhibit phase lags between 1 and 3 degrees (Supplementary Fig. 3). Accordingly, the Lorentz force is very small, and these magnetically driven jets are expected to be very weak compared with the velocity of the Jovian field. In other words, Jovian moons are inefficient induction pumps producing induced fields almost identical to the ones predicted by non-MHD studies.

Note that the dimensioned values indicated in Fig. 1 are specific to the set of parameters used here, still very far from the one relevant to real moons. Thus, to estimate the value of the Lorentz force $F$ acting on Europa's ocean, we next compare our direct numerical simulations (DNS) to observations. Following previous studies, we use Galileo measurements of Europa's induced magnetic field to constrain the value of the ocean's electrical conductivity. To wit, we extensively explored our parameter space and compared our results to the spacecraft mission's flyby E14. As shown in Fig. 2, we found that the best fits are obtained for $\sigma$ in the range $0.3-3 \text{ S m}^{-1}$ if the ocean depth is $h=147 \text{ km}$. This value, corresponding to a magnetic Reynolds number $Rm=\mu_0\omega\phi$ between 8 and 80, is in good agreement with previous predictions<sup>19</sup> and implies a salinity comparable to terrestrial oceans<sup>20</sup>.

Quantitative predictions clearly require to identify which term balances the time-averaged Lorentz force $F$ at large scale. On azimuthal averaging, thermal buoyancy does not contribute much, and the central question is to determine whether the viscous force or the nonlinear term balances $\bar{F}$. As shown in Methods, a purely viscous balance (ignoring nonlinear terms and global rotation) leads to velocity of the oceanic jet such that:

<footnote n="1">¹Laboratoire de Physique de l'Ecole Normale Superieure, ENS, Université PSL, CNRS, Paris, France. ²LERMA, CNRS, Paris, France. ³These authors contributed equally: Christophe Gissinger, Ludovic Petitdemange. *e-mail: christophe.gissinger@lps.ens.fr</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** | Velocity field and ohmic currents. **a–d**, Snapshot of the azimuthal $U_\phi$ and polar $U_\theta$ components of the velocity field (**a–c**) and of the ohmic currents $\sqrt{|J|^2}$ (**d**). The fields are shown in the $\phi,\theta$ plane (elliptical projection) for one of our global models of Europa, $\text{Ek}=10^{-7}$, $\Lambda=10^{-1}$, $\text{Pm}=10^{-4}$, $h=100$ km, corresponding to the highest possible value for $N$ (see text). $R_i$ is the radius of the inner sphere, and $h$ is the ocean’s thickness. Panel **b** shows that due to geostrophic constrain, the flow exhibits a more complex structure close to the ice shell.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2** | Galileo magnetometer data in the EPhiO coordinates (Europa-centred with x along the direction of corotation, y radially inward towards Jupiter, and z parallel to Jupiter's rotation axis), compared with results from DNS using $h=147$ km. The black curves show time evolution of magnetic field components $B_x$ (top) and $B_y$ (bottom) recorded by Galileo during the trajectory E14 (ref. ³⁶), while other solid lines correspond to our DNS with various values of electrical conductivity $\sigma=2 \times 10^{-2}$, $5 \times 10^{-2}$, $1.8 \times 10^{-1}$, 1.8, 3.6 and $18$ S m⁻¹. Best fits are obtained for electrical conductivities between 0.3 and $3$ S m⁻¹.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 3 | Magnitude of the oceanic jet.** a, Predictive scaling law for the magnitude of the oceanic jet. DNS computed for various values of Ek and the ocean thickness *h* collapse on the law $\bar{U}/c \propto N$, with $N = B_0^2/(\mu_0^2 \rho \nu \sigma c^2)$ (solid black line). Estimate of *N* leads to an oceanic jet of a few centimetres per second for Europa (black empty circle). Error bars indicate maximal and minimal bounds on the value of eddy viscosity. b, Corresponding Joule dissipation versus the electrical conductivity of the ocean.</figcaption></figure>

$$ \bar{U} \propto \frac{B_0^2}{\mu_0^2 \rho \nu \sigma \omega R_E} \qquad (3) $$

meaning that the normalized velocity $Q = \bar{U}/c$ depends on one dimensionless number only, $N = B_0^2/(\mu_0^2 \rho \nu \sigma c^2)$. This formula predicts a very strong equatorial jet of a few kilometres per hour when applied to Europa's parameters. In contrast, one may rather expect a fully inertial regime (ignoring the viscous term), in which a significant part of the injected power per area $P_{ohm} \sim 2B_0^2/\mu_0^3\sigma^2ch^2$ is evacuated through turbulent dissipation. In this case, boundary-layer theory predicts $\gamma P_{ohm} \sim \rho C_D U^3$, leading to a different scaling for the jet velocity $\bar{U} \sim (2\gamma B_0^2/\rho C_D \mu_0^3 \sigma^2 ch^2)^{1/3}$, where $C_D$ is the drag coefficient and $\gamma$ is the ratio between viscous and ohmic dissipation. This prediction rather leads to jet velocities of a few millimetres per second. Finally, an intermediate approach is to assume an eddy viscosity $\nu_t = \nu + C_D \bar{U} h$ in the laminar formula (3). Figure 3a shows a fairly good rescaling of all simulations when $Q$ is plotted as a function of $N$ using such an eddy viscosity. This suggests that molecular viscosity probably provides unrealistically high velocities. In contrast, because the global rotation produces a weakening of the turbulence intensity due to the two-dimensionalization of the large-scale flow$^{21}$, an eddy viscosity much smaller than the one observed in non-rotating turbulence is expected. In consequence, the viscous and turbulent scalings discussed above might be regarded as maximal and minimal bounds on the jet's magnitude. Under the assumption that Jovian moons lie between these two regimes, our simulations therefore predict that Europa have the most powerful jet, with a time-averaged azimuthal flow possibly reaching a few centimetres per second, while Ganymede should exhibit a few millimetres per second and Callisto a nearly negligible jet ($\bar{U} < 1 \text{ mm s}^{-1}$).

Because tidal and radiogenic heating generate hydrothermal plumes with velocities estimated at several centimetres per second, the electromagnetically driven jet described here could very well be negligible compared with thermally driven flow$^{13,14}$, especially because simple estimates of the Lorentz force gives very small values ($F_B \sim 10^{-13} \text{ N m}^{-3}$). To address this question, we now report in Fig. 4a a simulation including thermal buoyancy, in which both N and the convective Rossby number are such that the magnitudes of the magnetically driven jet and the thermally driven flows are similar to what presumably occurs in Europa's ocean. As expected$^{14}$,

geostrophic thermal plumes strongly dominate at small scale, with velocities around $10 \text{ cm s}^{-1}$ and typical diameters of 30 km. As the convective Rossby number $R_0 \sim 3 \times 10^{-2}$ is too small to generate a significant zonal wind, the same simulation with no magnetic forcing (not shown here) displays no large-scale component of the azimuthal velocity field at the equator. Europa's zonal flow shown in Fig. 4 is therefore entirely due to Jupiter's field. Because this magnetically driven jet is the main contribution to the axisymmetric time-averaged azimuthal velocity field ($\bar{U}_\theta \sim 2 \text{ cm s}^{-1}$), it is not suppressed by the vigorous buoyancy force, even with such a small Lorentz force. Note that a more complex situation arises if the convective Rossby number is of order 1, due to the generation of a strong thermal wind$^{12}$ (see our discussion in Methods and Supplementary Fig. 1).

We now discuss some consequences of the existence of a magnetically driven jet on Europa. First, even at $1 \text{ mm s}^{-1}$, such a zonal flow might play a central role in the potential development of life, by providing a time-independent westward transport of radiolytically produced oxidants and other biologically useful substances$^{22}$. Second, the inefficiency of the magnetic pumping ($U \ll c$) yields modest ohmic dissipation, expected to be around $10^7-10^8 \text{ W}$ on Europa (for $Rm \sim 8-80$, Fig. 3b). Radiogenic ($10^{11} \text{ W}$) or solid-body tidal heating ($10^{12} \text{ W}$) therefore represent a much larger contribution. Note, however, that skin effect concentrates ohmic heating inside a very thin layer close to the ice in the polar region (Fig. 4, right), such that the ohmic dissipation can locally reach between $10^{-2}$ and $10^{-1} \text{ mW m}^{-2}$. This represents a non-negligible fraction of the obliquity tidal heating $3 \times 10^{-1} \text{ mW m}^{-2}$ estimated in the ocean$^{11}$, which also peaks near Europa's poles where water plumes are preferentially generated$^{23}$.

More important, contrary to geostrophic buoyancy plumes or time-periodic tidally driven flows, the magnetically driven jet provides a constant unidirectional torque on Europa's ice shell. This new azimuthal force may have a direct influence on the reorientation of Europa's ice shell in the long term. It has been proposed that reorientation could take the form of a non-synchronous rotation produced by tidal torques$^{24}$, or true polar wander in the case of an asymmetric tidal heating in the ice shell$^{25}$. Non-synchronous rotation is also regularly invoked to explain some of the geological features observed at the surface of Europa$^{26}$, and some evidences of reorientation of the ice shell have been reported$^{27,28}$. The question of whether or not Europa's ice shell has reoriented in the past involves many complex phenomena
<figure><img src="image_6.png" /><figcaption>**Fig. 4** | Snapshot of the instantaneous azimuthal velocity field at $r=R_e+h/2$ from the seafloor, for $Ek=10^{-5}$, $Pr=12$, $Pm=10^{-3}$, $\Lambda=10^{-3}$ and $Ra=10^9$. This set of parameters corresponds to $N=2 \times 10^{-4}$, a value expected for Europa in the presence of eddy viscosity. Left and right panels show, respectively, $\phi$-averaged velocity field $U_\phi$ and electrical currents $\sqrt{|J|^2}$. Both the magnetically driven zonal flow and thermal plumes are geostrophic.</figcaption></figure>

(mainly based on tidal forces), but our results suggest that the magnetic field of Jupiter, by generating a net westward motion in the ocean, can influence any possible non-synchronous rotation of Europa. Taking into account this new effect may therefore shed a new light on the formation of Europa's global system of lineaments.

Note that a more complete modelling of subsurface oceans would certainly require tidal effects or topological features of the ice to be described, and more realistic Ekman numbers to be reached. In this respect, the upcoming space missions JUICE (Jupiter ICy moons Explorer)²⁹ and Europa Clipper³⁰ should help to constrain future models, for instance, by providing more precise estimates on the phase lag $\phi_i$ of the induced field. However, global DNS of Europa's interior, by clarifying the dynamics generated below the ice, may already be helpful in the design of these future spacecraft missions.

## Methods

**Numerical modelling.** We simulate Europa's interior as an electrically conducting fluid confined in a rotating spherical shell subject to an oscillatory background field $B_0$ generated by the host planet Jupiter. In space, such a field must satisfy $\nabla \times B_0 = 0$, and mainly consist of a dipole field with a stationary axial dipolar component and an oscillatory equatorial component. We use a spherical coordinate system, with $r$ the radial distance from Europa's centre, $\theta$ the polar angle measured from the North pole and $\phi$ the azimuth angle counted positive along the eastward prograde direction. In this system, the dimensionless non-stationary magnetic field of Jupiter $\tilde{B}_0$ writes:

$$ \tilde{B}_r^0 = \frac{2\sin \theta \cos(\omega\tilde{t} - \phi)}{\dot{r}_j^3} - \frac{3 \sin \theta \cos\phi \cos \omega\tilde{t}}{\dot{r}_j^5} + \frac{3\gamma\dot{r}^2}{\dot{r}_j^5} [\sin \omega\tilde{t} - \sin^2\theta \sin \phi \cos(\omega\tilde{t} - \phi)] $$

$$ \tilde{B}_{\theta}^{0} = \frac{2 \cos \theta \cos(\omega\tilde{t} - \phi)}{\dot{r}_{j}^{3}} - \frac{3 \cos \theta \cos \phi \cos \omega\tilde{t}}{\dot{r}_{j}^{5}} - \frac{3\gamma\dot{r}^{2}}{\dot{r}_{j}^{5}} [\sin \theta \cos \theta \sin \phi \cos(\omega\tilde{t} - \phi) + \gamma\dot{r} \cos \theta \cos(\omega\tilde{t} - \phi)] $$

$$ \tilde{B}_{\phi}^{0} = \frac{2\sin(\omega\tilde{t} - \phi)}{\dot{r}_{j}^{3}} + \frac{3\sin\phi \cos\omega\tilde{t}}{\dot{r}_{j}^{5}} + \frac{3\gamma\dot{r}^{2}}{\dot{r}_{j}^{5}} [\sin\theta \cos\omega\tilde{t} - \sin\theta \sin\phi \sin(\omega\tilde{t} - \phi) - \gamma\dot{r} \sin(\omega\tilde{t} - \phi)] $$

where $\tilde{t}$ is the time $\omega$ is the dimensionless rotation rate of the field in the framework of Europa, $\tilde{r}$ is the radius normalized by $R_E$, and the distance from Jupiter $\dot{r}_j$ is given by $\dot{r}_j = \sqrt{1 + \gamma^2\dot{r}^2 + 2\gamma\dot{r}^2 \sin \theta \sin \phi}$. The parameter $\gamma = R_E/R = 2 \times 10^{-3}$ is the ratio between Europa's radius and the distance Jupiter-Europa. Such a small value for $\gamma$ corresponds to an almost constant magnetic field rotating in the equatorial plane, although our simulations take into account the small variation of the field due to the fact that $\gamma$ is not exactly zero. In a few simulations (not reported here), we observed that the retrograde jet is slightly smaller when the constant vertical component of the magnetic field is removed, an effect that will be described elsewhere.

In Europa's conducting ocean, the magnetic field can be written as $\mathbf{B} = \mathbf{B}^0 + \mathbf{b}$, where $\mathbf{b}$ is the induced magnetic field. The oceanic salwater is an electrically conducting fluid (conductivity $\sigma$) with velocity $\mathbf{u}$ satisfying Maxwell equations:

$$ \frac{\partial \mathbf{B}}{\partial t} = - \nabla \times \text{Ek} \quad \mu_0 \mathbf{J} = \nabla \times \mathbf{B} \qquad (4) $$

By combining these equations with the modified Ohm's law $\mathbf{J} = \sigma(\text{Ek} + \mathbf{u} \times \mathbf{B})$, one obtains the induction equation which gives the evolution of the magnetic field and induced currents. When this equation is coupled to the Navier-Stokes equation and the problem is made dimensionless, one finally obtains MHD equations for a rotating conducting Boussinesq fluid:

$$ \frac{1}{\text{Ek}} \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} + 2\mathbf{e}_z \times \mathbf{u} + \nabla \pi \right) = \Delta \mathbf{u} + \frac{\Lambda}{\text{Pm}} (\nabla \times \mathbf{B}) \times \mathbf{B} \qquad (5) $$

$$ \frac{\partial \mathbf{b}}{\partial t} = \frac{\text{Ek}}{\text{Pm}} \Delta \mathbf{b} + \nabla \times (\mathbf{u} \times (\mathbf{B}_0 + \mathbf{b})) - \frac{\partial \mathbf{B}_0}{\partial t} \qquad (6) $$

$$ \nabla \cdot \mathbf{b} = 0, \quad \nabla \cdot \mathbf{u} = 0 \qquad (7) $$

where we introduce the Ekman number (Ek, ratio of viscous to Coriolis forces), the Elsasser number ($\Lambda$, ratio of Lorentz to Coriolis forces) and the magnetic Prandtl number ($\text{Pm}$, ratio between the viscosity and the magnetic diffusivity), defined as

$$ \text{Ek} = \nu / (\Omega_0 R_E^2), \quad \Lambda = B_0^2 / (\mu_0 \rho_0 \eta \Omega_0), \quad \text{Pm} = \nu / \eta \qquad (8) $$

As usual, $\nu, \eta = \mu_0\sigma, \rho_0\mu_0, \sigma$ denote the kinematic viscosity, the magnetic resistivity, the density, the magnetic permeability and the electrical conductivity, all assumed to be constant. $\Omega_0$ denotes the angular rotation of Europa, which serves as reference frame. In the above, we used the radius of the spherical domain of Europa $R_E$ as unit of length, and $1/\Omega_0$ as unit of time. The ocean's thickness is varied from 100 km to 147 km, which are the typical values generally discussed in the literature¹⁰. Both the ocean-mantle and ocean-ice boundaries are assumed to be electrically insulating and isothermal. No-slip boundary conditions for the velocity field are used, and internal heat sources that may be generated by tides or libration are neglected.

Alternatively, one may also define the Hartmann number $Ha = B_0 h \sqrt{\sigma/\rho\nu}$, which compares the Lorentz force to the viscous force, and the magnetic Reynolds number $Rm = \omega R_E h \mu_0 \sigma$, which measures the ratio between induction and diffusion. Supplementary Table 1 shows that our typical Ekman numbers are many orders of magnitude larger than the value expected in Jovian moons, but Europa's Hartmann and magnetic Reynolds numbers are correctly reproduced by our simulations. This means that our DNS correctly describe the electromagnetic processes occurring in the ocean, but probably underestimate the effect of global rotation on the flow structure. This underestimate of the global rotation also forces us to use larger values of the Elsasser number. Similarly, the weak magnetic Prandtl number used in the simulations implies that the small-scale turbulence of the ocean is not fully resolved, although large values of Pm may also reflects effects due to turbulent diffusivities. Finally, our Rayleigh numbers are chosen to select a convective Rossby number $Ro = \sqrt{\text{RaEk}^2/\text{Pr}}$ in the range estimated for Europa¹², where Ra is the Rayleigh number and Pr is the thermal Prandtl number. This ensures that the size and velocities of thermal plumes are correctly reproduced, and that Europa's ratio between the Lorentz and buoyancy forces is correctly reproduced. Interestingly, both the magnetic Ekman number $Ek_m = \text{Ek}/\text{Pm}$, which appears in the induction equation (6) and the modified Elsasser parameter $\Lambda_m = \Lambda/\text{Pm}$ in equation (5) are correctly reproduced in the numerical model. This similarity in the force balance may explain why our model successfully reproduces several aspects of Europa despite strong differences in Ek, Pm and $\Lambda$. In terms of dimensional units, apart from viscous and thermal diffusivities that are orders of magnitude too large compared with Europa's ocean, all other parameters are correctly reproduced by the simulations: ocean's thickness $h$ is either 100 or 147 km, gravitational acceleration is $g = 1.3 \text{ m s}^{-2}$, thermal expansion coefficient $\alpha = 3 \times 10^{-4} \text{K}^{-1}$, the global rotation rate is $\Omega = 2.1 \times 10^{-5} \text{s}^{-1}$, the magnetic field rotation rate is $\omega = 1.6 \times 10^{-4} \text{s}^{-1}$ and the electrical conductivity is $\sigma \sim 1 \text{ S m}^{-1}$.

The numerical benchmarked³¹ solver used to compute MHD equations is the semi-spectral code PaRoDy³². The code uses a poloidal/toroidal expansion and a
pseudo-spectral spherical harmonic expansion, while the radial discretization is based on finite differences on a stretched grid (allowing for a parallelization by radial domain decomposition). The time integration is performed using a Crank-Nicholson scheme for diffusion terms and an Adams-Bashforth scheme for other terms. For the lowest Ekman numbers reported in this paper, we use 288 points in the radial lowest and 305 spherical harmonic modes, corresponding to 460 points in the $\theta$ direction and 1,024 points in the azimuthal direction $\phi$.

**Galileo measurements.** If we assume that the ice crust is much thinner than the ocean depth and if we neglect magnetospheric plasma effects, the magnitude and phase of the induced field measured by the Galileo spacecraft only depends on the size of the conducting material (ocean's depth $h$), the size of the moon $R_E$, the phase speed $c=\omega R_E$ of the Jovian field, and the electrical properties of the ocean. In the literature, the ratio $\omega\mu_0\sigma R_E^2$ between the skin depth and the moon's radius is generally used to quantify the induction within the ocean. In the present case, we rather introduce the phase speed $c=\omega R_E$ of the magnetic field and use the magnetic Reynolds number $Rm=ch\mu_0\sigma$, which is more natural for describing magnetohydrodynamic effects in a channel of thickness $h$. Using this definition, $Rm$ can also be regarded as the ratio between some typical length scale $\sqrt{R_E}h$ and the skin depth over which electrical currents penetrate the ocean.

Figure 2 has been obtained by performing approximately 100 DNS with various values of Ek, $\Lambda$ and $Pm$, corresponding to different values of $Rm$. For each of these simulations, we solved the Laplace equation outside our spherical domain to compute the magnetic field beyond the ocean, simulated Galileo spacecraft trajectories and compared to space mission data. During its mission, Galileo made multiple flybys of Europa on several orbits around Jupiter (labelled E4, E11, E12, E14, E15, E17 and E19) if we restrict to encounters for which magnetometer data were acquired. Magnetic measurements from encounters E11, E12, E15, E17 and E19, exhibit short-scale fluctuations probably generated by disturbances due to the plasma currents from the Alfvén wings. Following previous studies³³, we therefore restricted our analysis to the remaining encounters E14 and E4, which are closer to the equatorial plane, outside of the Jovian current sheet. As the two flybys provide the same conclusion (best fit obtained for $\sigma \sim 1 \text{ S m}^{-1}$), we only show here the E14 flyby. Similarly to previous studies, we use the so-called EPhiO coordinate system. Figure 2 was obtained with simulations performed at a fixed value of the aspect ratio (corresponding to an ocean depth of $h=147$ km). Note that local magnetospheric plasma disturbances, by reducing the time offset between imposed and induced fields, are expected to improve the agreement with observations³⁴.

**Simple reduced model.** The behaviour observed in our simulations can be understood by studying the equatorial region of Europa, where the magnitude of the oceanic jet is significant. In this region, Jupiter's magnetic field is well approximated by an homogeneous horizontal magnetic field rotating in the $\phi$ direction at constant speed. For simplicity, let us assume that both induced magnetic field and electrical currents can be regarded as travelling waves propagating in the azimuthal direction, such that the total magnetic field writes:

$$ \mathbf{B} = (B_r(r)\mathbf{e}_r + B_\phi(r)\mathbf{e}_\phi)e^{i(\phi-\omega t)} \quad (9) $$

which includes both Jovian and induced magnetic fields. As we focus on the equatorial plane, the $\theta$ dependance of the fields is ignored, and induced currents are assumed to be only along the $\theta$ direction. Similarly, we seek for a simplified velocity field $U=U_\phi(r)\mathbf{e}_\phi$. The ocean is confined in the spherical gap between $r=R_\phi$ and $r=R_\phi+h=R_E$, where $h$ is the thickness of the ocean. In this reduced model, we also ignore the non-time-varying axial component of Jupiter's field, as it would require to seek for more complicated velocity and magnetic fields.

When using the vector potential $A(r)$ such that $\mathbf{B}=\nabla \times (A(r)e^{i(\phi-\omega t)}\mathbf{e}_\phi)$, the induction equation becomes:

$$ -i\omega A + \frac{iA}{r \sin\theta} U(r) = \frac{1}{\mu_0 \sigma} \left[ \frac{1}{r^2} \frac{\partial(r^2 \partial_A)}{\partial r} - \frac{2A}{r^2 \sin^2\theta} \right] \quad (10) $$

At the boundaries $r=R_E$, the magnetic field $B_\phi$ is required to match the applied Jovian magnetic field, $B_\phi(R_\phi, R_E)=B_0$, where $B_0$ represents the magnetic field of Jupiter at Europa's location. This boundary condition is similar to what is generally done in such electromagnetically driven flows¹⁷,³⁴. At the lower boundary $r=R_0$, we use $B_\phi(R_\phi, R_E)=0$, since the time-varying magnetic field wave can not diffuse into the interior due to the large $Rm$.

In the limit of an ocean thickness that is small compared with the moon's radius and by focusing on the jet very close to the equatorial plane ($\theta \sim \pi/2$), the induction equation then becomes:

$$ \partial_{\tau} A = \frac{2}{R^2} \left[ 1 - i \frac{\mu_0 \sigma R}{2} (c - U(r)) \right] A \quad (11) $$

where $R=(R_\phi + R_E)/2$ is the mean radius of the shell and $c=\omega R$ is the phase speed of the travelling magnetic field of Jupiter. This equation is associated with the boundary condition $\partial_\tau A|_{R_E}=B_0$. Similarly, the electrical currents induced in the

channel are given by $j=\frac{1}{\mu_0}\nabla \times \mathbf{B}=\frac{1}{\mu_0}\Delta A$. Using equation (11), the electrical currents read:

$$ j = i \frac{\sigma A}{R} (c - U(r)) \mathbf{e}_{\theta} = \sigma (U(r) - c) B_i \mathbf{e}_{\theta} \quad (12) $$

where $B_i$ is the radial component of the magnetic field. A closed system of equations is obtained by writing the $\phi$-component of the time-averaged Lorentz force:

$$ F = -\frac{1}{2} \mathrm{IRe} \{ j^+ B_r \} = \frac{\sigma}{2R^2} (c - U(r)) |A(r)|^2 \quad (13) $$

where $j^+$ denotes the complex conjugate of $j$. By using this expression of the Lorentz force in the time-averaged Navier-Stokes equation, one finally obtains the additional equation:

$$ \partial_{\tau} U = - \frac{\sigma}{2 \rho \nu R^2} (c - U(r)) |A|^2 \quad (14) $$

**Derivation of the scaling law.** Equations (11) and (14) can be easily solved numerically and compared to our DNS. However, most of the results reported here can be understood without such a numerical integration. Under the so-called block velocity assumption¹⁷,³⁴, we suppose that the velocity $U$ is constant across the gap, except in some boundary layers that we will ignore for now. In this case, equation (11) can be integrated:

$$ [\partial_{\tau} A]_{R_1}^{R_0} = \frac{2}{R^2} \left( 1 - i \frac{\mu_0 \sigma R}{2} (c - U) \right) \int_{R_1}^{R_E} A dr \quad (15) $$

By using boundary conditions for the potential vector, we obtain an expression for the spatially averaged magnetic field $\vec{B}_i = -\int iA/(hR)dr$:

$$ \vec{B}_i = -\frac{1}{h} \int_{R_1}^{R_E} \left( \frac{iA}{R} \right) dr = \frac{-iB_0}{2\chi \left( 1 - i \frac{\mu_0 \sigma R}{2} (c - U) \right)} \quad (16) $$

By combining this averaged expression of $B_i$ with (12) and (13), we obtain a simple expression for the Lorentz force acting on the ocean:

$$ \vec{F} = \frac{\sigma B_0^2 c (1-Q)}{8\chi^2 + 2Rm^2(1-Q)^2} \mathbf{e}_\phi \quad (17) $$

where $\chi=h/R$ is the aspect ratio and $Q=U/c$ is the velocity of the ocean in the midplane normalized by the travelling speed of the Jovian field. We have also introduced the magnetic Reynolds number $Rm=c\mu_0\sigma h$. A simple balance between the Lorentz force and the viscous force therefore suggests:

$$ U \propto \frac{\sigma B_0^2 c (1-Q)}{\rho \nu (8\chi^2 + 2Rm^2(1-Q)^2)} \quad (18) $$

Finally, note that in Jovian moons, it is obvious that the ocean always flows with a typical velocity much smaller than the propagation speed of the Jovian field $c=230 \text{ m s}^{-1}$, independent of the source of motions. It is thus reasonable to take the limit $Q=U/c \ll 1$ in the above expression. In addition, Galileo measurements suggest electrical conductivity around $\sigma \sim 1 \text{ S m}^{-1}$, such that the magnetic Reynolds number is much larger than unity. The previous expression therefore simplifies to:

$$ U \propto \frac{B_0^2}{\mu_0^2 \rho \nu \sigma c} \quad (19) $$

Alternatively, this expression can also be written in dimensionless form $Q \propto N$, where $N=B_0^2/(\mu_0^2 \rho \nu \sigma c^2)$ is a new dimensionless control parameter for the ocean's evolution. As shown in Fig. 3a, all our DNS collapse on this predictive scaling law. Note that the coefficient of proportionality can not be deduced from the reduced model derived here, as it involves details on velocity field fluctuations, geometry of the considered channel and so on. However, the rescaling of our DNS on a single curve provides an estimate of this coefficient, our data suggesting $Q \sim 0.2N$ (solid black line).

Even if a fairly good agreement between the model and the DNS is obtained here, note that the present model ignores several effects, such as the global rotation of the planet or other driving forces. In addition, as discussed in the main part of the article, the huge turbulent fluctuations expected at large Reynolds number may strongly change the dissipation, such that the $\nu$ must be understood as an eddy viscosity. Because the Ekman number is far too large and the velocity fluctuations remain moderate in our simulations, it is relatively difficult to estimate the value of eddy viscosity in Europa's ocean. However, rotating tank experiments³¹ have recently shown that global rotation, by producing a two-dimensionalization of
the flow, strongly reduces friction. Due to the Proudman–Taylor constraint, only the three-dimensional (3D) part of the flow participates in the eddy viscosity $v_1 \sim C_D U_{3D} H$, leading to values significantly larger than the non-rotating case. Although it is speculative, we argue that such a turbulent viscosity might lead to flows of a few centimetres per second or so.

**Comparison with tidally driven and thermally driven flows.** *Tidal effects.* While orbit eccentricity produces flexing of the solid parts of the moon and provides most of the heat that sustains the ocean*, the flow response to tidal forcing can involve either eccentricity or obliquity. In the first case, the tidal dissipation strongly depends on the ocean thickness *h*, and is associated with time-dependent cellular flows with velocities inversely proportional to *h*, leading to small velocities of a few millimetres per second for Europa. When obliquity of the orbit is taken into account, flow velocities of about 9 cm s⁻¹ may be generated, although such large flows require resonant excitation of Rossby waves¹¹.

**Thermal buoyancy.** A radiogenic flux of roughly 4 mW m⁻² from the silicate interior also contributes to the heating of the subsurface ocean, and generates hydrothermal plumes ascending from the seafloor. These plumes merely depend on the ratio between fluid momentum transport and Coriolis forces, measured by the so-called natural Rossby number $Ro = (B^3)^{1/4}/h$, in which *B* is the buoyancy flux emitted by the seafloor source and *f* = 2Ω sin(θ) is the Coriolis parameter. While some authors have suggested a seafloor heat flux around $Q \sim 0.1$ GW (ref. ³⁵), others rather used $Q \sim 1-10$ GW (ref. ⁷), leading to Rossby numbers in the range $Ro \sim 0.01-0.1$. The corresponding velocities associated with hydrothermal plumes are predicted¹³,¹⁴ to lie between 0.9 and 5 cm s⁻¹. In a more recent paper¹², it was suggested that the Rossby number may reach even larger values, $Ro > 1$. In this extreme case, one may expect quasi-3D turbulent convection in the ocean-producing retrograde zonal flows of about 250 cm s⁻¹, comparable to the velocities produced by our MHD mechanism.

In conclusion, the magnetically driven jet, with velocities comprised between 1 and 100 cm s⁻¹, is a non-negligible source of motion in Europa's ocean. Moreover, if the convective Rossby number is smaller than one, this may be the only contribution to the time-averaged axisymmetric zonal flow, which may significantly affect the non-synchronous rotation of the moon.

**Numerical simulations including thermal buoyancy.** The above discussion suggests that the mechanical flow response (of a few millimetres per second) to tidal forces may be neglected, whereas hydrothermal plumes due to radiogenic heating from the mantle remain crucial. Therefore, a full numerical model of Jovian moons should at least combine magnetohydrodynamical effects and thermal buoyancy due to radiogenic heating.

In some of our DNS, previous MHD equations have therefore been coupled to thermal Boussinesq convection with fixed isothermal boundary conditions, introducing the Rayleigh number $Ra = \alpha g \Delta T^3 / (\nu \kappa)$ and the thermal Prandtl number $Pr = \nu / \kappa$, where $\Delta T$ is the unstable superadiabatic gradient imposed through the ocean, $\kappa$ is the thermal diffusivity, $g$ is the gravity and $\alpha$ is the thermal expansion coefficient. We can also define a convective Rossby number $Ro_c = \sqrt{RaEk^2/Pr}$, similar to the natural Rossby number but based on the temperature contrast. Because a complete thermal model of Europa is beyond the scope of the present paper, we only discuss the effect of a thermal gradient imposed between inner and outer spheres and ignore internal heat source as produced by tidal heating.

Typical values for Europa have been estimated by various authors¹,²,¹³, suggesting $Pr \sim 12$ and $Ra \sim 10^{20}-10^{23}$, corresponding to convective Rossby numbers between 0.01 and 1. As usual, numerical limitations do not allow us to reproduce the low viscosity of Europa's ocean. We present results obtained for larger values of Ekman and smaller Rayleigh numbers, $Ek = 10^{-5}$ and Rayleigh numbers ranging from $Ra \sim 10^7$ to $Ra \sim 10^{10}$, such that we reproduce the range of Rossby numbers, $Ro \sim 0.01-1$ relevant to Europa.

Supplementary Fig. 1 shows snapshots of the azimuthal velocity field at mid-height from the seafloor obtained for $Ek = 10^{-5}$, $A = 10^{-1}$ and $Pm = 10^{-3}$ and $Pr = 12$, for 5 different values of the Rayleigh number. For numerical convenience, note that contrary to Fig. 4, this set of parameters corresponds to the upper bound of the magnitude of the jet (viscous scaling), although the same conclusion can be drawn for both regimes. It shows an increase of small-scale convection patterns as the Rayleigh number is increased (from quasi-2D columnar plumes to 3D convection). At small Rossby numbers ($Ro < 1$), thermal convection only weakly modifies the large-scale average structure of the MHD jet. Interestingly, for $Ro > 0.4$, thermal convection also generates a retrograde zonal flow near the equatorial region of the outer sphere, which even reinforces the magnetically driven jet.

As long as the convective Rossby number is smaller than unity, Supplementary Fig. 1f shows that the large-scale flow is always dominated by the magnetically induced retrograde jet. Because the Rossby number for Europa is believed to lie in this limit $Ro < 1$, most of the results reported in the paper are almost unaffected by the presence of thermal buoyancy. Note, however, that our simulations predict that the small-scale part of the flow may be dominated by rotating thermal plumes with typical velocities of a few centimetres per second. Finally, if the convective Rossby number is as large as suggested by Soderlund et al.¹² ($Ro \sim 1$), a retrograde jet of

several hundreds of centimetres per second is generated, nearly three times larger than the magnetically driven jet. Similarly, Supplementary Fig. 2 shows that Joule heating is only weakly modified by the presence of thermal convection.

## Data availability

The data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.

Received: 29 May 2018; Accepted: 30 January 2019;
Published online: 11 March 2019

## References

1. Neubauer, F. M. Oceans inside Jupiter's moons. *Nature* **395**, 749–750 (1998).

2. Khurana, K. K. et al. Induced magnetic fields as evidence for subsurface oceans in Europa and Callisto. *Nature* **395**, 777–780 (1998).

3. Sparks, W. B. et al. Transient water vapor at Europa's south pole. *Science* **343**, 171–174 (2014).

4. Roth, L. et al. Active cryovolcanism on Europa. *Astrophys. J. Lett.* **839**, L18 (2017).

5. Pappalardo, R. T. et al. A Europan ocean? The (circumstantial) geological evidence. In *Proc. Europa Ocean Conf.* 59–60 (San Juan Capistrano Research Institute, 1996).

6. Pappalardo, R. T. et al. Does Europa have a subsurface ocean? Evaluation of the geological evidence. *J. Geophys. Res.* **104**, 24015–24055 (1999).

7. Thomson, R. E. & Delaney, J. R. Evidence for a weakly stratified Europan ocean sustained by seafloor heat flux. *J. Geophys. Res.* **106**, 355–365 (2001).

8. Ross, M. N. & Schubert, G. et al. Tidal heating in an internal ocean model of Europa. *Nature* **325**, 133–134 (1987).

9. Spohn, T. & Schubert, G. Oceans in the icy Galilean satellites of Jupiter. *Icarus* **161**, 456–467 (2003).

10. Bills, B. G. Free and forced obliquities of the Galilean satellites of Jupiter. *Icarus* **175**, 233–245 (2005).

11. Tyler, R. H. Strong ocean tidal flow and heating on moons of the outer planets. *Nature* **456**, 770–773 (2008).

12. Soderlund, K. M., Schmidt, B. E., Wicht, J. & Blankenship, D. D. Ocean-driven heating of Europa's icy shell at low latitudes. *Nat. Geosci.* **7**, 16–19 (2014).

13. Goodman, J. C., Collins, G. C., Marshall, J. & Pierrehumbert, R. T. Hydrothermal plume dynamics on europa: implications for chaos formation. *J. Geophys. Res.* **109**, E03008 (2004).

14. Goodman, J. C. & Lenferink, E. Numerical simulations of marine hydrothermal plumes for Europa and other icy worlds. *Icarus* **221**, 970–983 (2012).

15. Vance, S. & Brown, J. M. Layering and double-diffusion style convection in Europa's ocean. *Icarus* **177**, 506–514 (2005).

16. Colburn, D. S. & Reynolds, R. T. Electrolytic currents in Europa. *Icarus* **63**, 39–44 (1985).

17. Gailitis, A. & Lielausis, O. Instability of homogeneous velocity distribution in an induction-type MHD machine. *Magnetohydrodynamics* **11**, 69–79 (1976).

18. Reddy, K. S., Fauve, S. & Gissinger, C. Instabilities of MHD flows driven by traveling magnetic fields. *Phys. Rev. Fluids* **3**, 063703 (2018).

19. Schilling, N., Neubauer, F. M. & Saur, J. Time-varying interaction of Europa with the Jovian magnetosphere: constraints on the conductivity of Europa's subsurface ocean. *Icarus* **192**, 41–55 (2007).

20. Hand, K. P. & Chyba, C. F. Empirical constraints on the salinity of the Europian ocean and implications for a thin ice shell. *Icarus* **189**, 424–438 (2007).

21. Campagne, A. et al. Turbulent drag in a rotating frame. *J. Fluid. Mech.* **794**, R5 (2016).

22. Greenberg, M. Transport rates of radiolytic substances into Europa's ocean: implications for the potential origin and maintenance of life. *Astrobiology* **10**, 3 (2010).

23. Sparks, W. B. et al. Probing for evidence of plumes on Europa with HST/STIS. *Astrophys. J.* **829**, 121 (2016).

24. Greenberg, M. & Weidenshilling, S. How fast do Galilean satellites spin? *Icarus* **58**, 186–196 (1984).

25. Ojakangas, G. W. & Stevenson, D. J. Polar wander of an ice shell on Europa. *Icarus* **81**, 242–270 (1989).

26. Helfenstein, P. & Parmentier, E. M. Patterns of fracture and tidal stresses due to non-synchronous rotation: implications for fracturing on Europa. *Icarus* **61**, 175–184 (1985).

27. Geissler, P. et al. Evidence for non-synchronous rotation of Europa. *Nature* **391**, 368–370 (1998).

28. Schenk, P., Matsuyama, I. & Nimmo, F. True polar wander on Europa from global-scale small-circle depressions. *Nature* **453**, 368–371 (2008).
29. Grasset, O. et al. JUpiter ICy moons Explorer (JUICE): an ESA mission to orbit Ganymede and to characterise the Jupiter system. *Planet. Space Sci.* **78**, 1–21 (2013).

30. Phillips, C. B. & Pappalardo, R. T. Europa Clipper mission concept: exploring Jupiter's ocean moon. *EOS* **95**, 165–167 (2014).

31. Christensen, U. R. et al. A numerical dynamo benchmark. *Phys. Earth Planet. Inter.* **128**, 25–34 (2001).

32. Dormy, E., Cardin, P. & Jault, D. MHD flow in a slightly differentially rotating spherical shell, with conducting inner core, in a dipolar magnetic field. *Earth Planet. Sci. Lett.* **160**, 15–30 (1998).

33. Zimmer, C., Khurana, K. K. & Kivelson, M. G. Subsurface oceans on Europa and Callisto: constraints from Galileo magnetometer observations. *Icarus* **147**, 329–347 (2000).

34. Gissinger, C., Rodriguez-Imazio, P. & Fauve, S. Instabilities in electromagnetically-driven flows, part I. *Phys. Fluids* **28**, 034101 (2016).

35. Vance, S. & Goodman, J. C. in *Europa* (eds Pappalardo, R. T., McKinnon, W. M. & Khurana, K. K.) 459–482 (Univ. Arizona Press, 2009).

36. Kivelson, M. G. et al. Europa's magnetic signature: Report from Galileo's pass on 19 December 1996. *Science* **276**, 1239–1241 (1997).

## Acknowledgements

The GO-J-MAG-3-RDR-HIGHRES-V1.02 dataset was obtained from the Planetary Data System. This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01)

of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche.

## Author contributions

C.G. conceived the presented idea and developed the theory. L.P. performed the numerical simulations. C.G. and L.P. performed the analysis of the results and contributed to the final manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** is available for this paper at https://doi.org/10.1038/s41550-019-0713-3.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Correspondence and requests for materials** should be addressed to C.G.

**Journal peer review information:** *Nature Astronomy* thanks Jason Goodman and the other anonymous reviewer(s) for their contribution to the peer review of this work.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2019