Hydrogen Engulfment into Sub-Neptune Cores through Magma Ocean Convective Surface Renewal

Abstract: A convective surface renewal explanation for ingassing of hydrogen from H2-rich envelopes into molten cores of growing planets is presented, with particular application to sub-Neptunes. The ingassing of hydrogen occurs at a moving front comprising the upper boundary of the magma ocean that engulfs hydrogen from the envelope above as it moves. The mechanism of engulfment is surface renewal as the upper diffusive boundary layer is diluted by convection. Surface renewal results in ingassing even where the stabilizing buoyancy of the low-density hydrogen neutralizes the destabilizing thermally-driven buoyancy. It is ultimately driven by cooling. Beginning during the gaseous disk phase, when hydrogen is effectively available in infinite supply, sub-Neptunes can acquire enough hydrogen to permit global core-envelope equilibrium. The structures are therefore consistent with hydrogen-rich, supercritical magma oceans overlain by hydrogen-rich envelopes. The process of surface renewal has implications for other magma ocean - atmosphere interactions in general.

• The paper develops a Danckwerts surface-renewal model showing that convective overturn in a cooling magma ocean can transport envelope hydrogen into sub-Neptune cores without giant impacts.
• The model predicts buoyancy-limited, cooling-controlled ingassing, with a fiducial 6-Earth-mass planet reaching about 3.4% core hydrogen by mass and retaining roughly 50% of its hydrogen there at 1 Gyr.
• The results are largely insensitive to hydrogen diffusivity but depend strongly on accretion heating, with cold-start conditions reducing peak core hydrogen to about 0.3% by mass and leaving major questions about mixing and long-term retention.

Overview and motivation

This paper by Edward D. Young (2608.14518) addresses a specific kinetic gap in sub-Neptune interior models: while thermochemical equilibrium calculations predict that molten silicate cores should dissolve large quantities of hydrogen from overlying H$_2$-rich envelopes, the physical mechanism by which hydrogen actually enters the core has not been established. The paper develops a Danckwerts-type convective surface-renewal model in which the upper boundary of a magma ocean acts as a moving engulfment front that consumes envelope hydrogen as it advances. The central claim is that substantial hydrogen ingassing is an intrinsic consequence of planetary growth—driven solely by top-down cooling of a convecting magma ocean—and requires no exogenous forcing such as giant impacts.

The work builds on recent ab initio molecular dynamics results showing that silicate, iron, and hydrogen are fully miscible above a composition-dependent binodal at sub-Neptune interior conditions (2608.14518). Under supercritical conditions, the melt-envelope interface is defined by bulk composition rather than a saturation limit; below the binodal crest, the melt-side limb of the solvus sets the interfacial concentration. This thermodynamic framework supplies the driving force for ingassing; the present paper supplies the transport kinetics.

Surface renewal versus stagnant-film transport

The paper contrasts three mass-transfer closures. Whitman two-film theory fixes a boundary-layer thickness $\delta$, giving $k = D/\delta \propto D$. Higbie penetration theory replaces the fixed film with a single transient contact of duration $t_c$, yielding $k(t) = \sqrt{D/(\pi t)}$. Danckwerts surface renewal generalizes this to an exponential distribution of surface ages with mean renewal rate $s$, giving $k_{\rm L} = \sqrt{Ds}$. The renewal rate is identified with the convective overturn frequency $s = v_{\rm conv}/d_{\rm ocean}$, using the large-eddy closure with the integral scale set by the ocean depth.

A notable methodological point concerns the choice between large-eddy and small-eddy (Kolmogorov) closures. For nominal magma-ocean parameters ($L \sim 10^6$–$10^7$ m, $\delta$0 m$\delta$1 s$\delta$2, overturn speeds 0.1–1 m s$\delta$3), $\delta$4–$\delta$5, so the small-eddy closure would give transfer coefficients larger by factors of $\delta$6–$\delta$7. The author argues that because delivery from depth to near-surface small eddies occurs in series, the effective coefficient follows a harmonic sum and is controlled by the slower whole-ocean overturn. This identification of the renewal length with the ocean depth is a modeling choice justified physically rather than derived; it effectively caps the mass-transfer coefficient at the large-eddy value.

Buoyancy limitation: a self-regulating flux ceiling

Because dissolved hydrogen lowers melt density substantially—at 6000 K and 3.5 GPa, adding 4% by mass H$\delta$8 to MgSiO$\delta$9 reduces density from 2500 to 1350 kg m$k = D/\delta \propto D$0—ingassing generates a stabilizing compositional buoyancy flux that opposes the thermal buoyancy driving convection. The paper shows that the ingassing flux is bounded:

$k = D/\delta \propto D$1

where $k = D/\delta \propto D$2 is the escaping heat flux. A stationary state ($k = D/\delta \propto D$3) is inconsistent with a top-cooled magma ocean: the convection that would be cancelled by compositional buoyancy is precisely what renews the surface and delivers hydrogen. The self-consistent flux satisfies $k = D/\delta \propto D$4 with $k = D/\delta \propto D$5, and for the strongly driven conditions considered here ($k = D/\delta \propto D$6), the system operates close to the buoyancy ceiling throughout the engulfment window.

Two consequences follow directly. First, since $k = D/\delta \propto D$7, ingassing is cooling-limited rather than diffusivity-limited; a non-cooling ocean ceases to engulf. Second, atmospheric escape thins the envelope, raises $k = D/\delta \propto D$8, and accelerates ingassing. The negative feedbacks also imply that the flux self-regulates just below the cap.

