Dynamic Europa ocean shows transient Taylor columns and convection driven by ice melting and salinity

The deep (~100 km) ocean of Europa, Jupiter’s moon, covered by a thick icy shell, is one of the most probable places in the solar system to find extraterrestrial life. Yet, its ocean dynamics and its interaction with the ice cover have received little attention. Previous studies suggested that Europa’s ocean is turbulent using a global model and taking into account non-hydrostatic effects and the full Coriolis force. Here we add critical elements, including consistent top and bottom heating boundary conditions and the effects of icy shell melting and freezing on ocean salinity. We find weak stratification that is dominated by salinity variations. The ocean exhibits strong transient convection, eddies, and zonal jets. Transient motions organize in Taylor columns parallel to Europa’s axis of rotation, are static inside of the tangent cylinder and propagate equatorward outside the cylinder. The meridional oceanic heat transport is intense enough to result in a nearly uniform ice thickness, that is expected to be observable in future missions.

Europa's ocean dynamics have been studied using a variety of models and mechanisms 1,[16][17][18][19][21][22][23][24] . It has been suggested that localized ocean convection plumes may underlie the observed surface patterns of Europa 16,17,22 . On Earth, due to the very low oceanic aspect ratio (depth over horizontal scale, ∼ 10 −3 ), only the vertical component of the Coriolis force is relevant. However, the aspect ratio of Europa's ocean is much higher (∼ 1/16), and thus the horizontal components of the Coriolis force must be included and have been suggested to result in convection plumes that are parallel to the axis of rotation 1,21,23,25 . Scaling arguments were used to suggest the existence of alternating zonal jets 1 , and tidal forcing was proposed to lead to Rossby-Haurwitz waves and thus to oceanic tidal dissipation 19 . Tides can also excite internal waves 26 and libration-driven elliptical instability can also drives ocean motions 27 . A recent study of Europa's ocean 23,25 used a global model, taking into account elements such as non-hydrostatic effects and the full Coriolis force, to study the ocean dynamics, and reported a wide low-latitude eastward jet, a high-latitude westward jet, and a rich eddy field. However, the model was adopted from core convection applications and therefore neglected salinity and ice freezing and melting effects that are shown below to dominate those of temperature; it also used upper and lower boundary conditions of prescribed temperature.
Here we show that a more self-consistent formulation, of prescribed bottom heat flux, and a top boundary condition that represents the full interaction with the icy shell and the resulting heat and fresh water fluxes, lead to a very different ocean temperature distribution. Our resolution is higher than that used previously by an order of magnitude, and the viscosity accordingly lower, allowing interesting small scale features to appear.

