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

Europa’s subsurface ocean, approximately 100 km deep beneath a thick icy shell, is a prime target in the search for extraterrestrial life. Prior modeling has indicated turbulent ocean dynamics influenced by non-hydrostatic processes and full Coriolis effects. However, earlier studies typically prescribed temperature at the boundaries and neglected salinity and explicit ice–ocean interaction, potentially missing key processes. This study aims to formulate a more self-consistent ocean–ice coupling using prescribed geothermal heat flux at the bottom and a three-equation ice–ocean interface at the top, explicitly including freezing/melting and salinity effects. The goal is to understand how these processes shape Europa’s stratification, circulation, Taylor columns, eddies, convection, jets, and meridional heat transport, and to assess implications for ice thickness patterns that could be tested by future missions.

Europa’s ocean dynamics have been explored with various models and mechanisms, including localized convection plumes possibly linked to surface chaos terrain, tidal forcing leading to Rossby–Haurwitz waves and tidal dissipation, internal waves, and libration-driven instabilities. Due to Europa’s relatively large ocean aspect ratio (~1/16), horizontal Coriolis components are dynamically important, leading to expectations of rotation-axis-aligned convection plumes and alternating zonal jets from scaling arguments. Recent global models using non-hydrostatic dynamics and full Coriolis reported broad low-latitude eastward jets, high-latitude westward jets, and eddy fields, but were adapted from core-convection applications, neglected salinity and freezing/melting, and used prescribed temperature at upper and lower boundaries. Magnetically driven equatorial jets have also been proposed. The present work builds on these by incorporating salinity and a self-consistent surface/bottom flux formulation within an ocean GCM framework at substantially higher resolution.

• Model: MITgcm configured in non-hydrostatic z-spherical coordinates with full Coriolis force (including horizontal components), free surface, and fully nonlinear equation of state relating temperature, salinity, and pressure.
• Configurations: (1) 2D latitude–depth (70°S–70°N, depth 100 km), 1/12° meridional resolution (~2.3 km), 50 vertical levels (2 km); (2) 3D near-global meridional and 30° zonal span with periodic zonal boundaries, 1/24° horizontal resolution, 100 vertical levels (25 m near top to ~1164 m bottom). Time steps: 7200 s (2D), 400 s (3D). Simulations run to statistical steady state.
• Boundary conditions: Bottom boundary with prescribed geothermal heat flux (0.0496 W m⁻², within estimated 5–200 mW m⁻² range); no diffusive salt flux at bottom. Top boundary is an ice–ocean interface represented via the shelf-ice package implementing a three-equation formulation for turbulent heat/salt exchange, latent heat of fusion, and conductive heat flux through ice. Ice thickness is assumed uniform (consistent with results showing efficient poleward heat transport). No-slip at top and bottom boundaries.
• Forcing and parameters: Surface ice temperature from prior studies; mean ocean salinity set to 50 g kg⁻¹ (ppt) as default, with sensitivity tests. Density variations are computed via a nonlinear equation of state; however, linearized approximations are accurate given small T and S excursions.
• Subgrid mixing: Explicit horizontal/vertical eddy viscosity and tracer diffusivity chosen for numerical stability and as conventional in ocean modeling; 2D: ν_h=50 m² s⁻¹, κ_h=5 m² s⁻¹; 3D: ν_h=20 m² s⁻¹, κ_h=2 m² s⁻¹; vertical ν_v=10⁻³ m² s⁻¹, κ_v=10⁻⁴ m² s⁻¹. Gent–McWilliams and KPP parameterizations not used; convection resolved non-hydrostatically.
• Diagnostics and analysis: Identification and characterization of Taylor columns via decomposition of meridional and vertical velocities into components parallel and perpendicular to rotation axis. Spacing analyzed with scaling involving rotation and viscosity. Tracking of meridional propagation of columns. Estimation of resolved eddy diffusivity/viscosity from autocorrelation-based methods and velocity deviations. Energetics: available potential energy (APE) and kinetic energy (KE) estimates following established methods. Meridional heat transport computed including both sensible and latent heat contributions from freezing/melting. Sensitivity studies varying mean salinity (e.g., 5–50 g kg⁻¹), ice thickness (e.g., 5 and 15 km), and bottom heat flux.

