The search for extraterrestrial life has focused on Europa, one of Jupiter's Galilean moons, due to its deep (~100 km) subsurface ocean beneath a thick icy shell. Evidence for this ocean includes an induced magnetic field, indications of ice tectonics, and possible water vapor plumes. Previous studies using global models suggested a turbulent Europa ocean, considering non-hydrostatic effects and the Coriolis force. However, these models lacked crucial elements such as consistent boundary conditions representing the interaction between the ocean and icy shell, and the effects of ice melting and freezing on ocean salinity. This study addresses these limitations to provide a more comprehensive understanding of Europa's ocean dynamics and its interaction with the ice shell. The goal is to generate a more realistic model incorporating previously neglected factors and refining predictions about Europa’s ocean, which has implications for the potential habitability of this moon.

Prior research on Europa's ocean dynamics utilized various models and mechanisms. Some proposed localized convection plumes underlying surface patterns, while others emphasized the importance of including horizontal components of the Coriolis force due to the ocean's relatively high aspect ratio (depth over horizontal scale). Scaling arguments predicted alternating zonal jets and Rossby-Haurwitz waves, along with tidal forcing effects on internal waves and libration-driven elliptical instability. Recent global models, while incorporating non-hydrostatic effects and the full Coriolis force, lacked crucial details, notably prescribed temperature boundary conditions rather than accounting for the heat and freshwater fluxes associated with ice-ocean interaction and the influence of salinity variations.

This research employed the MITgcm, a high-resolution ocean General Circulation Model (GCM), in both 2D (latitude-depth) and 3D configurations. The 2D simulations, while computationally less expensive, provided insights into critical elements. The model incorporated all components of the Coriolis force, full non-hydrostatic dynamics, and a prescribed heat flux at the ocean bottom. Crucially, a three-equation formulation was used at the top boundary to represent the interaction between the ocean and icy shell, explicitly considering the effects of melting, freezing, and heat diffusion through the ice. A salinity of 50 ppt (g/kg) was initially chosen, with sensitivity analyses performed to explore the effects of varying salinity. The model assumed a uniform ice shell thickness, a self-consistency of which was later demonstrated. The 3D simulation had a much higher resolution (1/24 of a degree compared to ~1° in previous studies) and spanned 30° longitude with periodic boundary conditions. Explicit eddy diffusion and viscosity coefficients were used to represent unresolved subgrid-scale mixing processes. The 3-equation formulation from the MITgcm shelf-ice package was used to compute the freshwater and heat fluxes across the ice-ocean boundary. This formulation incorporated the effects of temperature, salinity, and pressure on ice freezing and melting and heat diffusion through the ice. The bottom boundary condition specified the geothermal heat flux and a no-diffusive bottom flux for salt. The model ran until statistical steady state was reached in both configurations. Sensitivity analyses were done using varying mean salinity, bottom heating, and ice thickness to assess the robustness of the results.

The 2D simulations revealed that the ocean is weakly stratified, with salinity variations dominating density variations. Surprisingly, the coldest water was found at low latitudes, explained by the dynamics of Taylor columns. These columns, parallel to the rotation axis, efficiently transport heat from the bottom to the surface within the tangent cylinder (radius of Europa's core), but not outside, leading to freezing and higher salinity at low latitudes. The spacing between the Taylor columns was less than 20 km and explained through scaling arguments relating rotation rate and viscosity. The columns exhibited equatorward propagation outside the tangent cylinder. The 3D simulations confirmed the presence of Taylor columns, but also revealed a rich turbulent eddy flow and convection plumes driven by brine rejection during freezing and geothermal heating. The convection plumes were often perpendicular to the Taylor columns at low latitudes. The 3D simulations showed a zonal velocity of several cm s⁻¹, significantly lower than previously reported. Superrotation at the equator was attributed to eddy fluxes of zonal momentum and thermal-wind balance. The 3D simulations strengthened stratification by converting potential energy into kinetic energy. The meridional heat transport calculated was much larger than previous estimates, largely due to the inclusion of latent heat from freezing and melting. The ratio of available potential energy (APE) to kinetic energy (KE) was significantly smaller than on Earth, similar to that of the Snowball Earth ocean. A uniform ice thickness was predicted due to efficient meridional ocean heat transport, potentially observable in future missions.

The findings challenge previous assumptions about Europa's ocean. The dominance of salinity gradients in determining density stratification, the structure and propagation of Taylor columns, and the interplay between Taylor columns and convection plumes provide a more nuanced understanding of ocean dynamics. The lower zonal velocities and the smaller APE/KE ratio suggest a different energy balance compared to Earth’s ocean. The efficient meridional heat transport explains the predicted uniform ice thickness, a prediction that could be tested by future missions. The model’s predictions regarding salinity variations and the overall eddy diffusivity and viscosity could inform future estimates of tidal heating and the role of internal waves.

This study presents a more comprehensive and realistic model of Europa's ocean dynamics, highlighting the crucial role of salinity in shaping the ocean's structure and circulation. The findings suggest a complex interaction of Taylor columns, eddies, and convection, resulting in a weakly stratified ocean with nearly uniform ice thickness. Future missions to Europa can test the model's predictions regarding ice thickness and salinity variations. Further research could focus on refining the model by incorporating more detailed representations of tidal forcing and ice-shell dynamics.

The model used a simplified representation of the icy shell, assuming uniform thickness, and did not include detailed ice flow dynamics. The 3D simulation's zonal extent (30°) might limit the full representation of large-scale eddy dynamics. While sensitivity analyses were performed, uncertainties remain regarding the exact values of parameters like geothermal heat flux and mean ocean salinity. The model's resolution, while high compared to previous studies, might still be insufficient to fully capture all scales of turbulent motion.