Results
The model. We use a very high-resolution ocean General Circulation Model (GCM), the MITgcm 28,29 to investigate the ocean dynamics of Europa, first in a 2d (latitude-depth) configuration, and then in a near pole-to-pole 3d geometry. While the 2d simulations lack several important physical processes, these simulations provide invaluable insight into several critical elements that cannot be addressed in 3d, mostly due to computational cost. We include all components of the Coriolis force, and use the full, non-hydrostatic dynamics. We use a prescribed heat flux as a boundary condition at the bottom rather than prescribing the temperature. This allows the temperature, and in particular the vertical temperature gradient (stratification) to be determined by the model. We follow the modern oceanographic literature and use a three equation formulation 30 (Methods, subsection 3-equation top boundary condition formulation) of the interaction between the icy shell and the ocean temperature and salinity fields, which takes into account the effects of freezing and melting of the icy shell, and diffusion of heat through the ice, on the temperature and salinity. The icy shell is assumed of uniform thickness, an assumption that we show below to be self-consistent with the calculated ocean meridional heat fluxes that were shown previously 31 to lead to a uniform ice thickness. Estimates of the mean salinity of Europa's ocean vary widely 32 , and we choose a value that is close to the lower end of estimates, of 50 ppt (g/kg). We later analyze the sensitivity to this choice.
2d model results: Stratification, salinity and Taylor columns. The 2d (latitude-depth) simulations ( Fig. 1) show that the bottom geothermal heating results in a (potential) temperature at depth that is higher than near the ice-ocean interface by a mere 0.01 • C (Fig. 1a), suggesting that the ocean is well mixed. In addition, note several interesting features. First, surprisingly, the coldest water is at low latitudes, despite the much warmer low-latitude ice surface temperature 33,34 , in contradiction to the findings of previous studies of Europa's ocean 23,25 . This is explained below as an effect of the Taylor columns discussed there. The ocean is stably stratified at high latitudes and unstably at low latitudes (Fig. 1c), as opposed to the globally unstable stratification imposed in the above mentioned previous studies. Water density variations are dominated by salinity variations, which dwarf the effects of temperature variations (β∆S/α∆T 1, where α and β are the temperature and salinity expansion coefficients).
Previous studies suggest that the salinity may in fact be even higher than assumed here 32 . In that case, the salinity gradients due to melting and freezing are expected to be even larger, as salinity rate of change is proportional to the fresh water forcing times the mean salinity (e.g., in the limit of a fresh ocean evaporation does not lead to salinity changes). We therefore focus on the sensitivity of our results to lower mean salinity values and show below that for a wide range of parameters the idea that salinity dominates density variations is robust (Supplementary Figs. 1-3).
The zonal velocity (Fig. 2a,d) is westward in the low-latitude upper ocean and eastward elsewhere, with a typical velocity of a few cm per second. The deep equatorial zonal velocity is positive (eastward), indicating a superrotation, further discussed below.
Prominent arc-like structures appear in all fields (Figs. 1,2), which reflect features parallel to the rotation axis when plotted in spherical geometry (Fig. 2e). These are Taylor columns with ocean velocity nearly independent of the direction parallel to the rotation axis, and expected for an ocean with a nearly uniform density. While these columns were previously anticipated based on scaling arguments 1 , simulated and attributed to convection 25 , and seen in simulations of magnetically-driven ocean circulation 24 , their detailed dynamics, structure, role in setting the large-scale temperature and salinity structure, and their spacing and propagation have not been studied.
The velocity along the Taylor columns fluctuates as one moves from Europa's center outward, between being positive and negative. Accordingly, the heat advection changes sign as well. In the region inside of the tangent cylinder that is aligned with the rotation axis and has the radius of Europa's rocky core, corresponding to latitudes less than about 20 •21 , the columns intersect the ocean bottom and the ice base, and their along-column motions effectively transfer the bottom geothermal heat to the ocean surface. However, within the tangent cylinder, there is no such effective bottom-to-surface heat transport mechanism, as the Taylor columns do not intersect Europa's ocean bottom there, and this results in the colder ocean surface in the equatorial regime seen in Fig. 1a. This leads to freezing there, and thus to brine rejection and to the higher salinity as seen in Fig. 1b.
We find that the spacing between the Taylor columns is less than 20 km (0.75 degree latitude, Fig. 3a). In order to analyze the Taylor columns, we project the model's meridional v and vertical w velocity components on the directions parallel and perpendicular to the axis of rotation. The velocity parallel to the axis of rotation, u par = w sin φ + v cos φ where φ is the latitude, shows clear Taylor columns in which it is independent of the direction parallel to the axis of rotation ( Fig. 2b,e), in accordance with the Taylor-Proudman theorem. This parallel velocity is significantly smaller than the zonal velocity (Fig. 2a,d) and significantly larger than the velocity perpendicular to the axis of rotation in the latitude-depth plane u per = w cos φ − v sin φ (Fig. 2c).
In the zonal momentum budget of the 2d model, the two Coriolis terms dominate the others, so that the momentum balance is 2Ωw cos φ − 2Ωv sin φ ≈ 0, where Ω is Europa's rotation rate. This leads to u per ≈ 0, explaining the observation that u per u par . The parallel and zonal velocities are symmetric with respect to the equator, while the perpendicular is anti-symmetric, vanishing at the equator.
The distance between the Taylor columns can be estimated using scaling arguments (Methods, subsection The spacing between the Taylor columns). In the zonal and meridional dominant momentum balances, the sum of the two dominant Coriolis terms is balanced by parameterized horizontal viscosity, and one can form a length scale from the two relevant parameters, the horizontal viscosity coefficient ν h (m 2 s −1 ) and the rotation rate Ω (s −1 ), to find that the relevant length scale is C ν h sin(φ)/Ω where C ≈ 14 is an empirical constant found from the numerical results (Fig. 3a). The horizontal viscosity coefficient we used (50 m 2 s −1 ) represents parameterized viscosity (Methods, subsection Eddy coefficients, subgrid-scale representation).
Further verification of the scaling for the distance between columns is obtained below in the 3d runs, where the effective resolved eddy viscosity is found to be larger (300 m 2 s −1 , see Methods, subsection Estimating eddy coefficients), and the column distance is indeed larger ( Supplementary Fig. 6). While this particular spacing is likely sensitive to model assumptions, the qualitative dynamical insights obtained should be valid. Additional simulations suggest that the distance between the Taylor columns in the high latitudes scales like the square root of the ocean depth. The existence of Taylor columns in Europa, and the corresponding zonal jet structure was predicted by ref. 1 to be related to the Rhines scale, although our findings regarding the spacing between the columns is different from the Rhines scale scaling predicted in that work. These Columns seem to also appear in one of the simulations of 25 that was characterized by low viscosity, although no detailed analysis was provided.
In absence of dissipation the Taylor columns were predicted to be at a fixed latitude, as the potential vorticity (q = (2Ω + ζ)/h, where h is the column height, which depends on latitude) is preserved 1 . Yet we find the columns to show prominent equatorward propagation outside of the tangent cylinder (Fig. 3b,c) whose mechanism would require further elucidation in future work.
No propagation is visible within the tangent cylinder (Fig. 3b,c). We also note the oscillatory variations along the maximum/minimum lines (Fig. 3b,c). While the discussion in this subsection clearly explains the structure and spacing of the (2d) Taylor columns, the necessarily-over simplified eddy viscosity formulation used may affect the simulation. Below we show, based on 3d simulations, that in fact the eddy parameterized coefficient is much larger than the one used in the 2d results, lending credibility to the 2d results. in particular what is the specific instability mechanism involved, requires further study.
The temperature field is clearly turbulent, showing richly complex eddy filaments (Fig. 4a, see supplementary animations). The above relation between viscosity and Taylor column spacing, together with the fact that the Taylor column spacing is larger in the 3d simulation, suggests that the effective eddy viscosity due to resolved eddies is about 15 times larger than the small explicit parameterized viscosity used in the 3d runs for numerical stability, following common ocean modeling practice. This is consistent with an explicit estimate of the eddy coefficients calculated from the 3d runs ( Supplementary Fig. 12). As a result of the eddies and convection plumes, the Taylor columns are less persistent along the direction parallel to the rotation axis than in the 2d simulations (Figs. 4, 5, and Supplementary Figs. [8][9][10][11]. The existence of waves and eddies in the 3d simulation also affects mean flows. The zonal velocity is typically several cm s −1 (Fig. 5b), 1-2 orders of magnitude smaller than the previously reported velocities 23,25 ; thus, our estimate for Europa's ocean kinetic energy (see below) is several orders magnitude smaller than that of these previous studies. Note in particular differences in the vertical structure of the zonal jets along the equator, between the 2d ( zonal velocity with depth is consistent with the thermal-wind relation (i.e., u z = g 2Ωaρ 0 sin(φ) ρ φ ). In the 3d case, we find an additional term that cannot be neglected (i.e., While the stratification is very weak and the ocean well-mixed ( Fig. 1a-c), as in the 2d case, the extent of unstable water column with heavy water above light water is more limited in the 3d case (compare Fig. 1a-c and Figs. 4c, 5a). This is because the eddies in the 3d simulation strengthen the stratification, by converting potential energy into kinetic energy, as was suggested to be the case for Earth's snowball ocean 35,36 . The characteristic time of convection may be estimated via the buoyancy frequency, indicates statically unstable/stable water column. We find typical buoyancy and convection time and with previous higher resolution regional runs 22 .
Refs. 1, 37 suggested the possibility of double diffusion in Europa's ocean. We find in the low latitudes in the 3d model and in the high latitudes of the 2d model, warm salty water under a surface layer of 1-2 grid points that is cold and fresh, and where the stratification is stable. While this is a scenario that can, in principle, lead to double diffusion and therefore to an enhanced vertical mixing, the surface layer is hardly resolved numerically and our results therefore do not seem to provide definite evidence for or against the idea that double diffusion may play a role in Europa's ocean.
One can get further insight into the eddy field from an energetic point of view. The oceanic available potential energy (APE) is the potential energy that may be converted into kinetic energy (KE). The ratio between the APE and the KE provides a measure of the efficiency of kinetic energy extraction from the stratification, and an indication of the source of eddy kinetic energy.
For present-day Earth, the ratio between the APE and the KE ocean is over 33,000 38  Europa's ocean, convective plumes and barotropic instability of the zonal jets play a more prominent role in the generation of ocean macro-turbulence relative to baroclinic instability.
A previous study of the dynamics of the icy shell 31 showed that an efficient meridional ocean heat transport can lead to a uniform shell thickness. The geothermal heat flux entering the ocean from below is larger than the heat escaping through the ice in the tropics and smaller at high latitudes, due to the meridional ice surface temperature gradient 33,34 . This would lead to melting at low latitudes and freezing at high latitudes, leading to ice thickness gradients balanced by ice flow 31 .
However, an efficient poleward ocean heat transport can carry the excess heat meridionally, and thus overcome the tendency toward meridional ice thickness gradients, and result in almost uniform ice thickness (Fig. 6b). The meridional heat transports of the 2d and 3d ocean simulations are shown by the solid lines in Fig. 6 to be poleward, and have maximum values of about 0.5 × 10 11 W and 1.5 × 10 11 W, correspondingly. These estimates of the meridional heat fluxes in a full ocean model are at least four times larger than the heat transport estimated by 31 , because they include the contribution due to latent heat of freezing at the equator and melting at the poles, not included in previous studies.
The meridional heat flux without the latent heat contribution is shown by the green solid curve.
This heat flux is determined, as explained above, by the geothermal and surface heat fluxes calculated for the assumed uniform thickness ice shell. The ocean has no difficulty transporting this heat flux in a way that is consistent with the uniform ice shell assumption. An ocean without an efficient meridional heat flux mechanism would have been heated in the tropics and cooled in the high latitude, not being able to reach a steady state. We conclude that the efficient ocean heat transport in our simulation is self-consistent with the assumption of a uniform ice thickness, justifying the use of a uniform thickness icy shell in this study. That the 3d meridional heat flux is somewhat larger than the 2d flux is a direct result of the larger latent heat due to freezing in the low latitudes in the 3d case. The difference between the two model configurations is not large, and is within the uncertainty of the internal ocean variability, as estimated for example via the difference between the two hemispheres in the 3d case. Spatial variations in tidal heating within the ice may still cause a range of surface heat fluxes and therefore ice thickness variations 33,40 , if the ice is sufficiently thick (thicker than chosen here based on 41 ) to allow convection.

Discussion
The 2d and 3d high resolution simulations of Europa's ocean analyzed here show a number of surprising results. While both the temperature and salinity are nearly uniform, salinity gradients, not considered previously, dominate the gradients in ocean water density, and the heaviest water hemisphere to the other, and higher latitudes at which they intersect the ocean bottom. The Taylor columns, which were previously expected not to propagate in latitude due to potential vorticity conservation, and not to occupy the low-latitude regime 1 , exhibit meridional propagation due to frictional effects and occupy the entire ocean. Their spacing was explained above in terms of the rotation rate and viscosity. The 3d simulation shows a rich turbulent eddy flow, as well as convective plumes due to freezing and brine rejection near the ice-ocean interface, and due to geothermal heating from below. The convection plumes are perpendicular to the Taylor columns at low latitudes. We found superrotation at the equator, attributed it to eddy fluxes of zonal momentum and thermal wind balance, and attempted to explain the reasons for the difference in its structure between the 2d and 3d simulations. The meridional heat flux deduced here is much larger than previously estimated 31 , due to the contribution of the latent heat of freezing that was not considered in previous studies. The ratio between the APE and the KE is significantly smaller than on present-day Earth, yet similar to that estimated for the Snowball Earth ocean 35,36,42 .
A few recently submitted manuscripts investigate complementary aspects of the role of salinity in icy satellites to those discussed here, although they do not deal with the eddy motions and Taylor columns analyzed here. Ref. 43 59 . Tidal heating of the icy shell is not included explicitly, and tidal heating dissipation in the ocean is believed to be negligible 41 .

3-equation top boundary condition formulation.
We use the 3 equations-formulation of 29,30 to calculate the freshwater and heat fluxes between the ice shell and the ocean. The formulation represents an unresolved boundary layer just under the ice where these exchanges occur.
According to these equations, the heat balance of the boundary layer is, where the c p , c p,I are the specific heat and water and ice, ρ, ρ I are the density of ice and ocean, γ T temperature, salinity and pressure 62 . Finally, our resolution in both 2d and 3d is significantly higher than previously used.
Sensitivity to mean salinity, bottom heating and ice thickness. Estimates of Europa's ocean salinity vary widely, from the ocean being nearly fresh to highly saline 63 . Importantly, the magnetometer on the Europa Clipper may be able to estimate the mean ocean salinity 46,47 . The mean salinity affects the freezing temperature of ice and the ocean dynamics, as density variations are found in this work to be driven mostly by salinity gradients rather than temperature gradients.
We In all the simulations, we find that the coldest water is in the upper ocean, outside the tangent cylinder, as in our default simulation analyzed in the paper itself ( Supplementary Fig. 1-5). The salinity is maximal outside the tangent cylinder (around the equator), except for the freshwater case (1st raw). The density is maximal outside the tangent cylinder for all simulation except the two lowest mean salinities (1st and 2nd rows) for which the anomaly of sea water leads to denser water at the bottom due to the bottom heating there. The top to bottom temperature difference is robustly at around 0.01 • C, even when the bottom heat flux is changed from our default value. In all simulations the flow is westward except the bottom equatorial region for which superrotation is observed as discussed in the paper. The Taylor columns structure is clearly visible in the meridional velocity and is similar in all simulations, consistent with the results presented in the main text. Moreover, the structure of the temperature, salinity, and density fields is similar to structure of the those presented in the main text (Fig. 1). We conclude that the sensitivity simulations indicate the robustness of the results of the default experiment analyzed in the paper.
The dominant ocean momentum balance. In order to explore the momentum balance of Europa's ocean, we use the output of the model that uses the full set of equations as explained above, but consider only those terms that are not negligible. We therefore consider the following equations, assuming zonal symmetry, and neglecting the small curvature terms, as well as vertical viscosity/diffusion terms which are very small due to the weakly stratified nature of Europa's ocean and due to the relatively small (in comparison to the horizontal one) vertical/diffusion viscosity coefficient. The momentum equations are then where φ, z, t indicate the latitude, depth, and time, u, v, w are the zonal, meridional and vertical velocities, p hd,nh is the hydrostatic/nonhydrostatic pressure, ρ is the density, a is the radius of Europa, Ω is the rotation frequency, and g is the gravity acceleration. The continuity equation assuming again zonal symmetry is, The different terms in the momentum and continuity equations are shown in Supplementary   Figs. 14-16.The Coriolis terms dominate the zonal horizontal momentum equation, followed by the horizontal viscosity term. In the meridional momentum equation, the balance is geostrophic: the horizontal Coriolis term balances the pressure gradient term, and the horizontal viscosity term is smaller yet not completely negligible.
While the model simulations shown in this work use a fully nonlinear equation of state 62 relating density to temperature, salinity and pressure, we note that a linearized equation can be written as  Supplementary Fig. 15a,b) and the balance between them yields, This is more easily understood by writing the geostrophic approximation in cylindrical coordinates, where now v r is the velocity perpendicular to the axis of rotation, written in terms of the spherical coordinate velocity field as v r = w cos φ − v sin φ and θ is the longitude. The 2d model configuration assumes zonal symmetry (no variations in θ), so that the last equation implies v r ≈ 0, exactly equivalent to (4) in spherical coordinates.
The simple relation (4), written as w≈v tan(φ), does not depend on any parameters and reflects the symmetry properties of v, w: The symmetry of w is opposite of the symmetry of v since tan(φ) is anti-symmetric about the equator. Finally, the velocity parallel to the axis of rotation can be expressed in terms of the meridional and vertical velocities v, w as v z = w sin φ + v cos φ and since w cos φ ≈ v sin φ, v z = v/ cos φ, implying that v z is symmetric as v.
The spacing between the Taylor columns. Based on Supplementary Figs. 14-16, the dominant terms in the zonal and meridional momentum equations near the top of the ocean, where we find viscosity to be non-negligible (2) are, The term 2Ω cos(φ)w in the first equation is small in the upper 10 km or so of the ocean due to the no-normal flow conditions ( Supplementary Fig. 15b). The viscosity term is generally smaller than the Coriolis term, especially in the interior ( Supplementary Fig. 16a,b) in which geostrophy holds, 2Ω sin(φ)ū = − 1 aρ 0 p φ , whereū denotes the geostrophic zonal velocity. Subtracting the geostrophic balance from the fuller momentum equation (6) and approximating the pressure with its geostrophic value throughout, we find, Assuming that we can neglect the meridional gradient of the geostrophic term is small, we can approximate u φ ≈ (u −ū) φ . This heuristic argument is supported by the smoother interior structure of the zonal velocity seen in Fig. 2a. Eqs. (6) can now be written as, whereũ = u −ū. Using a complex variable α =ũ + iv the above equations can be written in terms of a single differential equation that holds near the top of the ocean where viscosity is or Given that the wave number k 2 φ is a slowly varying function of latitude (relative to the meridional scale of the Taylor columns), it can now be used to estimate the spacing between the Taylor columns. The corresponding wavelength in spherical coordinates is λ = 2π/k φ . The distance of a given point, at a latitude φ, along the ocean surface from the axis of rotation is given by a cos φ.
The distance between two adjacent columns in the direction perpendicular to the axis of rotation, is therefore, Or, more explicitly This approximation reproduces the functional behavior shown in Fig. 3a, although a factor that seems to be about π is missing to be consistent with the numerical fit to the simulated distances.
The above arguments are admittedly heuristic at best, yet suggest that eddy viscosity may indeed be at the heart of the process that sets the Taylor column spacing. 1/3 so that 1/3 . One can see that this cannot apply in our case for two reasons. First, there is a jump by a factor of two in the Taylor column height in our case across the tangent cylinder, but no jump is seen in the Taylor column spacing. Second, the Taylor column height in our Europa simulations first increases with latitude from the equator to the tangent cylinder and then decreases with latitude for higher latitudes. Their scaling would predict that the horizontal scale is therefore not monotonous in latitude (i.e., the spacing between the Taylor columns increases from the equator towards the tangent cylinder and then decreases towards the higher latitudes), in contrast to our numerical findings (Fig. 3) and to our own scaling arguments (given above) of monotonically increasing Taylor spacing from the equator towards the high latitudes.
Estimating eddy coefficients. We estimate an horizontal eddy mixing coefficient, κ h , that effectively represents the effects of ocean macro-turbulence, using time series of the zonal velocity at multiple grid points 65,66 . We use the auto-correlation function, R(τ ), and the variance of the zonal velocity temporal anomaly as follows: Lagrangian zonal velocity, the overbar indicates mean over time while the prime indicates the temporal anomaly around this mean. When using the Eulerian velocity field, it is necessary to multiply κ h by a constant γ which we choose to be γ = 4 36,66 . The estimated eddy parameterized viscosity coefficient is larger or equal to the estimated diffusion coefficient 36,66 .
We find that the estimated diffusion coefficient depends on latitude, where for latitudes larger than   Fig. 15a,b), the most dominant terms in the meridional momentum equation are the Coriolis and the pressure gradient terms (Supplementary Fig. 15c,d) which nearly balance each other. The next dominant terms in the zonal and meridional momentum equations are the horizontal viscosity terms (Supplementary Fig. 16a,b).