Energy budget and double-diffusive layering

Operating near $k = D/\delta \propto D$9 raises two additional constraints. The mechanical power required to redistribute ingassed hydrogen through the mixed depth, compared against the entrainment-available fraction of convective kinetic energy, gives a dimensionless ratio $t_c$0. Using first-principles MgSiO$t_c$1 viscosities ($t_c$2–4 at 10,000–4000 K) and $t_c$3, the model yields $t_c$4: deep mixing operates marginally at the energy limit, with chemical free energy of dissolution pushing $t_c$5 slightly below unity. The result depends on the assumption $t_c$6; the energy limit would bind if $t_c$7, which the author argues is unlikely given envelope-insulated cooling rates.

Regarding double-diffusive layering, the relevant Lewis number is $t_c$8–8, between the $t_c$9 regime (no persistent layers; progressive erosion of stratification) and $k(t) = \sqrt{D/(\pi t)}$0 regime (transient layers eventually engulfed) found in prior simulations. The conclusion is that any layers are transient and slow but do not eliminate ingassing, though solutions with $k(t) = \sqrt{D/(\pi t)}$1 near unity could contain unresolved transient layering.

Model implementation and numerical results

The planet model integrates a four-component state vector (potential temperature, dissolved interior H$k(t) = \sqrt{D/(\pi t)}$2 mass fraction, envelope mass, engulfment front radius) across disk accretion and post-disk phases, with a Chabrier et al. non-ideal H$k(t) = \sqrt{D/(\pi t)}$3/He equation of state, Freedman opacities, and escape via the larger of XUV photoevaporation and core-powered mass loss. Rotation is included through a reduced vertical convective velocity $k(t) = \sqrt{D/(\pi t)}$4 where applicable.

For the fiducial case—a $k(t) = \sqrt{D/(\pi t)}$5 planet at 7-day orbital period accumulating 10% total hydrogen over a 3 Myr disk lifetime—the peak core hydrogen reaches approximately 3.4% by mass during the disk phase, with envelope base pressures peaking near 8 GPa. Ingassing continues for roughly another Gyr until the engulfment front meets the inward-migrating binodal surface at $k(t) = \sqrt{D/(\pi t)}$6 K and $k(t) = \sqrt{D/(\pi t)}$7 GPa, after which the interior equilibrates along the binodal. At peak (~1 Gyr), core hydrogen comprises ~3% of planetary mass and about 50% of total planetary hydrogen; by 4 Gyr this declines to ~2% and 30% respectively as quasi-static exchange returns hydrogen to the envelope. Contour calculations show peak core hydrogen depends primarily on planet mass and initial bulk hydrogen abundance, though the highest-hydrogen corner of the parameter space produces radii exceeding the observed ~$k(t) = \sqrt{D/(\pi t)}$8 radius cliff and is therefore not applicable to sub-Neptunes—possibly relevant instead to Neptune-like planets with hydrogen-rich magma oceans.

An important robustness result: because the system operates in the buoyancy-limited regime where $k(t) = \sqrt{D/(\pi t)}$9 is independent of $s$0, a ten-fold variation in hydrogen diffusivity shifts core H$s$1 concentrations by only ~1% (e.g., 0.08 absolute in the 3.4% fiducial value). This decouples the conclusions from uncertainties in first-principles diffusivities.

Limitations and open questions

The most consequential assumption is the accretional heating efficiency. With $s$2 (hot start), the interface remains hotter than the binodal during accretion, enabling full miscibility and large engulfed fractions. Halving $s$3 to 0.5 reduces peak core hydrogen by roughly half; a cold start ($s$4) yields only ~0.3% by mass. The mechanism therefore requires that the core-envelope interface exceed the binodal temperature during accretion if large hydrogen inventories are to form—an assumption whose validity depends on poorly constrained accretion-shock physics.

Second, the treatment assumes a well-mixed interior validated by a Péclet number $s$5, a single-reservoir magma ocean, and harmonic-series closure for eddy-scale transport; each is reasonable but not independently verified for these conditions. Third, the longevity of double-diffusive staircases at the upper end of the relevant Prandtl range remains unresolved. Finally, whether hydrogen can remain sequestered ("frozen") in molten cores after complete or near-complete envelope stripping—relevant to stripped-core super-Earths—is explicitly deferred to separate work, as is the question of whether subsequent degassing can lower hydrogen sufficiently to trigger silicate–iron differentiation.

Conclusion

The paper establishes convective surface renewal as a viable, internally driven kinetic route to hydrogen-rich sub-Neptune cores, governed by a Danckwerts transfer coefficient with renewal set by convective overturn and capped by a cooling-prescribed buoyancy ceiling. Because engulfment is most effective during the gaseous disk phase when hydrogen supply is effectively infinite, growing planets can reach global core-envelope equilibrium consistent with supercritical, hydrogen-rich magma oceans overlain by hydrogen-rich envelopes. The engulfment front complements the quasi-static, cooling-driven descent of the silicate–hydrogen binodal modeled for post-formation evolution, implying that self-consistent sub-Neptune evolution requires both stages. The framework extends naturally to other magma ocean–envelope exchanges, including silicate species and noble gases, and bears on observable atmospheres and radii of both sub-Neptunes and stripped super-Earth cores.