• Stratification and density control: The ocean is nearly well mixed with a small top-to-bottom potential temperature difference of ~0.01 °C; stratification is weak and dominated by salinity rather than temperature (βΔS/αΔT ≫ 1). Coldest and saltiest waters occur at low latitudes outside the tangent cylinder due to surface freezing and brine rejection combined with limited along-column heat connection to the bottom.
• Taylor columns: Prominent, transient Taylor columns align with the rotation axis across the ocean. Inside the tangent cylinder (~|lat| ≲ 20°), columns intersect both bottom and ice base and are stationary in latitude; outside, they do not intersect the bottom and propagate equatorward. Observed equatorward propagation outside the tangent cylinder has an indicative rate of ~0.18° per year. Column spacing in 2D is typically <20 km (~0.75° latitude) near the top, increasing with distance from the rotation axis and roughly scaling with sin(φ)/Ω; in 3D, isolated columns have widths/separations of about 20–50 km. Columns are associated with near-zero velocity perpendicular to the rotation axis consistent with Taylor–Proudman constraints.
• Circulation and jets: Zonal jets of a few cm s⁻¹ are present with westward flow in the low-latitude upper ocean and eastward flow elsewhere. Equatorial superrotation is found: in 2D primarily at depth, and in 3D at all depths, attributable to eddy momentum fluxes, Rossby waves, and thermal-wind-related balances. Zonal velocities here are 1–2 orders of magnitude smaller than some prior estimates.
• Convection: Two distinct convective plume types are present—downward plumes near the ice–ocean interface at low latitudes from brine rejection during freezing, and upward plumes from geothermal heating at the bottom. Near the equator, these plumes are nearly perpendicular to Taylor columns and are advected by mean flows. Typical buoyancy/internal wave and convective timescales (2π/|N|) exceed ~50 days.
• Eddies and effective mixing: The flow is richly turbulent with filamentary temperature structures. Based on Taylor column spacing differences and explicit estimates, effective horizontal eddy viscosity/diffusivity is O(200–300 m² s⁻¹) (larger at |lat|<40° and near the bottom), about an order of magnitude smaller than Earth’s tropical ocean values but similar to Snowball Earth estimates.
• Energetics: In the 3D domain, APE ≈ 2.3×10¹⁸ J, with APE/KE ≈ 190 (much smaller than present-day Earth’s >3.3×10⁴), implying macro-turbulence is energized more by convective plumes and barotropic instabilities than by baroclinic instability.
• Meridional heat transport and ice thickness: The modeled meridional heat transport is poleward in both hemispheres, with maxima of ~0.5×10¹¹ W (2D) and ~1.5×10¹¹ W (3D), significantly larger than earlier slab-ocean estimates when latent heat from freezing/melting is included. This transport efficiently redistributes heat to yield an almost uniform ice thickness, consistent with the assumption of uniform ice used in the simulations.
• Sensitivity: Across salinity, bottom heating, and ice thickness variations, key features remain robust: low-latitude cold/salty upper waters outside the tangent cylinder, small vertical temperature contrast (~0.01 °C), westward surface flow with equatorial superrotation at depth (2D), and pervasive Taylor columns. Very low mean salinities can alter density maxima positions due to seawater anomalies but do not eliminate columns.

The incorporation of realistic top and bottom flux boundary conditions and explicit ice–ocean salinity interactions leads to a Europa ocean regime where salinity gradients dominate density and drive unexpected structures: low-latitude cold, saline, dense surface waters and weak overall stratification. The simulated Taylor columns organize transient motions, differ dynamically inside versus outside the tangent cylinder, and—contrary to expectations based on potential vorticity conservation in inviscid settings—propagate meridionally outside the cylinder due to frictional effects. The 3D eddy-resolving simulations reveal that convection from brine rejection and geothermal heating coexists with and is largely orthogonal to low-latitude Taylor columns, indicating distinct mechanisms underpinning columns versus convection. The poleward ocean heat transport, amplified by latent heat effects, is sufficient to maintain nearly uniform ice thickness despite strong meridional gradients in surface ice temperature, supporting the modeling assumption of uniform ice and offering a potentially observable signature for upcoming missions. The relatively low APE/KE ratio compared to Earth suggests that barotropic processes and convective plumes are primary energy pathways for macro-turbulence in Europa’s ocean, with baroclinic instability playing a reduced role.

This work advances understanding of Europa’s ocean by implementing a self-consistent ice–ocean coupling with flux boundary conditions, salinity effects, and high-resolution non-hydrostatic dynamics. Key contributions include: (1) demonstration that salinity dominates stratification and density variations, (2) identification and characterization of transient, rotation-axis-aligned Taylor columns with latitude-dependent spacing and equatorward propagation outside the tangent cylinder, (3) documentation of coexisting convection plumes (from brine rejection and geothermal heating) largely orthogonal to columns at low latitudes, (4) quantification of modest zonal jet speeds and equatorial superrotation, (5) estimation of effective eddy mixing coefficients, and (6) evidence for strong poleward heat transport sufficient to yield nearly uniform ice thickness. These findings produce observationally testable predictions for Europa Clipper and JUICE, including near-uniform ice thickness and potentially detectable salinity structure via magnetic induction. Future work should elucidate the detailed instability mechanisms energizing Taylor columns, explore broader zonal extents and tidal forcing, refine ice dynamics coupling, and improve resolution of surface boundary layers to assess double diffusion.

• The ice shell is treated with a uniform thickness and without explicit ice dynamics (flow), using the shelf-ice package for thermodynamic exchanges only; spatially variable tidal heating in ice and dynamic ice-thickness evolution are not included.
• Ocean tidal dissipation is assumed negligible; explicit tidal forcing and libration-driven processes are not modeled.
• The 3D domain spans only 30° in longitude with periodic boundaries, which may limit representation of very large-scale zonal variability, though eddy scales are much smaller than the domain extent.
• Surface boundary layers are only marginally resolved, limiting definitive conclusions about double-diffusive processes in thin fresh–cold surface layers.
• Subgrid viscosities and diffusivities are held constant (non-KPP/GM), and effective eddy coefficients are diagnosed a posteriori; results may depend on these parameter choices.
• Geothermal heat flux, mean salinity, and ice thickness are uncertain; while sensitivity tests indicate robustness of main conclusions, quantitative details (e.g., column spacing) may vary with parameters and resolution.
• The analysis of Taylor column spacing relies on heuristic scaling arguments involving eddy viscosity; precise mechanisms and scaling require further study.