A scaling law for stagnant-lid convection with composite visco-elasto-plastic rheology: application to Europa's ice shell

Abstract

Scaling laws relating the Nusselt number (Nu) and the Rayleigh number (Ra) are the standard tool for translating numerical models of convection into predictions of surface heat flux and the thicknesses of the stagnant lid and the rheological sublayer. Existing scaling frameworks for icy moons normalise this relationship using a temperature-dependence parameter, γ, calibrated against viscous and incompressible models with composite creep rheology (diffusion, dislocation, grain boundary sliding, and basal slip). These calibrations, however, omit physics included in recent models of Europa’s ice shell: visco-elasto-plastic (VEP) deformation and a free surface, needed to capture realistic surface stress and topography, and compressible (anelastic liquid approximation, ALA) flow, required once shear and adiabatic heating become non-negligible. Whether this additional physics alters the Nu-Ra-γ scaling remains untested, even though heat flux and ice shell thickness estimates for Europa rely on that scaling and would need reassessment if this breaks down under more realistic rheology and boundary conditions. Therefore, we investigate whether the additional physics disturbs this scaling, using numerical simulations of convection in Europa’s ice shell at 60, 40 and 20 km, chosen to bracket independent estimates of Europa’s shell thickness. At each thickness, we ran every combination of four modelling choices: i) rheology (viscoplastic, VP, vs. visco-elasto-plastic, VEP), ii) surface boundary condition (free-slip vs. free-surface), iii) creep mechanism (with or without grain boundary sliding, GBS), and iv) compressibility (i.e., the Boussinesq approximation, BA, vs. the anelastic liquid approximation, ALA). We find that the 20-km cases do not convect within the 100 Myr model runtime. We also note the standardised form of the rheological-sublayer parameter γ_rh matches published expectations only for our diffusion–dislocation cases, and even there departs from the value we measure directly by up to 8.5 (51 per cent). Additionally, we find that γ_rh tracks ice shell thickness rather than cleanly measuring temperature-dependent viscosity for our models. Instead, we find that boundary-layer depth normalises the Nu–Ra relationship far better than either form of the γ_rh, and we report the dimensional depth of the rheological sublayer d_rh, because d_rh keeps the interior Rayleigh number, Ra_i, as an independent, statistically significant predictor of the Nusselt number whereas its nondimensional form d_rh/d_shell does not. Of the four additional physical ingredients, only the presence of GBS creep significantly affects the boundary-layer geometry. Elasticity and the top boundary have no detectable effect, and compressibility only has a small and inconsistently significant one. What does matter is the choice of normalising variable itself: the standardised γ_(rh-std) formulation is only dependent on the creep-mechanism contributions and cannot distinguish our two shell thicknesses, whereas the dimensional depth of the rheological sublayer, d_rh, can. More broadly, these results argue for revisiting how composite-rheology convection models of icy-satellite ice shells are normalised.